A STUDY ON FLUID-STRUCTURE INTERACTION OF SWIMMING BEAM USING IMMERSED BOUNDARY- LATTICE BOLTZMANN METHOD By Md Towhidur Rahman A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering โ€“ Master of Science 2021 ABSTRACT A STUDY ON FLUID-STRUCTURE INTERACTION OF SWIMMING BEAM USING IMMERSED BOUNDARY- LATTICE BOLTZMANN METHOD By Md Towhidur Rahman Fluid-structure interaction (FSI) study is of great importance to understand the hydrodynamic coupling of biological swimmers in surrounding environmental domain. Multiple numerical and experimental studies have taken place to capture the behavioral pattern from the environment, explore the physical phenomena and comprehension of dynamics to make contribution in real life applications. In this study, an immersed boundary-lattice Boltzmann method (IB-LBM) for fluid-structure interaction problems is presented. The impact of solid structure on to the surrounding fluid domain is dealt with by immersed boundary method (IBM), where the structure is assumed to be immersed into surrounding fluid and the effect of the immersed boundary are considered by exertion of Lagrangian force onto the surrounding fluid grid points as body force. The flow dynamics is determined by solving discrete lattice Boltzmann equation of a single relaxation time model. The structural dynamics are solved by the finite difference method. For solving the structural dynamics, inextensibility condition was applied. A staggered grid is used in the Lagrangian coordinate system, where tension force is defined on the interfaces (half-grids) and other variables are defined on the nodes. Tension force is calculated at the intermediate steps and used as inextensibility constraint to obtain filament position at the next time step. In the present study, a detailed derivation and corresponding discretization is done for multiple free-swimming cases for a thin flexible filament. The thin flexible filament is actuated by imposing oscillatory heaving and pitching motion at the leading edge with prescribed control parameters. The flow physics of the system is investigated and pressure on the surfaces of the flexible filament is obtained. The results obtained in this study shows consistency with previous publications. The presented computational modelling may be used in future with multiple obstacles in the domain, to investigate the surface pressure variation of the swimming flexible filament and generated data sets may contribute to optimization of control mechanism of the swimmer. ACKNOWLEDGMENTS My heartfelt gratitude to my primary advisor Dr. Tong Gao for patiently mentoring me and constantly guiding me to achieve progress towards this thesis. Without his ideas and directions this work would not have been possible. His support for developing my theoretical knowledge base, applying mathematical methods, and enhancing programming skill helped me to complete this work. I convey my cordial thanks to Dr. Xiaobo Tan for following this work closely from the beginning as the Primary Investigator in our research collaboration and providing feedback as a subject matter expert to improve our outcome. He consistently encouraged me to be in right track towards my research and investigate plausible outcomes during the period of problem formulations. My sincere thanks to Dr. Farhad Jaberi for mentoring me during my graduate studies. He has been phenomenal for my clear understanding on computational fluid dynamics. His kind gesture and proper instructions have been instrumental for my progress. I convey my hearties felicitations to Dr. Fang-Bao Tian for helping me through modelling of the simulations and debugging with the code. His recent works and benchmarking on IB-LB method has helped directly to shape my research. I am incredibly lucky to have very helpful team members, Dr. Andrew Hess and Tejas Patel, who supported me during various computational analysis and strenuous debugging process. I convey my sincere thanks to Dr. Zhaowu Lin for guiding me through different aspect of simulation issues with my base code and providing me guideline for necessary modifications. I also thank Ibrahim Tahmid and RB Shahadat for their valuable contribution in the subject matters relevant to my research. I cordially thank my collaborators in the Soft Robotics group for their valuable suggestions. I thank all my committee members for their guidance, and my course teachers at Michigan State University for teaching me the skills I needed to complete my work. iii I would like to thank my family and friends for their support. My heartiest thanks to my wife, Saima Alam, for her constant support and valuable suggestions throughout my academic journey. Finally, I would like to dedicate my work to my beloved parents. Without their endless support, sacrifices and blessings, it would not be possible for me to be the person I am today. iv TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................... vi LIST OF FIGURES ........................................................................................................................ vii 1 Introduction ........................................................................................................................1 1.1 Fluid-Structure Interaction ...........................................................................................1 1.2 Immersed Boundary Method .......................................................................................3 1.3 Lattice Boltzmann Method...........................................................................................4 2 Literature Review ................................................................................................................6 3 Methodology .....................................................................................................................10 3.1 Problem formulation .................................................................................................10 3.2 Numerical Approach ..................................................................................................13 3.2.1 The lattice Boltzmann method for Navier-Stokes equations ...................................13 3.2.2 The immersed boundary method ...........................................................................15 3.2.3 Discretization of filament governing equation........................................................15 3.2.4 Time marching scheme ..........................................................................................18 4 Derivation and Discretization of Governing Equations and Boundary Conditions ...............20 4.1 Detailed Discretization of governing equation for the flexible filament ......................20 4.2 Discretization of governing equation and boundary conditions for the self-propelled fish-like swimming filament with variable mass ratio and bending rigidity ............................26 5 Results and Discussion.......................................................................................................39 5.1 A pitching flexible filament with fixed leading edge ...................................................39 5.2 Single swimming filament with heaving motion .........................................................43 5.3 Single swimming filament with heaving and pitching motion .....................................48 5.4 Self-propelled fish-like swimmer with variable mass ratio and bending rigidity .........53 6 Conclusion.........................................................................................................................59 APPENDIX .................................................................................................................................65 REFERENCES ..............................................................................................................................65 v LIST OF TABLES Table 5-1: Control Parameters for only heaving motion applied on the leading edge ................43 Table 5-2: Control Parameters for heaving and pitching motion applied on the leading edge ...48 Table 5-3: Control parameter values used in the simulation for self-propelled fish-like swimmer with variable mass ratio and bending rigidity ............................................................................54 Table 5-4: Coefficients in the exponential decay functions of ฮฒ(s) and ฮณ(s) ................................54 Table 5-5: Comparison between the kinematics of a red nose tetra fish Hemigrammus bleheri [73], fish-like self-propelled swimmer by Dai et al. [72] and self-propelled swimming filament using LB-IBM method ................................................................................................................56 vi LIST OF FIGURES Figure 3-1: D2Q9 lattice model ..................................................................................................13 Figure 5-1: Sketch of a pitching flexible filament in a uniform flow............................................39 Figure 5-2: (a) Drag coefficients, (c) lift coefficients, (e) X-coordinates and (g) Y-coordinates of a pitching filament for ๐พ๐ต = 0.125, and (b) drag coefficients, (d) lift coefficients, (f) X-coordinates and (h) Y-coordinates of the pitching filament for ๐พ๐ต = 0.0625. .............................................40 Figure 5-3: Flapping pattern of a pitching filament in a uniform flow for (a) ๐พ๐ต = 0.125 and (b) ๐พ๐ต = 0.0625. The patterns are plotted for 66 instants at rime interval of 0.125. ................41 Figure 5-4: Vorticity fields at four different instants during a flapping cycle of a 2D pitching filament for ๐พ๐ต = 0.125 (left column) and ๐พ๐ต = 0.125 (right column). The corresponding color bars are shown at the bottom of each column. .........................................................................42 Figure 5-5: Sketch of a model for locomotion of a swimming filament with only heaving motion applied on the leading edge. .....................................................................................................43 Figure 5-6: Self-propelled swimming filament with oscillatory heaving motion applied on the leading edge..............................................................................................................................44 Figure 5-7: Drag Coefficient (๐ถD) and lift coefficient (๐ถL) for a swimming filament with heaving motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). ..............................................................................................................44 Figure 5-8: Y-coordinate of the trailing point of a swimming filament with heaving motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25)...................................................................................................................................45 Figure 5-9: Vorticity at six different instants during a period of 2D swimming filament with only heaving motion applied on the leading edge. The corresponding color bars are shown at the bottom of each column. ............................................................................................................46 Figure 5-10: Surface pressure on the swimming filament, with heaving motion applied on the head, at (a) leading edge, (b) trailing edge, and (c) middle point. ..............................................47 Figure 5-11: Self-propelled swimming filament with heaving and pitching motion applied on the leading edge..............................................................................................................................48 Figure 5-12: Drag Coefficient (๐ถD) and lift coefficient (๐ถL) for a swimming filament with heaving and pitching motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25)......................................................................................49 vii Figure 5-13: Y-coordinate of the trailing point of a swimming filament with heaving and pitching motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). ..............................................................................................................49 Figure 5-14: Vorticity at six different instants during a period of 2D swimming filament with heaving and pitching motion applied on the leading edge. The corresponding color bars are shown at the bottom of each column. .......................................................................................50 Figure 5-15: Surface pressure on the swimming filament, with heaving and pitching motion applied on the head, at (a) leading edge, (b) trailing edge, and (c) middle point. .......................51 Figure 5-16: Schematic for computational model with heaving and pitching imposed on the leading edge..............................................................................................................................53 Figure 5-17: The comparison of swimming patterns between a red nose tetra fish Hemigrammus bleheri [73], fish-like self-propelled swimmer by Dai et al. [72] and self-propelled swimming filament using LB-IBM method ..................................................................................................55 Figure 5-18: (a) The time histories of instantaneous swimming speed of the central node of the filament, and (b) drag and lift coefficients of the swimming filament for control parameter and coefficients listed in Table 5-1 and Table 5-2. ............................................................................56 Figure 5-19: Vorticity fields (left column) and pressure fields (right column) at four different instants during a period of a self-propelled fish-like swimmer. The corresponding color bars are shown at the bottom of each column. .......................................................................................57 Figure 5-20: Surface pressure on the filament at (a) leading edge, (b) trailing edge, and (c) middle point. ........................................................................................................................................58 viii 1 Introduction The locomotion and behavior of swimming filament in uniform flow has gained attention in many industrial and biological systems. Different mathematical modelling has been introduced to determine physical parameters involved in this type of fluid-structure interaction problems. In order to mimic locomotion of fish, insects, sperm or micro-organism in fluid medium, 2D filaments and plates are often used. In real case, fish or micro-organism acts as flexible body and deformation occurs during locomotion. The deformations are usually generated by hydrodynamic forces, bending forces, elastic forces, and inertial forces. Also, this deformation dictates the further locomotive behavior. This flapping motion of fins or wings of different animals, including fish, insects etc. are generated for developing thrust or lift to propel or float in the surrounding fluid medium. A thin beam or filament with different control actuation such as heaving and pitching can be used for a simplified modelling to study animal locomotion. 1.1 Fluid-Structure Interaction Fluid-Structure Interaction is an interdisciplinary subject which is associated with the field of fluid dynamics. It is a Multiphysics coupling between laws describing fluid dynamics and structural mechanics. It is characterized by stable or oscillatory interaction between a dynamic deformable structure and a surrounding or internal flow field. If a flow field interacts with a solid object, the flow exerts stresses and strain on the solid structure and the resulting forces generate deformations. The range of deformation depends on the physical parameters, such as pressure, density and velocity of the flow field and material properties of the solid structure. In case of small deformations of the structure and slow variation in time, the flow behavior is not greatly affected by the deformation. Hence, the focus can be kept on stresses generated on the solid structure. But if the deformations are relatively large, it affects flow behavior. Flow velocity and pressure field will be changed, and the problem needs to be treated as a coupled one. This coupled multi-physics problem needs to be treated bidirectionally, where the velocity and pressure of flow field develop solid structural deformations, and this deformation significantly changes the velocity and pressure of surrounding flow field. 1 Fluid-Structure Interactions (FSI) are very significant consideration while designing engineering systems in a vast range of fields, such as, automobiles, spacecrafts, industrial mixer, engines etc. The interaction between bearings and gears, and lubricant, as tribological machine components, is another example of FSI. It also plays a very important role in cardiovascular mechanics, blood flow modelling and micro-organic locomotion. To estimate wall shear stress of blood vessels, the change in tube size due to change in blood pressure and flow velocity needs to be considered properly. These considerations are crucial in case of analyzing aneurysms, hence it has become a standard procedure to apply computational fluid dynamics for analyzing patient specific models. Fluid-Structure Interaction problems are usually very complex for analytical solving, and so they are analyzed by experiments and numerical simulation. With the gradual development of computational fluid and structural dynamics, numerical solution of FSI problems can be treated with great efficiency. Fluid-Structure Interaction can be classified based on the procedures employed for solution. Two main approached are (a) a monolithic approach in which the governing equations for both flow and solid deformation and displacement are solved simultaneously, with one solver, as one unified system, and (b) a partitioned approach in which the governing equations for fluid and solid are solved separately, with two distinct solvers, as two different but coupled system. A partitioned approach permits using independently tested solvers for fluids and solids. Because of this reason, in most engineering application, this approach is preferred. Within this method, the fluid-solid coupling can be applied through a strongly coupled approach or weakly coupled approach. Due to the instability issues of weekly coupled approach, in most cases the strongly coupled approach is preferable. For the last couple of decades, immersed boundary methods have been largely used for FSI applications. These methods are very useful for complex FSI problems, particularly in which it is quite difficult to approach through complex mesh generation. 2 1.2 Immersed Boundary Method In an immersed boundary method (IBM), the solid structure is assumed to be immersed into surrounding fluid and the forces are interacted between the boundary of fluid and solid. In this method, conforming meshes are not necessary, since transfer of interface force is only required. This methodology is useful for working with boundary conditions at fluid-fluid and fluid-solid interfaces based on meshed that are nonconformal to the shape of interfaces. The main advantage of IBM is its capacity to deal with complicated geometries and large boundary movement. For complicated geometries, the mesh generation is easier. With immersed boundary methods, mesh regeneration and mesh movements can be avoided for fluid-structure interaction and flows with non-stagnant boundaries. IBM can process mesh in simpler way than traditional body-conformal grid approach and hence, the computational efficiency could be higher. The immersed boundary method was first introduced and developed by Peskin to simulate cardiovascular mechanics and associated blood flow [34]. The uniqueness of this method was that the complete simulation took place on a Cartesian grid, which was nonconformal to heart geometry. To apply the impact of immersed boundary on the flow field, an unprecedented procedure was introduced. In fitted boundary methods, grid generation is consisted of two segments. Primarily, a surface grid is generated that discretely represents the geometry of the body. Then a structured or unstructured mesh is built with grid generating algorithm to fill the fluid domain. These grids are called body-fitted grids or body-conformal. In immersed boundary method, the grid system comes with much more simplicity. A cartesian grid covers the complete computational domain, including the immersed body. By means of subdividing cells, local grid refinements can be applied in boundary vicinity. In fitted boundary methods, applying the boundary conditions is very simple. As calculations are performed in grid points, information about the solution at the boundary can be easily extracted in case of necessity. This is generally not the case with immersed boundary method since there is no connection between the structure surface and the grid point coordinates. A force function is included as an extra term in the governing equations to impose 3 boundary condition. This can also be done by making change to the numerical stencils in the vicinity of the boundary. This evolving field of immersed boundary method is mostly concentrated on moving boundary flows and simulation of flow around complex geometries. Recent research in these methods is mostly focused toward improving efficiency and accuracy of the solutions. Applying adaptive grid refinement with immerse boundary methods is getting good attention, particularly for flows with high Reynolds number. Most of applications of these IB methods are currently found in multiphase and biological flows. Besides, these methods have great potential for increased application in fluid-structure interaction, complex turbulent flows, and multi-physics simulations. 1.3 Lattice Boltzmann Method Lattice Boltzmann Method (LBM) is a powerful computational tool for simulating the fluid. It was developed in the 1880s based on the statistical mechanics. The LBM can be classified as a mesoscopic approach, in contrast to macroscopic approaches, such as finite difference method and finite volume method, and microscopic approaches, such as molecular dynamics, in computational fluid dynamics. This approach is a particle-based kinetic method in which the behavior of the fictitious fluid particles is modelled, and physical quantities of flow field are determined. While in conventional CFD techniques, the Navierโ€“Stokes equations are solved directly, in the LBM, without solving the Navierโ€“Stokes equations directly, a fluid density on a lattice is simulated with streaming and collision (relaxation) processes by using the Lattice Boltzmann Equation (LBE). Furthermore, by using the Chapmanโ€“Enskog expansion, the Navierโ€“ Stokes equations can be derived from the Boltzmann equation The unique feature of LBE is how it can reduce the degrees of freedom associated with the velocity. Initially, the idea of simplified Boltzmann equations with discrete speeds were mostly meant to generate simpler model equations for rarefied gas dynamics. Then it evolved as a computational alternative for the numerical solution of the Navier-Stokes. LBE as a computational solver for flow problems, has the distinctive property of space-time locality. Quantitative information moves along the straight lines defined by the particle velocities 4 associated with the lattice, instead of the flow speed defined space-time dependent material lines. LBM can also be recognized as a versatile approach as the model fluid can be made to contain usual fluid behavior like vapor/liquid coexistence. So, fluid systems such as liquid droplets can be straightforwardly simulated with LBM. Besides, this method is highly efficient in complex, coupled flow with heat transfer and chemical reactions. The LBM was designed to run efficiently on massive parallel architectures, which contains a wide range from inexpensive embedded Field Programmable Gate Arrays (FPGAs) and Digital Signal Processor (DSPs) to GPUs, heterogeneous clusters, and supercomputers. It enables complex physics and sophisticated algorithms with much more efficiency that leads to better qualitative understanding. The method is obtained from molecular description of a fluid and can associate physical parameters from the interaction between molecules. LBM is also efficient in simulating fluids in complex environments such as porous media, while conventional CFD methods are hard to work with complex boundary structures. For these significant properties, the Lattice Boltzmann Method is widely used in fields of fluid dynamics and associated area of research including material science and biological applications. Though the popularity of LBM in simulating complex fluid systems is increasing, this efficient novel method has some limitations. High-Mach number flows in aerodynamics are difficult to simulate using LBM. The absence of a consistent thermo-hydrodynamic scheme makes it less convenient. However, LBM methods, with Navierโ€“Stokes based CFD, have been successfully coupled with thermal-specific solutions to provide efficiency in heat transfer simulation. For multiphase flow system, generally the interface thickness is larger and the density ratio across the interface is smaller compared to real fluids. But despite of limitations, fast advancements of this method, along with evolving field of wider applications, have been a proof of its great potential in computational physics and numerical simulations. 5 2 Literature Review Study on interaction between flexible bodies and surrounding fluid has drawn significant interest from biomechanics, soft robotics, and material science and engineering. System that involves fluid structure interaction (FSI) is common in nature which includes but not limited to flapping motions of wings/fins by animals such as fishes, birds, and insects. These flapping motions are used for generating thrust to propel or lift to maintain aloft position. The same principle is used by engineers and researchers to develop autonomous air and underwater vehicles in microscale. With significant motivation from the biological FSI examples, there has been rapid growth in understanding efficient control system, structural design, and maneuvering. This phenomenon of fluid and structural interaction, involving a significant amount of momentum exchange between structure and fluid domain, is challenging to model numerically. This is due to different complex geometries and free movement off boundaries. In these systems, the swimming motions are achieved by both-way forcing on the boundary. The flexible solid structure exerts force on the enveloped fluid domain and the fluid also acts on the structure by pressure gradients and viscous shear stresses. To ensure efficient propulsion, the structure needs to act in specific swimming gates [1]. The employment of certain pattern is more important in case of low Reynolds number, where viscous effect is dominant [2]. Proper exchange of momentum between flexible body and fluid, in specific pattern, can lead to a sustainable oscillatory propulsion system analogous to the flapping of a flag. Study on undulatory fish swimming [3, 5] has been conducted focusing on muscle activation and strain pattern, which was further investigated numerically at moderate Reynolds number [4], where vortex shedding is found to be consistent with experimental flow visualization. Flapping plates are mostly used for mimicking fin movement in fluid through exertion of forced by fish fins and insect wings. Fins and wings are flexible complex geometries [6, 7, 10]. Deformation of fins and wings during flapping motion is generally resulted from hydrodynamic forces, bending forces and initial forces. Mechanical properties such as flexural stiffness of wings have been studied to understand the relation between its flexibility and venation pattern [8, 9, 11]. These mechanical 6 properties help to control shape and stiffness of fins for ensuring proper response to external hydrodynamic forces [12]. To understand the motion of flexible filament in fluid, different theoretical [1, 18], experimental [12- 17] and numerical [19- 28] methods have been used. In the work by Zhang et al. [29], the flexible filament motion in a flowing soap was visualized as a two-dimensional model of flapping flag problem. Stretched- straight state and self-sustained state was observed for a single filament, and four distinct dynamic states were found for two side-by-side filaments. Later, numerical simulations on the filament-fluid flow interaction were carried out. For both single and double side- by- side filaments, simulations were modeled by Peskin and coworkers [31- 33], for comparing experimental data using a modified version of immersed boundary [IB] method. In this method, mass of the filament was taken into consideration, and significant impact of mass in flapping dynamics was found. A distributed Lagrangian multiplier/ fictitious domain method (DLM/FD) was applied by Yu [35] to model fluid/ flexible body interaction. Lin et al. [36] studied fluid- structure interaction of soft robotic swimmer by employing a fictitious domain/ active strain method. In the work of Farnell et al. [37], the filament motion was formulated, considering the flexible filament as an N-tuple pendulum, by Lagrange mechanics. Pressure difference was calculated across the filament and hydrodynamic force exerted on the filament was obtained from that. Huang et al. [27] proposed an improved version of immersed boundary method where the fractional step method and a staggered Cartesian grid system was adopted in the Navier- Stokes solver. In this work, direct numerical method was introduced to obtain locomotion of filament under inextensibility constraint. IB methods has gained much attention for fluid- structure interaction simulation as their grid generation requirement is much simplified. This method can handle complex flow geometries without the requirement for a boundary conforming grid, by applying a momentum forcing to the Navier- Stokes equation to model no slip condition [30]. For moving boundaries, it does not require the formation of grid. Eulerian variables are defined on a fixed Cartesian mesh, and they describe the fluid motion. The Lagrangian variables are defined on the freely moving mesh and they describe the immersed boundary motion. A smoothed approximation of Dirac delta function is employed to relate the Eulerian and Lagrangian variables. This method has also been applied 7 for flow modelling which involved bluff bodies, fluid-fluid interfaces, and sharp structure- fluid interfaces. Versatile approaches including discreet forcing IB method [61] and moving-least squared IB method [62] have been used to predict different FSI problems. In recent years, lattice Boltzmann method (LBM) has obtained popularity for fluid flow simulation [39, 40, 48, 55]. This method offers simple formulation and significant advantage of scalability on system of parallel processing. LBM is based on particle kinetics and can bypass discretization of momentum equation. LBM has been applied to study viscoelastic fluids [24, 43, 44], immiscible fluids [45], incompressible flows [46], and nearly incompressible flows [47]. Finite volume method was used to solve constitutive equations and IBM was applied for FSI modeling by Goyal and Derksen [43]. To study FSI problems, IBM and LBM has been coupled for gaining combined advantages [56- 60]. The standard LBM, that utilized a uniform Cartesian mesh, was combined with IBM for the first time by Feng and Michaelides [63, 64] to deal with rigid moving particles. Afterwards, different versions of the IB-LBM have been developed and applied to model flow simulation with rigid bodies [65, 66, 67]. F.B. Tian and his coworkers [38, 49, 50, 51] applied immersed boundary-lattice Boltzmann method (IB-LBM) for modelling interaction between elastic filament and fluid flow with different arrangements. Zhu [52] applied IB-LBM to study 3D FSI flows over a massless flexible sheet. Zhu et al. [53] discussed the use of 3D lattice Boltzmann model (D3Q19) within IB method to simulate viscous flow past a flexible sheet in a 3D channel. Sun et al. [54] used a hybrid method, which coupled the multi- relaxation- time (MRT)-LB and IB models to study hydrodynamic focusing of particles in microchannels. LBM uses a uniform Cartesian mesh and hence, distribution of grid points becomes an issue resulting in high computational cost with increased Reynolds number. Multiblock LBM was introduced to handle this issue, by overlapping fine-resolution Cartesian mesh with a coarse background [68, 69]. Afterwards, multiblock LBM was coupled with IBM for simulating elastic body- incompressible viscous fluid interaction [70] and IB- LBM was employed on non-uniform Cartesian mesh to use grid points efficiently for 3D incompressible flows [71]. In this work, we present immersed boundary-lattice Boltzmann method for solving fluid- structure interaction problems to obtain data on flow dynamics and corresponding structural 8 surface pressure data. We derived and discretized thoroughly the specific mixed-type boundary conditions for combination of horizontally unconstrained and vertically imposed oscillatory heaving and/or pitching motion on the leading edge of a thin flexible filament. The inextensibility condition of the structure has been studied extensively and applied in the computational work. The modelling and results obtained from this study has the potential to be used in generating database for pressure-based control optimization. 9 3 Methodology 3.1 Problem formulation The model problem of a self-propelled flexible filament was considered (Figure 1). The plate has a partially free leading edge where it is clamped and undergoes a harmonic oscillation in the vertical direction but can move freely horizontally. The filament is assumed to be a thin inextensible beam. The governing equations for the filament, written in Lagrangian form, are ๐œ•2๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ•2 ๐œ•2๐‘ฟ (1) ๐œŒl 2 = (๐‘‡ ) โˆ’ 2 (๐›พ 2 ) + ๐œŒl ๐’ˆ โˆ’ ๐‘ญ, ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ•๐‘ฟ (2) . =1, ๐œ•๐‘  ๐œ•๐‘  where, X = (X(s,t),Y(s,t)) refers to the position vector of the filament, s is the Lagrangian coordinate along the filament, T is the tension coefficient along the filament axis, ๐›พ is the bending rigidity, ๐‘ญ is the hydrodynamic force acted on the filament by the surrounding fluid. The fluid flow is assumed to be laminar and incompressible and is governed by the Navier-Stokes (N-S) equation ๐œ•๐’– (3) ๐œŒ0 ( + ๐’–. โˆ‡๐’–) = โˆ’โˆ‡๐‘ + ๐œ‡โˆ‡2 ๐’– + ๐’‡ , ๐œ•๐‘ก โˆ‡. ๐’– = ๐ŸŽ , (4) where u = (u,v) is the fluid velocity vector, ๐œŒ0 is the fluid density, p is the pressure, ๐œ‡ is the dynamic viscosity, and f is the Eulerian forcing applied to enforce no-slip boundary condition along the immersed boundary. Non-dimensionalized form of Eqs. (1)-(4) can be obtained by using characteristic scales: 10 Reference length of filament ๐ฟ๐‘“ Far-field (reference) velocity ๐‘ˆโˆž Time ๐ฟ๐‘“ /๐‘ˆโˆž 2 Pressure ๐œŒ0 ๐‘ˆโˆž 2 Momentum force ๐œŒ0 ๐‘ˆโˆž /๐ฟ๐‘“ 2 Lagrangian force ๐œŒl ๐‘ˆโˆž /๐ฟ๐‘“ 2 Tension force ๐œŒl ๐‘ˆโˆž Bending rigidity 2 2 ๐œŒl ๐‘ˆโˆž ๐ฟ๐‘“ Reynolds number, Re ๐œŒl ๐‘ˆโˆž ๐ฟ๐‘“ /๐œ‡ 2 Froude number, Fr g๐ฟ๐‘“ /๐‘ˆโˆž Mass ratio, ๐›ฝ ๐œŒ1 /(๐œŒ0 ๐ฟ๐‘“ ) As a result of non-dimensionalizing, Eqs. (1) and (3) gets new form respectively as, ๐œ•2๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ•2 ๐œ•2๐‘ฟ ๐’ˆ (5) ๐›ฝ 2 = (๐‘‡ ) โˆ’ 2 (๐›พ 2 ) + ๐›ฝ๐น๐‘Ÿ โˆ’ ๐‘ญ , ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐‘” ๐œ•๐’– 1 2 (6) + ๐’–. โˆ‡๐’– = โˆ’โˆ‡๐‘ + โˆ‡ ๐’–+๐’‡, ๐œ•๐‘ก ๐‘…๐‘’ A Poisson equation for the tension coefficient T(s,t) can be derived from Eqs. (2) and (5) as ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ 1 ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ•๐‘ฟ ๐œ• (7) โ‹… 2 (๐‘‡ ) = ( โ‹… ) โˆ’ ๐›ฝ โ‹… โˆ’ โ‹… ( ๐‘ญ โˆ’ ๐‘ญ) , ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  2 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐’ƒ ๐œ•2 ๐œ•2๐‘ฟ where, ๐‘ญ๐’ƒ = โˆ’ ๐œ•๐‘  2 (๐›พ ๐œ•๐‘  2 ) is the bending force. The first term on the right-hand side of Eq. (7) is theoretically zero, but this is kept for correction of numerical error in the computation due to inextensibility constraints. In the current study, the thin filament represents a swimming body which is actuated by prescribed heaving and pitching motion at the leading edge. Different vertical displacements and 11 rotational angles was imposed. Meanwhile, the horizontal displacement of the leading edge was not restricted to allow free swimming. The boundary conditions at the leading edge of the self-propelled filament are considered as mixed-type, which accommodates both the vertical forced oscillation, clamped condition and horizontally unconstrained condition: (๐‘Œ)๐‘ =0 = ๐‘ฆ(๐‘ก) , (8) ๐œ•๐‘ฟ (9) ( ) = (1,0) , ๐œ•๐‘  ๐‘ =0 ๐œ•3๐‘‹ (10) ( 3) =0, ๐œ•๐‘  ๐‘ =0 ๐‘‡๐‘ =0 = 0 , (11) where X and Y are respectively horizontal and vertical components of displacement vector X. For this mixed-type boundary condition, the forces (tension, bending and hydrodynamic force) are treated accordingly along vertical and horizontal direction. For simply supported (pinned) case at the leading edge, the boundary conditions: (๐‘ฟ)๐‘ =0 = ๐‘ฟ0 , (12) ๐œ•2๐‘ฟ (13) ( ) = (0,0) , ๐œ•๐‘  2 ๐‘ =1 For the trailing edge of the flexible filament, free end condition was employed: ๐œ•2๐‘ฟ (14) ( 2) = (0,0) , ๐œ•๐‘  ๐‘ =1 ๐œ•3๐‘ฟ (15) ( ) = (0,0) , ๐œ•๐‘  3 ๐‘ =1 ๐‘‡๐‘ =1 = 0 , (16) 12 3.2 Numerical Approach In this section, numerical methods are described for solving governing equation of fluid and solid body. The Navier-Stokes equations are solved by applying LBM and immersed boundary method was used to deal with moving boundary and hydrodynamic interactions. 3.2.1 The lattice Boltzmann method for Navier-Stokes equations In this study, the discreet lattice Boltzmann equation of a single relaxation time model was applied to govern the kinematics of the fluid [39,40,66], 1 ๐‘’๐‘ž (17) ๐‘“๐‘– ( ๐ฑ + ๐ž๐‘– ฮ”๐‘ก, ๐‘ก + ฮ”๐‘ก) โˆ’ ๐‘“๐‘– ( ๐ฑ, ๐‘ก) = โˆ’ [๐‘“๐‘– ( ๐ฑ, ๐‘ก) โˆ’ ๐‘“๐‘– ( ๐ฑ, ๐‘ก)] + ฮ”๐‘ก๐น๐‘– , ๐œ ๐‘’๐‘ž where ๐‘“๐‘– ( ๐ฑ, ๐‘ก) is the distribution function, ๐‘“๐‘– ( ๐ฑ, ๐‘ก) is the equilibrium distribution function, ๐ฑ is the position of fluid particle, ๐ž๐‘– is the velocity along direction i, t is the time, ฮ”๐‘ก is the time step, ๐น๐‘– refers to body force term affecting the distribution function, and ๐œ is non-dimensional relaxation time. The relaxation time can be found from its relationship with dynamic viscosity of the solvent in Navier-Stokes equations by [38,39] ๐œ‡๐‘  ๐Ÿ (18) = (๐œ โˆ’ ) ๐‘๐‘ 2 ฮ”๐‘ก , ๐œŒ 2 where ๐‘๐‘ 2 = 1 / โˆš3 is sound speed, ๐œ‡๐‘  is the solvent dynamic viscosity. Figure 3-1: D2Q9 lattice model 13 Here the two-dimensional nine-speed model, briefly D2Q9 model, is shown in figure 3-1, and the nine discreet velocities of the particle are given by, ๐ž0 = (0,0) , (19) ๐œ‹(๐‘–โˆ’1) ๐œ‹(๐‘–โˆ’1) ฮ”๐‘ฅ (20) ๐ž๐‘– = (๐‘๐‘œ๐‘  , ๐‘ ๐‘–๐‘› ) ฮ”๐‘ก , for i = 1 to 4, 2 2 ๐ž๐‘– = (๐‘๐‘œ๐‘  ๐œ‹(๐‘–โˆ’9/2) , ๐‘ ๐‘–๐‘› ๐œ‹(๐‘–โˆ’9/2) โˆš2 ฮ”๐‘ฅ ) ฮ”๐‘ก , for i = 5 to 8, (21) 2 2 where, ฮ”๐‘ฅ is the lattice spacing. In equation (17), the force term ๐น๐‘– and the equilibrium ๐‘’๐‘ž distribution function ๐‘“๐‘– were calculated using [66] ๐Ÿ ๐ž๐‘– โˆ’ ๐ฎ ๐ž๐‘– . ๐ฎ (22) ๐น๐‘– = (1 โˆ’ ) ๐‘ค๐‘– ( 2 + 4 ๐ž๐‘– ) . ๐Ÿ , 2๐œ ๐‘๐‘  ๐‘๐‘  ๐‘’๐‘ž ๐ž๐‘– . ๐ฎ ๐ฎ๐ฎ โˆถ (๐ž๐‘– ๐ž๐‘– โˆ’ ๐‘๐‘ 2 ๐ˆ) (23) ๐‘“๐‘– = ๐‘ค๐‘– ๐œŒ (1 + + ), ๐‘๐‘ 2 2๐‘๐‘ 4 where u = (u,v) is velocity of the fluid, ๐‘๐‘ 2 = 1 / โˆš3 is sound speed, f is the body force acting on the fluid, and ๐‘ค๐‘– are weight factors. For D2Q9 lattice model, weight factors are given by ๐‘ค0 = 4/9, ๐‘ค1 = ๐‘ค2 = ๐‘ค3 = ๐‘ค4 = 1/9, and ๐‘ค5 = ๐‘ค6 = ๐‘ค7 = ๐‘ค8 = 1/36. The particle density distribution is computed by the series of equations stated, and once it is known, the fluid density, velocity and pressure was calculated by [38,39] ๐œŒ = โˆ‘ ๐‘“๐‘– , (24) ๐‘– ฮฃ๐‘– ๐ž๐‘– ๐‘“๐‘– + 0.5๐Ÿฮ”๐‘ก (25) ๐ฎ= , ๐œŒ P = ๐œŒ๐‘๐‘ 2 . (26) 14 3.2.2 The immersed boundary method Immersed boundary method has been employed to quantify the interaction between fluid and structure. This method considers the boundary effect on fluid by exerting the Lagrangian force on to the fluid as body force. The following expression was used for this method [24,38], ๐’‡(๐ฑ, ๐‘ก) = โˆซ ๐‘ญ(๐’”, ๐‘ก)๐›ฟโ„Ž (๐ฑ โˆ’ ๐‘ฟ(๐‘ , ๐‘ก)) โ…†๐’” , (27) where ๐‘ญ is Lagrangian force density, ๐’‡(๐ฑ, ๐‘ก) is fluid body force density, ๐›ฟโ„Ž (๐ฑ โˆ’ ๐‘ฟ(๐‘ , ๐‘ก)) is a smoothed approximation of Dirac delta function. In two-dimension, the Dirac delta function was calculated by 1 ๐‘ฅ โˆ’ ๐‘‹ (๐’”, ๐‘ก) ๐‘ฆ โˆ’ ๐‘Œ (๐’”, ๐‘ก) (28) ๐›ฟโ„Ž (๐ฑ โˆ’ ๐‘ฟ(๐‘ , ๐‘ก)) = 2 ๐œ™( )๐œ™( ), โ„Ž โ„Ž โ„Ž where h refers to the mesh size. ๐œ™ was defined by using four-point delta function as 1 (3 โˆ’ 2|๐‘Ÿ| + โˆš1 + 4|๐‘Ÿ| โˆ’ 4๐‘Ÿ 2 0 โ‰ค |๐‘Ÿ| < 1 , (29) 8 ๐œ™(r) = {1 (5 โˆ’ 2|๐‘Ÿ| โˆ’ โˆšโˆ’7 + 12|๐‘Ÿ| โˆ’ 4๐‘Ÿ 2 1 โ‰ค |๐‘Ÿ| < 2 , 8 0, 2 โ‰ค |๐‘Ÿ| , The Lagrangian force ๐‘ญ was derived for structure with mass by the following equation ๐” โˆ’ ๐ฎ๐‘–๐‘ (30) ๐‘ญ=๐›ผ , ฮ”๐‘ก ๐ฎ๐‘–๐‘ (๐ฌ, ๐‘ก) = โˆซ ๐ฎ(๐ฑ, t)๐›ฟโ„Ž (๐ฑ โˆ’ ๐‘ฟ(๐‘ , ๐‘ก)) d๐ฑ , (31) ๐œ•๐‘ฟ(๐‘ , ๐‘ก) (32) = ๐ฎ๐‘–๐‘ (๐’”, ๐‘ก) , ๐œ•๐‘ก Where U is structural velocity, which was achieved by solving the dynamics of structure, ๐›ผ is a positive constant, and ๐ฎ๐‘–๐‘ is the velocity of the filament which was obtained through interpolation from the flow field, and filament position is calculated and updated explicitly. 3.2.3 Discretization of filament governing equation In the Lagrangian coordinate system, a staggered grid has been used. The indices start from the leading edge (๐‘– = 0), where a mixed boundary condition is applied, that is a blend of horizontally free and vertically forced oscillation and stops at the trailing end (๐‘– = ๐‘), which is treated as a 15 free end. All the variables are defined on the nodes, except for tension force defined on the interfaces (half nodes). For an arbitrary variable ฮต, first-order central, forward and backward difference approximations are described as follows, ๐ท๐‘ 0 ๐œ– = (๐œ–(๐‘  + ฮ”๐‘ /2) โˆ’ ๐œ– (๐‘  โˆ’ ฮ”๐‘ /2)) โˆ• ฮ”๐‘ ) , (33) ๐ท๐‘ + ๐œ– = (๐œ– (๐‘  + ฮ”๐‘ ) โˆ’ ๐œ– (๐‘ )) โˆ• ฮ”๐‘ ) , (34) ๐ท๐‘ โˆ’ ๐œ– = (๐œ– (๐‘ ) โˆ’ ๐œ– (๐‘  โˆ’ ฮ”๐‘ )) โˆ• ฮ”๐‘ ) , (35) and the second-order central difference approximation is denoted by ๐ท๐‘ + ๐ท๐‘ โˆ’ ๐œ– = (๐œ–(๐‘  + ฮ”๐‘ ) โˆ’ 2๐œ– (๐‘ ) + ๐œ– (๐‘  โˆ’ ฮ”๐‘ )) โˆ• ฮ”๐‘  2 ) , (36) where ๐ท๐‘  and ๐ท๐‘ ๐‘  denotes first and second order derivatives of arclength s. Same approximations apply for discretization along time (t). The bending force term is discretized as following (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–+1 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– + (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 (37) (๐‘ญ๐‘ )๐‘– = โˆ’[๐ท๐‘ + ๐ท๐‘ โˆ’ (๐›พ๐ท๐‘ ๐‘  ๐‘ฟ)]๐‘– = โˆ’๐›พ , ฮ”๐‘  2 for ๐‘– = 1, 2, . . . , ๐‘ โˆ’ 1 and ๐›พ is constant , (1) (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)0 (๐ท๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)0 (38) (๐น๐‘ ) = โˆ’๐›พ = , 0 ฮ”๐‘  ฮ”๐‘  2 (2) (๐ท๐‘ ๐‘ ๐‘  ๐‘Œ)1 โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘Œ)0 (๐ท๐‘ ๐‘  ๐‘Œ)2 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘Œ)1 + (๐ท๐‘ ๐‘  ๐‘Œ)0 (39) (๐น๐‘ ) = โˆ’๐›พ = , 0 ฮ”๐‘  ฮ”๐‘  2 (๐ท ๐‘ฟ) โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘ฟ)๐‘โˆ’1 (๐ท๐‘ ๐‘  ๐‘‹)๐‘ โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)๐‘โˆ’1 (40) (๐‘ญ๐‘ )๐‘ = โˆ’๐›พ ๐‘ ๐‘ ๐‘  ๐‘ = ๐›พ , ฮ”๐‘  ฮ”๐‘  2 where (๐ท๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)0 (41) (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)1 = , ฮ”๐‘  (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– = (๐ท๐‘ + ๐ท๐‘ โˆ’ ๐‘ฟ)๐‘– , ๐‘– = 1, 2, . . . , ๐‘ โˆ’ 1 (42) 16 The boundary conditions at mixed end (๐‘– = 0) and free end (๐‘– = ๐‘) are described as (๐ท๐‘ ๐‘  ๐‘ฟ)0 = (0,0) , for simply supported case (43) (๐ท๐‘ ๐‘  ๐‘ฟ)0 = [(๐ท๐‘  ๐‘ฟ)1/2 โˆ’ (0,0)]/(ฮ”๐‘ /2), for clamped case (44) (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘ = (0,0) , (45) (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)0 = (0,0) , (46) (๐ท๐‘ ๐‘ ๐‘  ๐‘ฟ)๐‘ = (0,0) , (47) In our validation case, we used the discretization of bending force considering ๐›พ as a function of s, and it took an extended form as (๐น๐‘ )๐‘– = โˆ’[๐ท๐‘ + ๐ท๐‘ โˆ’ (๐›พ๐ท๐‘ ๐‘  ๐‘ฟ)]๐‘– (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–+1 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– + (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 = โˆ’๐›พ ฮ”๐‘  2 (๐›พ๐‘–+1 โˆ’ ๐›พ๐‘–โˆ’1 ) ((๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– โˆ’ (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 ) 2 2 โˆ’2[ ] (48) ฮ”๐‘  2 (๐›พ๐‘–+1 โˆ’ 2๐›พ๐‘– + ๐›พ๐‘–โˆ’1 ) (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– โˆ’ ฮ”๐‘  2 for ๐‘– = 1, 2, . . . ๐‘ โˆ’ 1 and ๐›พ is a function of ๐‘  ๐œ• ๐œ•๐‘ฟ The tension force term ๐œ•๐‘  (๐‘‡ ๐œ•๐‘  ) is discretized as ๐‘‡๐‘–+1/2 (๐ท๐‘ 0 ๐‘ฟ)๐‘–+1/2 โˆ’ ๐‘‡๐‘–โˆ’1/2 (๐ท๐‘ 0 ๐‘ฟ)๐‘–โˆ’1/2 (49) [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘ฟ)]๐‘– = [๐ท๐‘ 0 (๐‘‡๐ท๐‘ 0 ๐‘ฟ)]๐‘– = , ฮ”๐‘  for ๐‘– = 1, 2, . . . , ๐‘ โˆ’ 1 , At leading edge and trailing edge, mixed end boundary condition and free end boundary condition are applied respectively 17 (๐‘‡๐ท๐‘  ๐‘‹)1/2 โˆ’ (๐‘‡๐ท๐‘  ๐‘‹)0 ๐‘‡1/2 (๐ท๐‘ 0 ๐‘‹)1/2 (50) [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘‹)]0 = = , ฮ”๐‘ /2 ฮ”๐‘ /2 ๐’ˆ(2) (51) [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘Œ)]0 โˆ’ [๐ท๐‘ ๐‘  (๐›พ๐ท๐‘ ๐‘  ๐‘Œ)]0 = (๐น (2) ) โˆ’ ๐›ฝ๐น๐‘Ÿ ( ) + ๐‘ฆฬˆ (๐‘ก) , 0 ๐‘” (๐‘‡๐ท๐‘  ๐‘ฟ)๐‘โˆ’1/2 โˆ’ (๐‘‡๐ท๐‘  ๐‘ฟ)๐‘ ๐‘‡1/2 (๐ท๐‘ 0 ๐‘ฟ)๐‘โˆ’1/2 (52) [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘ฟ)]๐‘ = = , ฮ”๐‘ /2 ฮ”๐‘ /2 where ๐œ•2๐‘Œ ๐‘ฆฬˆ (๐‘ก) = ( ) ๐œ•๐‘ก 2 0 3.2.4 Time marching scheme Eq. (5) can be written in discretized form as ๐‘ฟ๐‘›+1 ๐‘– โˆ’ 2๐‘ฟ๐‘›๐‘– + ๐‘ฟ๐‘›โˆ’1 ๐‘– ๐’ˆ (53) ๐›ฝ 2 = [๐ท๐‘  (๐‘‡ ๐‘›+1/2 ๐ท๐‘  ๐‘ฟ๐‘›+1 )]๐‘– + (๐‘ญโˆ—๐‘ )๐‘– โˆ’ (๐‘ญ๐‘› )๐‘– + ๐›ฝ๐น๐‘Ÿ ( ) , ฮ”๐‘ก ๐‘” for ๐‘– = 0, 1, 2, . . . , ๐‘ , where ๐‘› represents the ๐‘›th time step, ฮ”๐‘ก represents the increment of time. The bending force term is calculated explicitly as ๐‘ญโˆ—๐‘ = ๐‘ญ๐‘ (2๐‘ฟ๐‘› โˆ’ ๐‘ฟ๐‘›โˆ’1 ) (54) The discretized form of Eq. (7) is expressed as 1 (55) (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘–+ 1 . [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘–+ 2 1 = ๐ท๐‘ก+ ๐ท๐‘กโˆ’ (๐ท๐‘ 0 ๐‘ฟ๐‘› . ๐ท๐‘ 0 ๐‘ฟ๐‘› )๐‘–+1 โˆ’ (๐ท๐‘ 0 ๐‘ผ๐‘› . ๐ท๐‘ 0 ๐‘ผ๐‘› )๐‘–+1 2 2 2 โˆ’ (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . [๐ท๐‘ 0 (๐‘ญโˆ—๐‘ โˆ’ ๐‘ญ๐‘› )]๐‘–+1 , 2 2 for ๐‘– = 0, 1, 2, . . . , ๐‘ โˆ’ 1 , here, the tension force term is calculated implicitly, and ๐‘‡ ๐‘›+1/2 is tension force coefficient at the intermediate time step. ๐‘ฟโˆ— denotes the predicted displacement vector of points in the Lagrangian coordinate system at time step ๐‘› + 1 by applying temporal extrapolation from those points at 18 time steps ๐‘› โˆ’ 1 and ๐‘›, in the form of ๐‘ฟโˆ— = 2๐‘ฟ๐‘› โˆ’ ๐‘ฟ๐‘›โˆ’1 . The first order forward and backward difference operators are represented by ๐ท๐‘ก+ and ๐ท๐‘กโˆ’ . The first term on the right-hand side of Eq. (55) can be discretized as 1 + โˆ’ 0 ๐‘› 0 ๐‘› (๐ท๐‘ 0 ๐‘ฟ . ๐ท๐‘ 0 ๐‘ฟ)๐‘›+1 โˆ’ 2(๐ท๐‘ 0 ๐‘ฟ . ๐ท๐‘ 0 ๐‘ฟ)๐‘› + (๐ท๐‘ 0 ๐‘ฟ . ๐ท๐‘ 0 ๐‘ฟ)๐‘›โˆ’1 ๐ท๐‘ก ๐ท๐‘ก (๐ท๐‘  ๐‘ฟ . ๐ท๐‘  ๐‘ฟ ) = 2 2ฮ”๐‘ก 2 (56) 1 โˆ’ 2(๐ท๐‘ 0 ๐‘ฟ . ๐ท๐‘ 0 ๐‘ฟ)๐‘› + (๐ท๐‘ 0 ๐‘ฟ . ๐ท๐‘ 0 ๐‘ฟ)๐‘›โˆ’1 = , 2ฮ”๐‘ก 2 where the inextensibility condition (๐ท๐‘  ๐‘ฟ . ๐ท๐‘  ๐‘ฟ)๐‘›+1 = 1 has been used. (๐ท๐‘  ๐‘ฟ . ๐ท๐‘  ๐‘ฟ)๐‘› and (๐ท๐‘  ๐‘ฟ . ๐ท๐‘  ๐‘ฟ)๐‘›โˆ’1 terms are penalized as numerical errors introduced in the previous steps. To solve for ๐‘‡ ๐‘›+1/2 , tension force coefficient at intermediate time step, Eq. (55) has been applied, where discretized form of bending and tension force were used, and the boundary conditions were employed for node, ๐‘– = 0 and ๐‘– = ๐‘ โˆ’ 1. After solving ๐‘‡ ๐‘›+1/2 , it was applied as an inextensibility constraint in Eq. (53) to solve for ๐‘ฟ๐‘›+1 . 19 4 Derivation and Discretization of Governing Equations and Boundary Conditions 4.1 Detailed Discretization of governing equation for the flexible filament In this section, the structural governing equations and boundary conditions are discretized. The derivation of Poisson equation for ๐‘‡ is described in Appendix. Eq. (55) is the discretized form of Eq. (7), which is written as 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘–+ 1 . [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘–+ 2 1 = ๐ท๐‘ก+ ๐ท๐‘กโˆ’ (๐ท๐‘ 0 ๐‘ฟ๐‘› . ๐ท๐‘ 0 ๐‘ฟ๐‘› )๐‘–+1 โˆ’ (๐ท๐‘ 0 ๐‘ผ๐‘› . ๐ท๐‘ 0 ๐‘ผ๐‘› )๐‘–+1 2 2 2 โˆ’ (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . [๐ท๐‘ 0 (๐‘ญโˆ—๐‘ โˆ’ ๐‘ญ๐‘› )]๐‘–+1 , 2 2 for ๐‘– = 0, 1, 2, . . . , ๐‘ โˆ’ 1 . 1 The term [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 can be discretized as ๐‘–+ 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘–+1 โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘– [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+1/2 ๐ท๐‘  ๐‘ฟโˆ— ))]๐‘–+1/2 = ฮ”๐‘  1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘–+3 โˆ’ 2๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 +๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–+ 2 ๐‘–โˆ’ 2 2 2 2 = 2 ฮ”๐‘  ๐‘›+ 1 ๐‘›+ 1 ๐‘›+ 1 (57) ๐‘‡ 32 [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] โˆ’ 2๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] + ๐‘‡ 12 [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ] ๐‘–+ ๐‘–+ ๐‘–โˆ’ 2 2 2 = ฮ”๐‘  3 here, ๐‘ซ๐‘  is the first-order central difference operator defined as following, ๐‘ฟ๐‘–+1 โˆ’ ๐‘ฟ๐‘– (๐ท๐‘ 0 ๐‘ฟ) 1 = , ๐‘–+ 2 ฮ”๐‘  ๐‘ฟ๐‘– โˆ’ ๐‘ฟ๐‘–โˆ’1 (๐ท๐‘ 0 ๐‘ฟ) 1 = , (58) ๐‘–โˆ’ 2 ฮ”๐‘  ๐‘ฟ๐‘–+2 โˆ’ ๐‘ฟ๐‘–+1 (๐ท๐‘ 0 ๐‘ฟ) 3 = , ๐‘–+ 2 ฮ”๐‘  20 The RHS of Eq. (57) is derived by applying the following steps, 1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘–+3 โˆ’๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 ๐‘–+ 2 ๐‘–+ 2 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘–+1 = ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] โˆ’ ๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] ๐‘–+ ๐‘–+ 2 2 = ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 โˆ’ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–โˆ’ 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘– = 2 2 (59) ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ 2 โˆ— ๐‘‡ 1 [๐‘ฟ๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] โˆ’ ๐‘‡ 2 1 [ ๐‘ฟ๐‘– โˆ— โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ] ๐‘–+ ๐‘–โˆ’ 2 2 = ฮ”๐‘  2 So, we can express the LHS of Eq. (55) as 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘–+ 1 . [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘–+ 2 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘–+3 โˆ’ 2๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 +๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–+ 2 ๐‘–โˆ’ 2 2 2 2 = (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . 2 ฮ”๐‘  2 1 1 1 ๐‘›+ 2( 0 โˆ—) 0 โˆ—) ๐‘›+ 2 0 โˆ— 0 โˆ— = (๐‘‡ ๐ท ๐‘ฟ 1 . ( ๐ท ๐‘ฟ 3 โˆ’ 2๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 . (๐ท๐‘  ๐‘ฟ )๐‘–+1 ฮ”๐‘  2 ๐‘–+32 ๐‘  ๐‘–+ 2 ๐‘  ๐‘–+ 2 ๐‘–+ 2 2 2 ๐‘›+ 1 (60) 2 0 โˆ— +๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 . (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–โˆ’1 ) ๐‘–โˆ’ 2 2 2 1 1 ๐‘›+ 2[ โˆ— = (๐‘‡ ๐‘ฟ โˆ’ ๐‘ฟโˆ—๐‘– ] . [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] ฮ”๐‘  3 ๐‘–+32 ๐‘–+1 1 1 ๐‘›+ ๐‘›+ 2 โˆ— โˆ’ 2๐‘‡ 1 [๐‘ฟ๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] . [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] + ๐‘‡ 2 1 [๐‘ฟ๐‘–+1 โˆ— โˆ’ ๐‘ฟโˆ—๐‘– ] . [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ]) ๐‘–+ ๐‘–โˆ’ 2 2 At mixed end, ๐‘– = 0, where vertical forced oscillation and horizontal unconstrained condition is combined. Along the horizontal direction, the LHS of Eq. (55) is discretized as, 21 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘‹ โˆ— ))1 โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘‹ โˆ— ))0 [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+1/2 ๐ท๐‘  ๐‘‹ โˆ— ))]1/2 = ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ 2( 2( ๐‘‡3 ๐ท๐‘ 0 ๐‘‹ โˆ— )3 โˆ’ 3๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 2 2 2 2 = ฮ”๐‘  2 ๐‘›+ 1 ๐‘›+ 1 (61) ๐‘‡3 2 [๐‘‹2โˆ— โˆ’ ๐‘‹1โˆ— ] โˆ’ 3๐‘‡1 2 [๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— ] 2 2 = ฮ”๐‘  3 The RHS of Eq. (61) is derived by applying the following steps, 1 1 ๐‘›+ ๐‘›+ 2( 2( ๐‘‡3 ๐ท๐‘ 0 ๐‘‹ โˆ— )3 โˆ’ ๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 1 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘‹ โˆ— )) = 2 2 1 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ 2[ โˆ— ๐‘‡3 ๐‘‹2 โˆ’ ๐‘‹1โˆ— ] โˆ’ ๐‘‡1 2[ โˆ— ๐‘‹1 โˆ’ ๐‘‹0โˆ— ] 2 2 = , ฮ”๐‘  2 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ 2( 2( 2( ๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 โˆ’ ๐‘‡0 ๐ท๐‘ 0 ๐‘‹ โˆ— )0 ๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘‹ โˆ— ))0 = 2 2 = 2 2 (62) ฮ”๐‘ /2 ฮ”๐‘ /2 1 ๐‘›+ 2๐‘‡1 2 [๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— ] 2 = , ฮ”๐‘  2 Along the vertical direction, the LHS of Eq. (55) is discretized as, (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ— ))1 โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ— ))0 [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+1/2 ๐ท๐‘  ๐‘Œ โˆ— ))]1/2 = ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )3 โˆ’ ๐‘‡1 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 2 2 2 2 = 2 +๐‘… ฮ”๐‘  ๐‘›+ 1 ๐‘›+ 1 (63) 2[ โˆ— ๐‘‡3 ๐‘Œ2 โˆ’ ๐‘Œ1โˆ— ] โˆ’ ๐‘‡1 2[ โˆ— ๐‘Œ1 โˆ’ ๐‘Œ0โˆ— ] 2 2 = +๐‘…, ฮ”๐‘  3 22 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ—))0 here ๐‘… = โˆ’ , and the RHS of Eq. (63) is derived by applying the following steps, ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )3 โˆ’ ๐‘‡1 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 1 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘Œ โˆ— )) = 2 2 1 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2 [๐‘Œ2โˆ— โˆ’ ๐‘Œ1โˆ— ] โˆ’ ๐‘‡1 2 [๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— ] 2 2 = , (64) ฮ”๐‘  2 ๐’ˆ(2) [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘Œ)]0 = [๐ท๐‘ ๐‘  (๐›พ๐ท๐‘ ๐‘  ๐‘Œ)]0 + (๐น (2) ) โˆ’ ๐›ฝ๐น๐‘Ÿ ( ) + ๐‘ฆฬˆ (๐‘ก) , 0 ๐‘” ๐œ• 2๐‘Œ where ๐‘ฆฬˆ (๐‘ก) = ( ๐œ•๐‘ก 2 ) . This ๐‘… term is free of T and hence, is shifted to the RHS of Eq. (55) and is 0 omitted in the LHS for the mixed end discretization. Then we get, 1 (๐ท๐‘ 0 ๐‘ฟโˆ— )1 . [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))]1 = 2 2 1 1 ๐‘›+ 2 ( 0 โˆ—) ( 0 โˆ—) = (๐‘‡ [ ๐ท๐‘  ๐‘‹ 1 ๐ท๐‘  ๐‘‹ 3 + (๐ท๐‘ 0 ๐‘Œ โˆ— )1 (๐ท๐‘ 0 ๐‘Œ โˆ— )3 ] ฮ”๐‘  2 32 2 2 2 2 1 ๐‘›+ 2 โˆ’ ๐‘‡1 [3(๐ท๐‘ 0 ๐‘‹ โˆ— )1 (๐ท๐‘ 0 ๐‘‹ โˆ— )1 + (๐ท๐‘ 0 ๐‘Œ โˆ— )1 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 ]) 2 2 2 2 2 (65) 1 1 ๐‘›+ = 3 (๐‘‡3 2 [(๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— )(๐‘‹2โˆ— โˆ’ ๐‘‹1โˆ— ) + (๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )(๐‘Œ2โˆ— โˆ’ ๐‘Œ1โˆ— )] ฮ”๐‘  2 1 ๐‘›+ 2 โˆ’ ๐‘‡1 [3(๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— )(๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— ) + (๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )(๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )]) 2 23 At free end, ๐‘– = ๐‘, (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘ โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘โˆ’1 [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+1/2 ๐ท๐‘  ๐‘ฟโˆ— ))]๐‘โˆ’1/2 = ฮ”๐‘  1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ โˆ’2๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 โˆ’ (๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 โˆ’ ๐‘‡ 32 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’3 ) ๐‘โˆ’ 2 ๐‘โˆ’ 2 ๐‘โˆ’ 2 2 2 2 = 2 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’3 โˆ’ 3๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 ๐‘โˆ’ 2 2 ๐‘โˆ’ 2 2 (66) = ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘โˆ’1 โˆ’ ๐‘ฟโˆ—๐‘โˆ’2 ] โˆ’ 3๐‘‡ 12 [๐‘ฟโˆ—๐‘ โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] ๐‘โˆ’ ๐‘โˆ’ 2 2 = ฮ”๐‘  3 The RHS of Eq. (66) is derived by applying the following steps, 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ 2( 2 2 ๐‘‡๐‘ ๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘ โˆ’ ๐‘‡ 0 โˆ— 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 ๐‘‡ 0 โˆ— 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 ๐‘โˆ’ 2 ๐‘โˆ’ 2 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘ = = โˆ’ ฮ”๐‘  ฮ”๐‘  2 2 1 ๐‘›+ 2๐‘‡ 12 [๐‘ฟโˆ—๐‘ โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] ๐‘โˆ’ 2 = โˆ’ , ฮ”๐‘  2 1 1 (67) ๐‘›+ ๐‘›+ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘‹ โˆ— )๐‘โˆ’1 โˆ’ ๐‘‡ 32 (๐ท๐‘ 0 ๐‘‹ โˆ— )๐‘โˆ’3 1 ๐‘โˆ’ 2 ๐‘โˆ’ 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘‹ โˆ— )) = 2 2 ๐‘โˆ’1 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡ 12 [๐‘‹๐‘โˆ— โˆ’ โˆ— ๐‘‹๐‘โˆ’1 ] โˆ’ ๐‘‡ 32 [๐‘‹๐‘โˆ’1 โˆ— โˆ— โˆ’ ๐‘‹๐‘โˆ’2 ] ๐‘โˆ’ ๐‘โˆ’ 2 2 = 2 . ฮ”๐‘  24 The LHS of Eq. (55) becomes, 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘โˆ’ 1 . [๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘โˆ’ 2 1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’3 โˆ’ 3๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 ๐‘โˆ’ 2 ๐‘โˆ’ 2 2 2 = (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 . 2 ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 . (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’3 โˆ’ 3๐‘‡ 0 โˆ— 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 . (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 (68) ๐‘โˆ’ 2 2 ๐‘โˆ’ 2 2 2 2 = ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘ โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] . [๐‘ฟโˆ—๐‘โˆ’1 โˆ’ ๐‘ฟโˆ—๐‘โˆ’2 ] โˆ’ 3๐‘‡ 2 1 [๐‘ฟ๐‘ โˆ— โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] . [๐‘ฟโˆ—๐‘ โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] ๐‘โˆ’ ๐‘โˆ’ 2 2 = ฮ”๐‘  3 The bending for is determined by (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–+1 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– + (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 (๐‘ญ๐‘ )๐‘– = โˆ’[๐ท๐‘ + ๐ท๐‘ โˆ’ (๐›พ๐ท๐‘ ๐‘  ๐‘ฟ)]๐‘– = โˆ’๐›พ , ฮ”๐‘  2 for ๐‘– = 1, 2, . . . , ๐‘ โˆ’ 1 and ๐›พ is constant , (1) (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘‹)0 (๐ท๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)0 (๐น๐‘ ) = โˆ’๐›พ = , 0 ฮ”๐‘  ฮ”๐‘  2 (69) (2) (๐ท๐‘ ๐‘ ๐‘  ๐‘Œ)1 โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘Œ)0 (๐ท๐‘ ๐‘  ๐‘Œ)2 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘Œ)1 + (๐ท๐‘ ๐‘  ๐‘Œ)0 (๐น๐‘ ) = โˆ’๐›พ = , 0 ฮ”๐‘  ฮ”๐‘  2 (๐ท๐‘ ๐‘ ๐‘  ๐‘ฟ)๐‘ โˆ’ (๐ท๐‘ ๐‘ ๐‘  ๐‘ฟ)๐‘โˆ’1 (๐ท๐‘ ๐‘  ๐‘‹)๐‘ โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)๐‘โˆ’1 (๐‘ญ๐‘ )๐‘ = โˆ’๐›พ = ๐›พ , ฮ”๐‘  ฮ”๐‘  2 25 4.2 Discretization of governing equation and boundary conditions for the self-propelled fish- like swimming filament with variable mass ratio and bending rigidity The governing equations for the motion of the flexible filament can be written in dimensionless form as, ๐œ•2๐‘ฟ 1 ๐œ• ๐œ•๐‘ฟ ๐œ•2 ๐œ•2๐‘ฟ โ‡’ 2 = [ (๐‘‡(๐‘ ) ) โˆ’ 2 (๐›พ(๐‘ ) 2 ) โˆ’ ๐‘ญ ] ๐œ•๐‘ก ๐›ฝ (๐‘ ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  (70) ๐œ•๐‘ฟ ๐œ•๐‘ฟ . =1, ๐œ•๐‘  ๐œ•๐‘  where ๐›ฝ(๐‘ ) and ๐›พ(๐‘ ) are exponential decay functions of ๐‘ , expressed as following, ๐‘š ๐›ฝ (๐‘ ) = ๐‘Ž1 ๐‘’ โˆ’๐‘1 ๐‘  (71) ๐‘› ๐›พ (๐‘ ) = ๐‘Ž2 ๐‘’ โˆ’๐‘2 ๐‘  The tension force is determined for our validation case in the same manner, by using the constraint of inextensibility constraint. The Poisson equation for ๐‘‡ can be derived as following, ๐œ•๐‘ฟ ๐œ•๐‘ฟ Differentiating ( ๐œ•๐‘  โ‹… ๐œ•๐‘  ) = 1, twice with respect to ๐‘ก, we get, ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ( โ‹… )=0, ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ โ‡’ [ ( โ‹… )] = 0 ๐œ•๐‘ก ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ โ‡’2 [ โ‹… ( )] = 0 ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ โ‡’ 2[ โ‹… 2( )+ ( )โ‹… ( ) ] = 0 ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ โ‡’2 โ‹… 2( )+2 ( )โ‹… ( ) = 0 ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  26 ๐œ•๐‘ฟ ๐œ• ๐œ• 2 ๐‘ฟ ๐œ•2๐‘ฟ ๐œ• 2 ๐‘ฟ โ‡’2 โ‹… ( )+2 โ‹… =0 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก 2 ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 1 ๐œ• ๐œ•๐‘ฟ ๐œ•2 ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ โ‡’2 โ‹… ( ( (๐‘‡ ๐‘  ( ) ) โˆ’ 2 (๐›พ ๐‘  ( ) ) โˆ’ ๐‘ญ) ) + 2 โ‹… =0 ๐œ•๐‘  ๐œ•๐‘  ๐›ฝ (๐‘ ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  2 ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐œ•2๐‘ฟ ๐œ• 2 ๐‘ฟ โ‡’2 โ‹… (๐›ฝ(๐‘ )โˆ’1 ( (๐‘‡(๐‘ ) ))) + 2 โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ)) + 2 โ‹… =0 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•2 ๐œ•2๐‘ฟ ๐œ•2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ where ๐‘ญ๐‘ = โˆ’ ๐œ•๐‘  2 (๐›พ ๐œ•๐‘  2 ). Since, ๐œ•๐‘ก 2 ( ๐œ•๐‘  โ‹… ๐œ•๐‘  ) = 0, we can write, ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ 2 โ‹… (๐›ฝ(๐‘ )โˆ’1 ( (๐‘‡(๐‘ ) ))) + 2 โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ)) + 2 โ‹… ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ = 2( โ‹… ) ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐œ•2๐‘ฟ ๐œ•2 ๐‘ฟ โ‡’ โ‹… (๐›ฝ(๐‘ )โˆ’1 ( (๐‘‡(๐‘ ) ))) + โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ)) + โ‹… ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  1 ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ = ( โ‹… ) 2 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  Then the Poisson equation for ๐‘‡(๐‘ ) can be written as ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ โ‹… (๐›ฝ(๐‘ )โˆ’1 ( (๐‘‡(๐‘ ) ))) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  (72) 1 ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ•๐‘ฟ ๐œ• = ( โ‹… )โˆ’ โ‹… โˆ’ โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ)) 2 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  27 Now the LHS of Eq. (72) can be expanded as ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ•2 ๐œ•๐‘ฟ โ‹… (๐›ฝ(๐‘ )โˆ’1 ( (๐‘‡(๐‘ ) ))) = โ‹… (๐›ฝ(๐‘ )โˆ’1 ( 2 (๐‘‡(๐‘ ) ))) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  (73) ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ•๐›ฝ (๐‘ )โˆ’1 + โ‹… ( (๐‘‡(๐‘ ) )) ( ), ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  which can be written in the discretized form as, 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘–+ 1 . [๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘–+ 2 1 1 = (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )] 1 (74) 2 ๐‘–+ 2 where the predicted position ๐‘ฟโˆ— = 2๐‘ฟ๐‘› โˆ’ ๐‘ฟ๐‘›โˆ’1 was used. 1 The term [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 can be discretized as ๐‘–+ 2 1 [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 ๐‘–+ 2 1 1 (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘ฟโˆ— )) โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘ฟโˆ— )) = ๐›ฝ (๐‘ )โˆ’11 ๐‘–+1 ๐‘– ๐‘–+ 2 ฮ”๐‘  1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘–+3 โˆ’ 2๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 +๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–+ 2 ๐‘–โˆ’ 2 = ๐›ฝ (๐‘ )โˆ’11 2 2 2 2 ๐‘–+ 2 ฮ”๐‘  (75) 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] โˆ’ 2๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] + ๐‘‡ 12 [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ] ๐‘–+ ๐‘–+ ๐‘–โˆ’ = ๐›ฝ (๐‘ )โˆ’11 2 2 2 . ๐‘–+ 2 ฮ”๐‘  3 28 The RHS of Eq. (75) is derived by applying the following steps, 1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— ๐‘‡ 3 (๐ท๐‘  ๐‘ฟ )๐‘–+3 โˆ’๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 ๐‘–+ 2 ๐‘–+ 2 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘–+1 = ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] โˆ’ ๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] ๐‘–+ ๐‘–+ 2 2 = ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ 2 0 โˆ— 2 0 โˆ— ๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–+1 โˆ’๐‘‡ 1 (๐ท๐‘  ๐‘ฟ )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–โˆ’ 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘– = 2 2 (76) ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] โˆ’ ๐‘‡ 12 [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ] ๐‘–+ ๐‘–โˆ’ 2 2 = 2 ฮ”๐‘  1 Then the term [(๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )] 1 can be discretized as ๐‘–+ 2 1 [(๐ท๐‘ + (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )] 1 ๐‘–+ 2 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ (๐‘‡ 32 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+3 โˆ’ 2๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 + ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–โˆ’1 ) (๐›ฝ(๐‘ )โˆ’1 ๐‘–+1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 ๐‘– ) ๐‘–+ 2 ๐‘–+ 2 ๐‘–โˆ’ 2 2 2 2 = 2ฮ”๐‘  2 1 1 1 (77) ๐‘›+ ๐‘›+ ๐‘›+ (๐‘‡ 32 [๐‘ฟโˆ—๐‘–+2 โˆ’ ๐‘ฟโˆ—๐‘–+1 ] โˆ’ 2๐‘‡ 12 [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] + ๐‘‡ 12 [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ]) (๐›ฝ(๐‘ )โˆ’1 ๐‘–+1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 ๐‘– ) ๐‘–+ ๐‘–+ ๐‘–โˆ’ 2 2 2 = 2ฮ”๐‘  3 So the LHS of Eq. (74) can be expressed as 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘–+ 1 . [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘–+ 2 1 1 = (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )] 1 2 ๐‘–+ 2 29 1 1 ๐‘›+ 2 ( 0 โˆ—) 0 โˆ—) ๐›ฝ(๐‘ )โˆ’1๐‘–+1 โˆ’1 ๐›ฝ (๐‘ )โˆ’1๐‘– = (๐‘‡ ( ๐ท ๐‘  ๐‘ฟ 1 . ( ๐ท ๐‘  ๐‘ฟ 3 ) ( + ๐›ฝ ( ๐‘  ) 1 โˆ’ ) ฮ”๐‘  2 ๐‘–+32 ๐‘–+ 2 ๐‘–+ 2 2 ๐‘–+ 2 2 1 ๐‘›+ โˆ’1 โˆ’ ๐‘‡ 12 ((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 ) (๐›ฝ (๐‘ )โˆ’1 ๐‘–+1 + 2๐›ฝ(๐‘ ) 1 โˆ’ ๐›ฝ (๐‘ )๐‘– ) โˆ’1 ๐‘–+ ๐‘–+ 2 2 2 2 (78) 1 ๐‘›+ ๐›ฝ (๐‘ )โˆ’1๐‘–+1 ๐›ฝ (๐‘ )โˆ’1 ๐‘– + ๐‘‡ 1 2 ((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–+1 . (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘–โˆ’1 ) ( + ๐›ฝ (๐‘ )โˆ’11 โˆ’ )) ๐‘–โˆ’ 2 2 2 2 ๐‘–+ 2 2 1 1 ๐‘›+ 2 ([ โˆ— โˆ—] [ โˆ— โˆ— ]) ๐›ฝ (๐‘ )โˆ’1 ๐‘–+1 โˆ’1 ๐›ฝ (๐‘ )โˆ’1 ๐‘– = (๐‘‡ ๐‘ฟ ๐‘–+1 โˆ’ ๐‘ฟ ๐‘– . ๐‘ฟ ๐‘–+2 โˆ’ ๐‘ฟ ๐‘–+1 ( + ๐›ฝ ( ๐‘  ) 1 โˆ’ ) ฮ”๐‘  3 ๐‘–+32 2 ๐‘–+ 2 2 1 ๐‘›+ โˆ’1 2 โˆ— โˆ’ ๐‘‡ 1 ([๐‘ฟ๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] . [๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ]) (๐›ฝ (๐‘ )โˆ’1 ๐‘–+1 + 2๐›ฝ(๐‘ ) 1 โˆ’ ๐›ฝ(๐‘ )๐‘– ) โˆ’1 ๐‘–+ ๐‘–+ 2 2 1 ๐‘›+ ๐›ฝ (๐‘ )โˆ’1 ๐‘–+1 โˆ’1 ๐›ฝ (๐‘ )โˆ’1 ๐‘– + ๐‘‡ 12 ([๐‘ฟโˆ—๐‘–+1 โˆ’ ๐‘ฟโˆ—๐‘– ] . [๐‘ฟโˆ—๐‘– โˆ’ ๐‘ฟโˆ—๐‘–โˆ’1 ]) ( ( +๐›ฝ ๐‘  1โˆ’ ) )) ๐‘–โˆ’ 2 2 ๐‘–+ 2 2 At mixed end, ๐‘– = 0, discretization of the LHS of Eq. (74) can be written as, 1 (๐ท๐‘ 0 ๐‘ฟโˆ— )1 . [๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))]1 2 2 1 1 = (๐ท๐‘ 0 ๐‘ฟโˆ— )1 . [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )]1 (79) 2 2 The RHS terms of (79) can be discretized horizontally as 1 1 ๐‘›+ ๐‘›+ 2( 2( ๐‘‡3 ๐ท๐‘ 0 ๐‘‹ โˆ— )3 โˆ’ 3๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 1 2 2 [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘‹ โˆ— ))]1 = ๐›ฝ (๐‘ )โˆ’1 1 2 2 2 ฮ”๐‘  2 2 1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2[ โˆ— ๐‘‹2 โˆ’ ๐‘‹1โˆ— ] โˆ’ 3๐‘‡1 2[ โˆ— ๐‘‹1 โˆ’ ๐‘‹0โˆ— ] (80) = ๐›ฝ (๐‘ )โˆ’1 1 2 2 2 ฮ”๐‘  3 30 1 [(๐ท๐‘ + (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘‹ โˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )]1 2 1 1 ๐‘›+ ๐‘›+ 2 (๐‘‡1 2 (๐ท๐‘ 0 ๐‘‹ โˆ— )1 โˆ’ ๐‘‡0 2 (๐ท๐‘ 0 ๐‘‹ โˆ— )0 ) (๐›ฝ(๐‘ )1โˆ’1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 0 ) 2 2 = ฮ”๐‘  2 1 ๐‘›+ 2( 2 (๐‘‡1 ๐ท๐‘ 0 ๐‘‹ โˆ— )1 ) (๐›ฝ(๐‘ )1โˆ’1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 0 ) 2 2 (81) = ฮ”๐‘  2 1 ๐‘›+ 2 (๐‘‡1 2 [๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— ]) (๐›ฝ(๐‘ )1โˆ’1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 0 ) 2 = ฮ”๐‘  3 By adding Eqs. (80) and (81), the following form can be written, 1 1 [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘‹ โˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘‹ โˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )]1 2 1 1 ๐‘›+ = 2 (๐‘‡3 2 ((๐ท๐‘ 0 ๐‘‹ โˆ— )3 ) (๐›ฝ (๐‘ )โˆ’1 1 ) ฮ”๐‘  2 2 2 1 ๐‘›+ โˆ’1 (82) โˆ’ ๐‘‡1 2 ((๐ท๐‘ 0 ๐‘‹ โˆ— )1 ) (2๐›ฝ(๐‘ )โˆ’1 โˆ’1 0 + 3๐›ฝ(๐‘ ) 1 โˆ’ 2๐›ฝ (๐‘ )1 )) 2 ๐‘–+ 2 2 1 1 ๐‘›+ = (๐‘‡ 2[ โˆ— ๐‘‹๐‘–+2 โˆ’ ๐‘‹๐‘–+1 โˆ— ] (๐›ฝ(๐‘ )โˆ’1 1 ) ฮ”๐‘  3 32 2 1 ๐‘›+ โˆ’1 2[ โˆ— โˆ’ ๐‘‡1 ๐‘‹๐‘–+1 โˆ’ ๐‘‹๐‘–โˆ— ] (2๐›ฝ(๐‘ )โˆ’1 โˆ’1 0 + 3๐›ฝ(๐‘ ) 1 โˆ’ 2๐›ฝ(๐‘ )1 )) . ๐‘–+ 2 2 31 1 Now the term [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘Œ โˆ— ))]1 can be discretized vertically as 2 1 [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘Œ โˆ— ))]1 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ— ))1 โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ— ))0 = ๐›ฝ (๐‘ )โˆ’1 1 2 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ 2( 2( ๐‘‡3 ๐ท๐‘ 0 ๐‘Œ โˆ— )3 โˆ’ ๐‘‡1 ๐ท๐‘ 0 ๐‘Œ โˆ— )1 2 2 = ๐›ฝ (๐‘ )โˆ’1 1 2 2 + ๐‘…๐‘  (83) 2 ฮ”๐‘  2 1 1 ๐‘›+ ๐‘›+ 2[ โˆ— ๐‘‡3 ๐‘Œ2 โˆ’ ๐‘Œ1โˆ— ] โˆ’ ๐‘‡1 2[ โˆ— ๐‘Œ1 โˆ’ ๐‘Œ0โˆ— ] = ๐›ฝ (๐‘ )โˆ’1 1 2 2 + ๐‘…๐‘  , 2 ฮ”๐‘  3 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘Œ โˆ— ))0 here ๐‘…๐‘  = โˆ’๐›ฝ (๐‘ )โˆ’1 1 , and the RHS of Eq. (83) was obtained by applying 2 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )3 โˆ’ ๐‘‡1 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 1 2 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘Œ โˆ— )) = 2 2 1 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡3 2 [๐‘Œ2โˆ— โˆ’ ๐‘Œ1โˆ— ] โˆ’ ๐‘‡1 2 [๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— ] 2 2 = , (84) ฮ”๐‘  2 ๐›ฝ (๐‘ )โˆ’1 โˆ’1 1 [๐ท๐‘  (๐‘‡๐ท๐‘  ๐‘Œ )]0 = ๐›ฝ (๐‘  )1 {[๐ท๐‘ ๐‘  (๐›พ๐ท๐‘ ๐‘  ๐‘Œ )]0 + (๐น (2) )0 + ๐‘ฆฬˆ (๐‘ก)} , 2 2 ๐œ• 2๐‘Œ where ๐‘ฆฬˆ (๐‘ก) = ( ๐œ•๐‘ก 2 ) . This ๐‘…๐‘  term is free of T and hence, is shifted to the RHS and is omitted in 0 1 the LHS for the mixed end discretization. And the term [(๐ท๐‘ + (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘Œ โˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )]1 can be 2 discretized vertically as 32 1 [(๐ท๐‘ + (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘Œ โˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )]1 2 1 1 ๐‘›+ ๐‘›+ 2 (๐‘‡1 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 โˆ’ ๐‘‡0 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )0 ) (๐›ฝ(๐‘ )1โˆ’1 โˆ’ ๐›ฝ(๐‘ )โˆ’1 0 ) 2 2 = ฮ”๐‘  2 1 1 (85) ๐‘›+ ๐‘›+ 2 (๐‘‡1 2 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 ) (๐›ฝ (๐‘ )1โˆ’1 โˆ’ ๐›ฝ(๐‘ )โˆ’10 ) 2 (๐‘‡1 2 [๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— ]) (๐›ฝ(๐‘ )1โˆ’1 โˆ’ ๐›ฝ (๐‘ )โˆ’1 0 ) 2 2 2 = = ฮ”๐‘  2 ฮ”๐‘  3 By adding Eqs. (83) and (85), the following form can be written, 1 1 [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘Œ โˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘Œ โˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )]1 2 1 1 ๐‘›+ = (๐‘‡ ( ๐ท๐‘  ๐‘Œ 3 ) (๐›ฝ (๐‘ )โˆ’1 2 ( 0 โˆ—) 1 ) ฮ”๐‘  2 32 2 2 1 ๐‘›+ โˆ’1 2 โˆ’ ๐‘‡1 ((๐ท๐‘ 0 ๐‘Œ โˆ— )1 ) (2๐›ฝ(๐‘ )โˆ’10 + ๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ(๐‘ )1 )) โˆ’1 2 ๐‘–+ 2 2 1 1 ๐‘›+ = 2 3 2 [๐‘‹๐‘–+2 (๐‘‡ โˆ— โˆ’ ๐‘‹๐‘–+1โˆ— ] (๐›ฝ (๐‘ )โˆ’11 ) ฮ”๐‘  2 2 (86) 1 ๐‘›+ โˆ’1 โˆ’ ๐‘‡1 2 [๐‘‹๐‘–+1โˆ— โˆ’ ๐‘‹๐‘–โˆ— ] (2๐›ฝ(๐‘ )โˆ’10 + ๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ (๐‘ )1 )) โˆ’1 ๐‘–+ 2 2 Then the LHS of Eq. (74) can be expressed as 1 (๐ท๐‘ 0 ๐‘ฟโˆ— )1 . [๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))]1 2 2 1 1 = (๐ท๐‘ 0 ๐‘ฟโˆ— )1 . [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )]1 2 2 33 1 1 ๐‘›+ = (๐‘‡ ( ๐ท๐‘  ๐‘‹ 1 ๐ท๐‘  ๐‘‹ 3 + (๐ท๐‘ 0 ๐‘Œ โˆ— )1 (๐ท๐‘ 0 ๐‘Œ โˆ— )3 ) (๐›ฝ (๐‘ )โˆ’1 2 ( 0 โˆ—) ( 0 โˆ—) 1 ) ฮ”๐‘  2 32 2 2 2 2 2 1 ๐‘›+ โˆ’1 โˆ’ ๐‘‡1 2 (((๐ท๐‘ 0 ๐‘‹ โˆ— )1 (๐ท๐‘ 0 ๐‘‹ โˆ— )1 ) (2๐›ฝ (๐‘ )โˆ’1 โˆ’1 0 + 3๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ(๐‘ )1 ) 2 2 ๐‘–+ 2 2 (87) โˆ’1 + ((๐ท๐‘ 0 ๐‘Œ โˆ— )1 (๐ท๐‘ 0 ๐‘Œ โˆ— )1 ) (2๐›ฝ(๐‘ )โˆ’10 + ๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ(๐‘ )1 ))) โˆ’1 2 2 ๐‘–+ 2 1 1 ๐‘›+ = (๐‘‡ [ ๐‘‹1 โˆ’ ๐‘‹0โˆ— )(๐‘‹2โˆ— โˆ’ ๐‘‹1โˆ— ) + (๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )(๐‘Œ2โˆ— โˆ’ ๐‘Œ1โˆ— )] (๐›ฝ(๐‘ )โˆ’1 2 ( โˆ— 1 ) ฮ”๐‘  2 32 2 1 ๐‘›+ โˆ’1 2 โˆ’ ๐‘‡1 ([(๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— )(๐‘‹1โˆ— โˆ’ ๐‘‹0โˆ— )] (2๐›ฝ (๐‘ )โˆ’1 โˆ’1 0 + 3๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ(๐‘ )1 ) ๐‘–+ 2 2 โˆ’1 + [(๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )(๐‘Œ1โˆ— โˆ’ ๐‘Œ0โˆ— )] (2๐›ฝ (๐‘ )โˆ’1 โˆ’1 0 + ๐›ฝ (๐‘ ) 1 โˆ’ 2๐›ฝ (๐‘ )1 ))) ๐‘–+ 2 At free end, ๐‘– = ๐‘, discretization of the LHS of Eq. (74) can be written as, 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘โˆ’ 1 . [๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 2 ๐‘โˆ’ 2 1 1 = (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 . [๐›ฝ(๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )] 1 (88) 2 ๐‘โˆ’ 2 34 The RHS terms of Eq. (88) can be discretized as 1 [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— ))] 1 ๐‘โˆ’ 2 (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘ โˆ’ (๐ท๐‘ 0 (๐‘‡ ๐‘›+1/2 ๐ท๐‘ 0 ๐‘ฟโˆ— ))๐‘โˆ’1 = ๐›ฝ(๐‘ )โˆ’1 1 ๐‘โˆ’ 2 ฮ”๐‘  1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ โˆ’2๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 โˆ’ (๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 โˆ’ ๐‘‡ 32 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’3 ) ๐‘โˆ’ 2 ๐‘โˆ’ 2 ๐‘โˆ’ 2 = ๐›ฝ(๐‘ )โˆ’1 1 2 2 2 2 ๐‘โˆ’ 2 ฮ”๐‘  1 1 (89) ๐‘›+ ๐‘›+ ๐‘‡ 32 (๐ท๐‘ 0 ๐‘ฟโˆ— )3 โˆ’ 3๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )1 ๐‘โˆ’ 2 ๐‘โˆ’ 2 = ๐›ฝ(๐‘ )โˆ’1 1 2 2 2 ๐‘โˆ’ 2 ฮ”๐‘  1 1 ๐‘›+ ๐‘›+ ๐‘‡ 32 [๐‘ฟโˆ—๐‘โˆ’1 โˆ’ ๐‘ฟโˆ—๐‘โˆ’2 ] โˆ’ 3๐‘‡ 12 [๐‘ฟโˆ—๐‘ โˆ’ ๐‘ฟโˆ—๐‘โˆ’1 ] ๐‘โˆ’ ๐‘โˆ’ = ๐›ฝ(๐‘ )โˆ’1 1 2 2 ๐‘โˆ’ 2 ฮ”๐‘  3 1 [(๐ท๐‘ + (๐‘‡ ๐‘›+2 ๐ท๐‘ 0 ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ (๐‘ )โˆ’1 )] 1 ๐‘โˆ’ 2 1 1 ๐‘›+ ๐‘›+ 2 (๐‘‡๐‘ 2 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘ โˆ’ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 ) (๐›ฝ(๐‘ )โˆ’1 ๐‘ โˆ’ ๐›ฝ(๐‘ )โˆ’1 ๐‘โˆ’1 ) ๐‘โˆ’ 2 2 = ฮ”๐‘  2 1 ๐‘›+ 2 โˆ’1 2 (๐‘‡ 0 โˆ— 1 (๐ท๐‘  ๐‘ฟ )๐‘โˆ’1 ) (๐›ฝ(๐‘ )๐‘ โˆ’ ๐›ฝ(๐‘ )โˆ’1๐‘โˆ’1 ) ๐‘โˆ’ 2 2 (90) = โˆ’ ฮ”๐‘  2 1 ๐‘›+ 2 (๐‘‡ 12 [๐‘‹๐‘โˆ— โˆ— โˆ’ ๐‘‹๐‘โˆ’1 ]) (๐›ฝ(๐‘ )โˆ’1 โˆ’1 ๐‘ โˆ’ ๐›ฝ(๐‘ )๐‘โˆ’1 ) ๐‘โˆ’ 2 =โˆ’ ฮ”๐‘  3 35 Then the LHS of Eq. (74) can be written as 1 1 (๐ท๐‘ 0 ๐‘ฟโˆ— ) ๐‘โˆ’ 1 . [๐›ฝ โˆ’1 ๐ท๐‘ 0 (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) + (๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟโˆ— )) ๐ท๐‘ 0 (๐›ฝ โˆ’1 )] 1 2 ๐‘โˆ’ 2 1 1 ๐‘›+ 2 โˆ’1 = (๐ท๐‘ 0 ๐‘ฟโˆ— ) 1 . 2 (๐‘‡ 0 โˆ— 3 ((๐ท๐‘  ๐‘ฟ )๐‘โˆ’3 ) (๐›ฝ๐‘โˆ’1 ) ๐‘โˆ’ 2 ฮ”๐‘  ๐‘โˆ’ 2 2 2 1 ๐‘›+ โˆ’ ๐‘‡ 12 ((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 ) (2๐›ฝ๐‘โˆ’1 + 3๐›ฝ โˆ’1 1 โˆ’ 2๐›ฝ๐‘โˆ’1 โˆ’1 )) ๐‘โˆ’ 2 ๐‘โˆ’ 2 2 1 1 ๐‘›+ = 2 (๐‘‡ 32 (((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 ) . ((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’3 )) (๐›ฝ โˆ’1 1 ) ฮ”๐‘  ๐‘โˆ’ 2 2 2 ๐‘โˆ’ 2 1 ๐‘›+ (91) โˆ’ ๐‘‡ 12 (((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 ) . ((๐ท๐‘ 0 ๐‘ฟโˆ— )๐‘โˆ’1 )) (2๐›ฝ๐‘โˆ’1 + 3๐›ฝ โˆ’1 1 โˆ’ 2๐›ฝ๐‘โˆ’1 โˆ’1 )) ๐‘โˆ’ 2 2 ๐‘โˆ’ 2 2 1 1 ๐‘›+ = 3 (๐‘‡ 32 [(๐‘ฟ1โˆ— โˆ’ ๐‘ฟโˆ—0 ) . (๐‘ฟโˆ—2 โˆ’ ๐‘ฟ1โˆ— )] (๐›ฝ โˆ’1 1 ) ฮ”๐‘  ๐‘โˆ’ 2 ๐‘โˆ’ 2 1 ๐‘›+ โˆ’ ๐‘‡ 12 [(๐‘ฟ1โˆ— โˆ’ ๐‘ฟโˆ—0 ) . (๐‘ฟ1โˆ— โˆ’ ๐‘ฟโˆ—0 )] (2๐›ฝ๐‘โˆ’1 + 3๐›ฝ โˆ’1 1 โˆ’ 2๐›ฝ๐‘โˆ’1 โˆ’1 )) . ๐‘โˆ’ ๐‘โˆ’ 2 2 Now the 3rd term from RHS of Eq. (72) is ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐›ฝ (๐‘ )โˆ’1 (92) โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ)) = โ‹… (๐›ฝ(๐‘ )โˆ’1 (๐‘ญ๐‘ โˆ’ ๐‘ญ) + (๐‘ญ๐‘ โˆ’ ๐‘ญ) ), ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘† which can be discretized as (๐ท๐‘ 0 ๐‘ฟโˆ— ) 1 . [๐›ฝ (๐‘ )โˆ’1 ๐ท๐‘ 0 (๐‘ญ๐‘ โˆ’ ๐‘ญ) + (๐‘ญ๐‘ โˆ’ ๐‘ญ)๐ท๐‘ 0 (๐›ฝ(๐‘ )โˆ’1 )]๐‘–+1 , (93) ๐‘–+ 2 2 36 The bending force ๐‘ญ๐‘ , where ๐›พ(๐‘ ) is an exponential decay function of ๐‘ , can be derived as ๐œ•2 ๐œ•2๐‘ฟ ๐œ• ๐œ• ๐œ•2๐‘ฟ (๐‘ญ๐‘ )๐‘– = โˆ’ 2 (๐›พ(๐‘ ) 2 ) = โˆ’ ( (๐›พ(๐‘ ) 2 )) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ• ๐œ• 3 ๐‘ฟ ๐œ•๐›พ (๐‘ ) ๐œ• 2 ๐‘ฟ =โˆ’ (๐›พ(๐‘ ) 3 + ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  2 ๐œ• ๐œ• 4 ๐‘ฟ ๐œ•๐›พ (๐‘ ) ๐œ• 3 ๐‘ฟ ๐œ•๐›พ (๐‘ ) ๐œ• 3 ๐‘ฟ ๐œ• 2 ๐›พ (๐‘ ) ๐œ• 2 ๐‘ฟ (94) =โˆ’ (๐›พ(๐‘ ) 4 + + + ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  3 ๐œ•๐‘  ๐œ•๐‘  3 ๐œ•๐‘  2 ๐œ•๐‘  2 ๐œ• ๐œ•4๐‘ฟ ๐œ•๐›พ (๐‘ ) ๐œ• 3 ๐‘ฟ ๐œ• 2 ๐›พ (๐‘ ) ๐œ• 2 ๐‘ฟ =โˆ’ (๐›พ(๐‘ ) 4 + 2 + ), ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  3 ๐œ•๐‘  2 ๐œ•๐‘  2 which can be discretized as, (๐‘ญ๐‘ )๐‘– = โˆ’[๐ท๐‘ + ๐ท๐‘ โˆ’ (๐›พ๐ท๐‘ ๐‘  ๐‘ฟ)]๐‘– (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–+1 โˆ’ 2(๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– + (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 = โˆ’๐›พ๐‘– ฮ”๐‘  2 (๐›พ๐‘–+1 โˆ’ ๐›พ๐‘–โˆ’1 ) ((๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– โˆ’ (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘–โˆ’1 ) (95) 2 2 โˆ’ 2[ ] ฮ”๐‘  2 (๐›พ๐‘–+1 โˆ’ 2๐›พ๐‘– + ๐›พ๐‘–โˆ’1 ) (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘– โˆ’ ฮ”๐‘  2 At mixed end, ๐‘– = 0 (1) (๐ท๐‘ ๐‘  ๐‘‹)1 โˆ’ (๐ท๐‘ ๐‘  ๐‘‹)0 (โˆ’๐›พ3 + 4๐›พ2 โˆ’ 5๐›พ1 + 2๐›พ0 ) (๐ท๐‘ ๐‘  ๐‘‹)0 (๐น๐‘ ) = โˆ’๐›พ0 โˆ’ 0 ฮ”๐‘  2 ฮ”๐‘  2 (2) (๐ท๐‘ ๐‘  ๐‘Œ)2 โˆ’2(๐ท๐‘ ๐‘  ๐‘Œ)1 + (๐ท๐‘ ๐‘  ๐‘Œ)0 (๐›พ1 โˆ’ ๐›พ0 ) ((๐ท๐‘ ๐‘  ๐‘Œ)1 โˆ’ (๐ท๐‘ ๐‘  ๐‘Œ)0 ) (๐น๐‘ ) = โˆ’๐›พ0 โˆ’ 2 [ ] (96) 0 ฮ”๐‘  2 ฮ”๐‘  2 (โˆ’๐›พ3 + 4๐›พ2 โˆ’ 5๐›พ1 + 2๐›พ0 ) (๐ท๐‘ ๐‘  ๐‘Œ)0 โˆ’ ฮ”๐‘  2 37 At free end, ๐‘– = ๐‘ (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘โˆ’1 โˆ’ (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘ (2๐›พ๐‘ โˆ’ 5๐›พ๐‘โˆ’1 + 4๐›พ๐‘โˆ’2 โˆ’ ๐›พ๐‘โˆ’3 ) (๐ท๐‘ ๐‘  ๐‘ฟ)๐‘ (97) (๐‘ญ๐‘ )๐‘ = โˆ’๐›พ๐‘ โˆ’ ฮ”๐‘  2 ฮ”๐‘  2 Now, from Eq. (70) the following discretized form can be obtained, ๐‘ฟ๐‘›+1 ๐‘– โˆ’ 2๐‘ฟ๐‘›๐‘– + ๐‘ฟ๐‘›โˆ’1 ๐‘– ฮ”๐‘ก 2 1 = ([๐ท๐‘  (๐‘‡ ๐‘›+2 ๐ท๐‘  ๐‘ฟ๐‘›+1 )] + (๐‘ญโˆ—๐‘ )๐‘– โˆ’ (๐‘ญ๐‘› )๐‘– ) ๐›ฝ๐‘–โˆ’1 ๐‘– 1 1 ๐‘›+ ๐‘›+ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟ๐‘›+1 )๐‘–+1 โˆ’ ๐‘‡ 12 (๐ท๐‘ 0 ๐‘ฟ๐‘›+1 )๐‘–โˆ’1 ๐‘–+ 2 ๐‘–โˆ’ 2 2 2 = + (๐‘ญโˆ—๐‘ )๐‘– โˆ’ (๐‘ญ๐‘› )๐‘– ๐›ฝ๐‘–โˆ’1 ฮ”๐‘  ( ) 1 1 ๐‘›+ ๐‘›+ ๐‘‡ 12 (๐‘ฟ๐‘›+1๐‘–+1 โˆ’ ๐‘ฟ๐‘– ) โˆ’ ๐‘‡ 12 (๐‘ฟ๐‘›+1 ๐‘›+1 ๐‘– โˆ’ ๐‘ฟ๐‘›+1๐‘–โˆ’1 ) (98) ๐‘–+ ๐‘–โˆ’ 2 2 = + (๐‘ญโˆ—๐‘ )๐‘– โˆ’ (๐‘ญ๐‘› )๐‘– ๐›ฝ๐‘–โˆ’1 , ฮ”๐‘  2 ( ) for ๐‘– = 0, 1, 2, . . . , ๐‘ , where the bending force, ๐‘ญโˆ—๐‘ is calculated explicitly as ๐‘ญโˆ—๐‘ = ๐‘ญ๐‘ (2๐‘ฟ๐‘› โˆ’ ๐‘ฟ๐‘›โˆ’1 ) . Then finally we can obtain ๐‘ฟ๐‘›+1 from the following equation, 1 1 1 1 ๐‘›+ ๐‘›+ ๐‘›+ ๐‘›+ ๐‘‡ 12 ๐‘‡ 12 ๐‘‡ 12 ๐‘‡ 12 ๐‘–+ 2 ๐›ฝ๐‘– ๐‘–+ 2 ๐‘–โˆ’ 2 ๐‘–โˆ’ 2 โˆ’ ๐‘ฟ๐‘›+1 ๐‘–+1 + ๐‘ฟ๐‘›+1๐‘– + + โˆ’ ๐‘ฟ๐‘›+1 ๐‘–โˆ’1 ฮ”๐‘  2 ฮ”๐‘ก 2 ฮ”๐‘  2 ฮ”๐‘  2 ฮ”๐‘  2 ( ) ( ) ( ) (99) ๐›ฝ๐‘– = (2๐‘ฟ๐‘›๐‘– โˆ’ ๐‘ฟ๐‘›โˆ’1 ๐‘– ) + (๐‘ญโˆ—๐‘ )๐‘– โˆ’ (๐‘ญ๐‘› )๐‘– ฮ”๐‘ก 2 38 5 Results and Discussion 5.1 A pitching flexible filament with fixed leading edge The flexible deformation of a thin pitching filament was simulated in this case. A pitching motion, prescribed by ๐œƒ (๐‘ก) = ๐œƒ0 ๐‘ ๐‘–๐‘›(2๐œ‹๐‘“๐‘ก), was applied on the leading edge, which was at the origin, of the flexible filament (shown in figure 5-1). Here ๐œƒ0 is the pitching amplitude and ๐‘“ is the pitching frequency. The computational domain used in the simulation was [-2L,14L] ร— [-3L,3L], uniform Cartesian mesh grid was used, and the grid spacing was L/50. The time step was 0.001๐ฟ/๐‘ˆ in the IB-LB solver and the results were compared to the DSD/SST solver by Xu et al. [74], where the time step was 1/(500๐‘“). Figure 5-1: Sketch of a pitching flexible filament in a uniform flow The bending rigidities used in the simulation are ๐พ๐ต = 0.125 and ๐พ๐ต = 0.125. The stretching ๐‘“๐ฟ coefficient, ๐พ๐‘  = 500, Reynolds number, ๐‘…๐‘’ = 100, ๐‘“ โˆ— = = 0.6, and ๐œƒ0 = 30ยฐ were taken as ๐‘ˆ control parameters for this case. This case was earlier used to validate the code used for simulation, for fixed leading edge-flexible plate using IB-LB solver [74]. With minor modifications in the code, the present study also validated this case with IB-LB solver, before moving towards free swimming case. 39 (a) (b) (c) (d) (e) (f) (g) (h) Figure 5-2: (a) Drag coefficients, (c) lift coefficients, (e) X-coordinates and (g) Y-coordinates of a pitching filament for ๐พ๐ต = 0.125, and (b) drag coefficients, (d) lift coefficients, (f) X-coordinates and (h) Y-coordinates of the pitching filament for ๐พ๐ต = 0.0625. 40 The results from DSD/SST method [74] are consistent with the IB-LBM. The simulations were continued until stability is obtained by the trailing node of the pitching filament. Results for ๐พ๐ต = 0.125 were more consistent with maximum error margin of 3%, while results for ๐พ๐ต = 0.0625 had maximum error margins of around 11% for flapping amplitude. For average drag, the maximum error margin was 23% for ๐พ๐ต = 0.125, and 28% for ๐พ๐ต = 0.0625. The deformation patterns of the pitching filament are shown in figure 5-3. (a) (b) Figure 5-3: Flapping pattern of a pitching filament in a uniform flow for (a) ๐พ๐ต = 0.125 and (b) ๐พ๐ต = 0.0625. The patterns are plotted for 66 instants at rime interval of 0.125. 41 (a) (b) (c) (d) (e) (f) (g) (h) Figure 5-4: Vorticity fields at four different instants during a flapping cycle of a 2D pitching filament for ๐พ๐ต = 0.125 (left column) and ๐พ๐ต = 0.125 (right column). The corresponding color bars are shown at the bottom of each column. 42 5.2 Single swimming filament with heaving motion A single filament was actuated by a sinusoidal heaving motion, without any pitching, at the leading edge (schematic is shown in figure 5-5). The movement along horizontal direction was unconstrained, while the vertical displacement was imposed at leading edge in the following forms respectively: ๐‘ฆ(๐‘ก) = ๐‘ฆ0 + ๐ด๐‘๐‘œ๐‘ (2๐œ‹๐‘“๐‘ก + ๐œ™), where ๐‘ฆ0 refers to the equilibrium lateral position of heaving motion, ๐ด is the heaving amplitude, ๐‘“ is the actuation frequency and ๐œ™ is the phase angle. The governing equations for fluid flow can be written in the form of Eq. (6). In this particular case, ๐‘…๐‘’ is the flapping Reynolds number defined as ๐‘…๐‘’ = (๐‘ˆref ๐ฟ)/๐œ. Here, ๐‘ˆref is the reference velocity, L is the length of the filament and ๐œ is the kinematic viscosity of the fluid. The reference velocity is defined as ๐‘ˆref = 2๐œ‹๐œƒ0 ๐‘“๐ฟ, bending coefficient, ๐พ๐ต = 0.125 and stretching coefficient, ๐พ๐‘  = 500. The tension coefficient is calculated by a Poisson-type equation to apply inextensibility constraint. The simulation was modelled using IB-LB method. Figure 5-5: Sketch of a model for locomotion of a swimming filament with only heaving motion applied on the leading edge. Table 5-1: Control Parameters for only heaving motion applied on the leading edge Control parameters Values Flapping Reynolds number (๐‘…๐‘’ = ๐‘ˆref ๐ฟ/๐œ) 600 Dimensionless heaving amplitude (๐ด/๐ฟ) 0.1 43 The swimming pattern for heaving only case with the described control parameters is shown in figure 5-6. Figure 5-6: Self-propelled swimming filament with oscillatory heaving motion applied on the leading edge. The computational domain used in the simulation was [-8L,8L] ร— [-3L,3L], uniform Cartesian mesh grid was used, and the grid spacing was L/25. (a) (b) Figure 5-7: Drag Coefficient (๐ถD ) and lift coefficient (๐ถL ) for a swimming filament with heaving motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). 44 Figure 5-8: Y-coordinate of the trailing point of a swimming filament with heaving motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). 45 (a) (b) (c) (d) (e) (f) Figure 5-9: Vorticity at six different instants during a period of 2D swimming filament with only heaving motion applied on the leading edge. The corresponding color bars are shown at the bottom of each column. The pressure on the surface was calculated at equidistant nodes on both sides. While calculating pressure, the upper surface was considered to be in the outward direction, and the lower surface was considered to be in the outward direction. Discreet marker can be used on the surface nodes, according to required node numbers, to obtain the surface pressure. The surface pressure on the leading edge, mid-point and trailing edge for a swimming filament with heaving motion applied on the leading edge is shown in figure 5-10. 46 (a) (b) (c) Figure 5-10: Surface pressure on the swimming filament, with heaving motion applied on the head, at (a) leading edge, (b) trailing edge, and (c) middle point. The pressure on the upper and lower surface of the filament is sinusoidal and anti-symmetric. The surface pressure measured on nodes becomes smoother as they move towards the center of the filament length. 47 5.3 Single swimming filament with heaving and pitching motion A thin elastic filament was used to replicate undulatory motion of a fish. The swimmer was actuated by a sinusoidal heaving and pitching motion at the leading edge. The movement along horizontal direction was unconstrained, while the vertical displacement was imposed at leading edge in the following forms respectively: ๐‘ฆ(๐‘ก) = ๐‘ฆ0 + ๐ด๐‘๐‘œ๐‘ (2๐œ‹๐‘“๐‘ก) and ๐œƒ(๐‘ก) = ๐œƒ0 ๐‘๐‘œ๐‘ (2๐œ‹๐‘“๐‘ก + ๐œ™ โˆ’ ๐œ‹/2), where ๐‘ฆ0 refers to the equilibrium lateral position of heaving motion, ๐ด is the heaving amplitude and ๐œƒ0 is the pitching amplitude, ๐‘“ is the actuation frequency and ๐œ™ is the phase angle, ๐‘…๐‘’ is the flapping Reynolds number defined as ๐‘…๐‘’ = (๐‘ˆref ๐ฟ)/๐œ. Here, ๐‘ˆref is the reference velocity, L is the length of the filament and ๐œ is the kinematic viscosity of the fluid. The reference velocity is related to flapping frequency by ๐‘ˆref = 2๐œ‹๐œƒ0 ๐‘“๐ฟ, the bending coefficient, ๐พ๐ต = 0.125 and the stretching coefficient, ๐พ๐‘  = 500. The tension coefficient is calculated by a Poisson-type equation to apply inextensibility constraint. The simulation was modelled using IB-LB method. Table 5-2: Control Parameters for heaving and pitching motion applied on the leading edge Control parameters values Flapping Reynolds number (๐‘…๐‘’ = ๐‘ˆref ๐ฟ/๐œ) 600 Dimensionless heaving amplitude (๐ด/๐ฟ) 0.1 Pitching amplitude (๐œƒ0 ) 5ยฐ Phase angle (๐œ™) ๐œ‹/2 The swimming pattern for only heaving case with the listed parameters is shown in figure 5-11. Figure 5-11: Self-propelled swimming filament with heaving and pitching motion applied on the leading edge. 48 The computational domain used in the simulation was [-8L,8L] ร— [-3L,3L], uniform Cartesian mesh grid was used, and the grid spacing was L/25. (a) (b) Figure 5-12: Drag Coefficient (๐ถD ) and lift coefficient (๐ถL ) for a swimming filament with heaving and pitching motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). Figure 5-13: Y-coordinate of the trailing point of a swimming filament with heaving and pitching motion applied on the leading edge. The results were obtained at three different grid resolutions (L/100, L/50 and L/25). 49 (a) (b) (c) (d) (e) (f) Figure 5-14: Vorticity at six different instants during a period of 2D swimming filament with heaving and pitching motion applied on the leading edge. The corresponding color bars are shown at the bottom of each column. 50 The pressure on the surface was calculated at equidistant nodes on both sides. The surface pressure on the leading edge, middle-point and trailing edge for a swimming filament with heaving and pitching motion applied on the leading edge is shown in figure 5-15. (a) (b) (c) Figure 5-15: Surface pressure on the swimming filament, with heaving and pitching motion applied on the head, at (a) leading edge, (b) trailing edge, and (c) middle point. 51 The pressure on the upper and lower surface of the filament is observed as sinusoidal and anti- symmetric. The pressure data was generated at multiple nodes and a smooth fitted graph can be used for optimization. 52 5.4 Self-propelled fish-like swimmer with variable mass ratio and bending rigidity The numerical algorithm described in Chapter 4 was validated by previously established results from the work of Dai et al. [72] on self-propelled fish like swimmers. A thin elastic filament was used to replicate undulatory motion of a fish. The swimmer was actuated by a sinusoidal heaving and pitching motion at the leading edge. The movement along horizontal direction was unconstrained, while the vertical displacement was imposed at leading edge in the following forms respectively: ๐‘ฆ(๐‘ก) = ๐‘ฆ0 + ๐ด๐‘๐‘œ๐‘ (2๐œ‹๐‘“๐‘ก + ๐œ™) and ๐œƒ (๐‘ก) = ๐œƒ0 ๐‘๐‘œ๐‘ (2๐œ‹๐‘“๐‘ก + ๐œ™ โˆ’ ๐œ‹/2), where ๐‘ฆ0 refers to the equilibrium lateral position of heaving motion, ๐ด is the heaving amplitude and ๐œƒ0 is the pitching amplitude, ๐‘“ is the actuation frequency and ๐œ™ is the phase angle. The governing equations for fluid flow can be written in the form of Eq. (6). In this particular case, ๐‘…๐‘’ is the flapping Reynolds number defined as ๐‘…๐‘’ = (๐‘ˆref ๐ฟ)/๐œ. Here, ๐‘ˆref is the reference velocity, ๐ฟ is the length of the filament and ๐œ is the kinematic viscosity of the fluid. The reference velocity is defined as ๐‘ˆref = 2๐œ‹๐œƒ0 ๐‘“๐ฟ. The governing equation for flexible filament motion is expressed in Eq. (71), where, ๐‘ฟ denotes the position vector, ๐‘  denotes the Lagrangian coordinate, and ๐‘ญ refers to the Lagrangian forcing due to hydrodynamic interaction with fluid. The mass ratio ๐›ฝ and dimensionless bending rigidity ๐›พ are exponential decay function of ๐‘  (from leading edge to trailing edge), expressed in Eq. (72). The schematic of the computational model is shown in figure 5-16. Figure 5-16: Schematic for computational model with heaving and pitching imposed on the leading edge The swimming pattern on the filament model is obtained by solving the FSI problem numerically using LB-IBM. The control parameter values used in the simulation are listed in table 5-3 and the coefficients in the decay functions are listed on table 5-4. 53 Table 5-3: Control parameter values used in the simulation for self-propelled fish-like swimmer with variable mass ratio and bending rigidity Control parameters Values Flapping Reynolds number (๐‘…๐‘’ = ๐‘ˆref ๐ฟ/๐œ) 600 Dimensionless heaving amplitude (๐ด/๐ฟ) 0.005 Pitching amplitude (๐œƒ0 ) 5ยฐ Phase angle (๐œ™) ๐œ‹/2 Table 5-4: Coefficients in the exponential decay functions of ฮฒ(s) and ฮณ(s) ๐‘Ž1 0.1 ๐‘1 3.0 ๐‘š 2.0 ๐‘Ž2 1.0 ๐‘2 9.0 ๐‘› 2.0 The detailed derivative and discretization are stated in Appendix C. The comparison of swimming patterns between a red nose tetra fish Hemigrammus bleheri [73], fish-like self-propelled swimmer by Dai et al. [72] and self-propelled swimming filament using LB-IBM method is shown in figure 5-17. 54 (a) (b) (c) Figure 5-17: The comparison of swimming patterns between a red nose tetra fish Hemigrammus bleheri [73], fish- like self-propelled swimmer by Dai et al. [72] and self-propelled swimming filament using LB-IBM method The computational domain used in the simulation was [-30L,30L] ร— [-6L,6L], uniform Cartesian mesh grid was used, and the grid spacing was L/50. In figure 5-17, the superimposed instantaneous middle-line locomotion is plotted. Here, the spatial envelop (thick black lines) are fitted with analytical function of ๐ด(๐‘ฅ ) = ๐ด๐‘Ÿ ๐‘’๐‘ฅ๐‘[๐›ผ(๐‘ฅ โˆ’ 1)]. In (b), the kinematics are plotted for 55 16 instants at time interval ๐‘‡๐‘“ /16, where ๐‘‡๐‘“ was referred to flapping period. In (c), the kinematic are plotted at time interval t/๐‘‡๐‘Ÿ๐‘’๐‘“ = 0.125, where ๐‘‡๐‘Ÿ๐‘’๐‘“ was reference time. (a) (b) Figure 5-18: (a) The time histories of instantaneous swimming speed of the central node of the filament, and (b) drag and lift coefficients of the swimming filament for control parameter and coefficients listed in Table 5-3 and Table 5-4. Table 5-5: Comparison between the kinematics of a red nose tetra fish Hemigrammus bleheri [73], fish-like self- propelled swimmer by Dai et al. [72] and self-propelled swimming filament using LB-IBM method Kinematics Red nose tetra Dai et al. LB-IBM fish Tail Amplitude ๐ด 0.138 0.133 0.137 Strouhal number ๐‘†๐‘ก 0.503 0.325 0.500 Dimensionless cruising speed 0.473 0.740 0.475 ๐‘ˆ๐‘ 56 (a) (b) (c) (d) (e) (f) (g) (h) Figure 5-19: Vorticity fields (left column) and pressure fields (right column) at four different instants during a period of a self-propelled fish-like swimmer. The corresponding color bars are shown at the bottom of each column. 57 The pressure on the surface was calculated at equidistant nodes on both sides. The surface pressure on the leading edge, middle point and trailing edge is shown in figure 5-20. (a) (b) (c) Figure 5-20: Surface pressure on the filament at (a) leading edge, (b) trailing edge, and (c) middle point. The pressure on the upper and lower surface of the filament is anti-symmetric. The pressure on nodes around mid-region of the filament length is smoother than more fluctuating near the leading and trailing edges. 58 6 Conclusion This dissertation presents a numerical investigation of a freely swimming filament with actuating force applied on the leading edge. We focus on the detailed derivation and discretization of the governing equations and boundary conditions for flow and structural dynamics. A uniform Cartesian grid is used for fluid domain and a staggered grid is used in the Lagrangian coordinate system, where tension force is defined on the interfaces (half-grids) and other variables are defined on the nodes. The immersed boundary-lattice Boltzmann method (IB-LBM) is applied for solving the fluid-structural interaction problems. In the immersed boundary method, the structure is considered to be immersed into the fluid domain. A set of discreet markers were employed to express the immersed boundary geometry. The Lagrangian force is exerted on the neighboring fluid nodes as a body force, and thus the effects of the immersed boundary are accounted for. The fluid body force density is calculated from the Lagrangian force density and a smoothed approximation of Dirac delta function, where the Lagrangian force is determined from the structural velocity and interpolated velocity flow field. In the lattice Boltzmann method, the fluid is modelled to be consisting of fictive particles, which consecutive collision and propagation processes over a discreet lattice structure. The discreet lattice Boltzmann equation of a single time relaxation (BGK) is used to solve the Navier-Stokes equations. The density distribution function is calculated from non-dimensional relaxation time, discreet velocity components and equilibrium distribution function, and then this term is used to solve for fluid density, velocity, and pressure. Oscillatory heaving and pitching motions are imposed at the leading edge of a thin flexible filament with prescribed control parameters. We investigated the flow physics and hydrodynamic coupling of the structure. Surface pressure on the filament is determined from the flow and structural dynamics. The swimming gaits of the filament is analyzed and different output parameters including drag and lift coefficients, X and Y- coordinate of the trailing edge, and cruising speed of the filament are determined. The results obtained in this study shows good correspondence with the relevant observations. The results are generated for different grid spacings, and they show consistency. 59 The presented study and pressure surface calculation scheme may be used in future with multiple obstacles in the domain. One or multiple bluff body structures may be positioned in the fluid domain and the change in swimming filament dynamics may be investigated. The position of the obstacle and the swimming angle of the filament, with respect to the bluff body structure, can be varied, and corresponding pressure data sets may be generated for control optimization. The present study may be remodeled with different control parameters. The Reynolds number may be varied, and the swimming patterns and flow physics may be investigated further. The results may be further analyzed with experimental results in future for different geometry considerations. The numerical simulation used in this study may be extended to 3-D geometry and D3Q19 model may be used in LBM instead of the present D2Q9 model. Instead of a thin filament, multiple side-by-side filaments may be modelled and surface pressure data for each filament may be generated for understanding the interaction among themselves. The algorithm used in the present study may be modified further for more accurate results. The major outcome of this study is to understand the hydrodynamic interaction between solid structure and fluid. Another major contribution of this study is the mathematical modelling of self-propelled free-swimming flexible filament and incorporating the calculation of surface pressure data at desired Lagrangian nodes, that may be significantly useful for future soft robotics applications. 60 APPENDIX 61 APPENDIX Derivation of Poisson equation for ๐‘ป Eq. (7) can be derived from Eqs. (2) and (5) as following, ๐œ•๐‘ฟ ๐œ•๐‘ฟ ( โ‹… )=1, ๐œ•๐‘  ๐œ•๐‘  Differentiating twice with respect to ๐‘ก, we get, ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ( โ‹… )=0, ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  Multiplying both sides with constant mass ratio ๐›ฝ, ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ โ‡’๐›ฝ ( โ‹… )=0 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•๐‘ฟ โ‡’๐›ฝ [ ( โ‹… )] = 0 ๐œ•๐‘ก ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ โ‡’ 2๐›ฝ [ โ‹… ( )] = 0 ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ ๐œ• ๐œ•๐‘ฟ โ‡’ 2๐›ฝ [ โ‹… 2 ( ) + ( ) โ‹… ( ) ] = 0 ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ โ‡’2 โ‹… (๐›ฝ 2 ) + 2๐›ฝ โ‹… =0 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• ๐œ• ๐œ•๐‘ฟ ๐œ•2 ๐œ•2๐‘ฟ ๐’ˆ ๐œ•2 ๐‘ฟ ๐œ•2๐‘ฟ โ‡’2 โ‹… ( (๐‘‡ ) โˆ’ 2 (๐›พ 2 ) + ๐›ฝ๐น๐‘Ÿ โˆ’ ๐‘ญ ) + 2๐›ฝ โ‹… =0 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐‘” ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐œ• 2 ๐œ•2๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐’ˆ ๐œ•๐‘ฟ ๐œ• โ‡’2 โ‹… 2 (๐‘‡ ) โˆ’ 2 โ‹… ( 2 (๐›พ 2 )) + 2 โ‹… (๐›ฝ๐น๐‘Ÿ ) โˆ’ 2 โ‹… ( ๐‘ญ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐‘” ๐œ•๐‘  ๐œ•๐‘  ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ + 2๐›ฝ โ‹… =0 ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  62 ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ•2 ๐‘ฟ ๐œ•2๐‘ฟ ๐œ•๐‘ฟ ๐œ• โ‡’2 โ‹… 2 (๐‘‡ ) + 2๐›ฝ โ‹… +2 โ‹… ( ๐‘ญ โˆ’ ๐‘ญ) = 0 , ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐‘ ๐œ•2 ๐œ•2๐‘ฟ ๐œ•2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ where ๐‘ญ๐‘ = โˆ’ ๐œ•๐‘  2 (๐›พ ๐œ•๐‘  2 ). Since, ๐œ•๐‘ก 2 ( ๐œ•๐‘  โ‹… ๐œ•๐‘  ) = 0, we can write, ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ•2๐‘ฟ ๐œ•2๐‘ฟ ๐œ•๐‘ฟ ๐œ• ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ 2 โ‹… 2 (๐‘‡ ) + 2๐›ฝ โ‹… +2 โ‹… ( ๐‘ญ๐‘ โˆ’ ๐‘ญ) = 2 ( โ‹… ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ•๐‘ฟ ๐œ• 1 ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ โ‡’ โ‹… 2 (๐‘‡ ) + ๐›ฝ โ‹… + โ‹… ( ๐‘ญ๐‘ โˆ’ ๐‘ญ) = ( โ‹… ) ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  2 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ฟ ๐œ• 2 ๐œ•๐‘ฟ 1 ๐œ• 2 ๐œ•๐‘ฟ ๐œ•๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ• 2 ๐‘ฟ ๐œ•๐‘ฟ ๐œ• โ‡’ โ‹… 2 (๐‘‡ ) = ( โ‹… ) โˆ’ ๐›ฝ โ‹… โˆ’ โ‹… ( ๐‘ญ โˆ’ ๐‘ญ) , ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  2 ๐œ•๐‘ก 2 ๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘ก๐œ•๐‘  ๐œ•๐‘  ๐œ•๐‘  ๐‘ here the first term on the RHS is zero theoretically, but to correct numerical errors of the inextensibility constraint, this term is not dropped. 63 REFERENCES 64 REFERENCES [1] Lighthill, M. (1960). Note on the swimming of slender fish. Journal of Fluid Mechanics, 9(2), 305-317. doi:10.1017/S0022112060001110. [2] E. Purcell, Life at low Reynolds number, American Journal of Physics 45, 3 (1977); doi: 10.1119/1.10903. [3] Altringham JD, Ellerby DJ. Fish swimming: patterns in muscle function. J Exp Biol. 1999 Dec; 202(Pt 23):3397-403. PMID: 10562522. [4] Eldredge JD. Numerical simulations of undulatory swimming at moderate Reynolds number. Bioinspir Biomim. 2006 Dec;1(4):S19-24. doi: 10.1088/1748-3182/1/4/S03. Epub 2006 Dec 22. PMID: 17671314. [5] F.E. Fish and G.V. Lauder, Passive and active flow control by swimming fishes and mammals. Annual Review of Fluid Mechanics 2006 38:1, 193-224. [6] J. J. Videler, Fish Swimming (Chapman and Hall, New York, 1993). [7] A. K. Brodsky, The Evolution of Insect Flight (Oxford University Press, Oxford, 1994). [8] S. A. Combes and T. L. Daniel, โ€œFlexural stiffness in insect wings. I. Scaling and the influence of wing venation,โ€ J. Exp. Biol. 206, 2979โ€“2987 (2003). [9] S. A. Combes and T. L. Daniel, โ€œFlexural stiffness in insect wings. II. Spatial distribution and dynamic wing bending,โ€ J. Exp. Biol. 206, 2989โ€“2997 (2003). [10] R. W. Blake, Fish Locomotion (Cambridge University Press, Cambridge, 1983). [11] R. J. Wootton, โ€œInvertebrate paraxial locomotory appendages: Design, deformation and control,โ€ J. Exp. Biol. 202, 3333โ€“3345 (1999). [12] S. Alben, P. G. Madden, and G. V. Lauder, โ€œThe mechanics of active fin-shape control in ray- finned fishes,โ€ J. R. Soc., Interface 4, 243โ€“256 (2007). [13] R. Godoy-Diana, J.-L. Aider, and J. E. Wesfreid, โ€œTransitions in the wake of a flapping foil,โ€ Phys. Rev. E 77, 016308 (2008). [14] M. S. Triantafyllou, G. S. Triantafyllou, and R. Gopalkrishnan, โ€œWake mechanics for thrust generation in oscillating foils,โ€ Phys. Fluids 3, 2835โ€“2837 (1991). [15] J. H. Buchholz and A. J. Smits, โ€œThe wake structure and thrust performance of a rigid low- aspect-ratio pitching panel,โ€ J. Fluid Mech. 603, 331โ€“365 (2008). [16] J. M. Anderson, K. Streitlien, D. S. Barrett, and M. S. Triantafyllou, โ€œOscillating foils of high propulsive efficiency,โ€ J. Fluid Mech. 360, 41โ€“72 (1998). 65 [17] N. Vandenberghe, S. Childress, and J. Zhang, โ€œOn unidirectional flight of a free flapping wing,โ€ Phys. Fluids 18, 014102 (2006). [18] T. Y.-T. Wu, โ€œSwimming of a waving plate,โ€ J. Fluid Mech. 10, 321โ€“344 (1961). [19] X.-Y. Lu and Q. Liao, โ€œDynamic responses of a two-dimensional flapping foil motion,โ€ Phys. Fluids 18, 098104 (2006). [20] Z. J. Wang, โ€œVortex shedding and frequency selection in flapping fight,โ€ J. Fluid Mech. 410, 323โ€“341 (2000). [21] Z. J. Wang, โ€œTwo dimensional mechanism for insect hovering,โ€ Phys. Rev. Lett. 85, 2216 (2000). [22] J. Carling, T. L. Williams, and G. Bowtell, โ€œSelf-propelled anguilliform swimming: Simultaneous solution of the two dimensional Navier-Stokes equations and Newtonโ€™s laws of motion,โ€ J. Exp. Biol. 201 (1998). [23] H. Dong, R. Mittal, and F. M. Najjar, โ€œWake topology and hydrodynamic performance of low- aspect-ratio flapping foils,โ€ J. Fluid Mech. 566, 309โ€“343 (2006). [24] Jingtao Ma, Zhen Wang, John Young, Joseph C.S. Lai, Yi Sui, Fang-Bao Tian, โ€œAn immersed boundary-lattice Boltzmann method for fluid-structure interaction problems involving viscoelastic fluids and complex geometries,โ€ Journal of Computational Physics, Volume 415, 2020, 109487, ISSN 0021-9991. [25] Ru-Nan Hua, Luoding Zhu, and Xi-Yun Lu,โ€œLocomotion of a flapping flexible plateโ€, Physics of Fluids 25, 121901 (2013). [26] Xiaojue Zhu, Guowei He, Xing Zhang, โ€œNumerical study on hydrodynamic effect of flexibility in a self-propelled plunging foil,โ€ Computers & Fluids, Volume 97, 2014, Pages 1-20, ISSN 0045- 7930. [27] Wei-Xi Huang, Soo Jai Shin, Hyung Jin Sung, โ€œSimulation of flexible filaments in a uniform flow by the immersed boundary method,โ€ Journal of Computational Physics, Volume 226, Issue 2, 2007, Pages 2206-2228, ISSN 0021-9991. [28] Farnell DJ, David T, Barton DC. Numerical model of self-propulsion in a fluid. J R Soc Interface. 2005 Mar 22;2(2):79-88. doi: 10.1098/rsif.2005.0027. PMID: 16849167; PMCID: PMC1578256. [29] J. Zhang, S. Childress, A. Libchaber, M. Shelley, Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind, Nature 408 (2000) 835โ€“839. [30] C.S. Peskin, The immersed boundary method, Acta Numer. (2002) 479โ€“517. [31] L. Zhu, C.S. Peskin, Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method, J. Comput. Phys. 179 (2002) 452โ€“468. [32] L. Zhu, C.S. Peskin, Interaction of two flapping filaments in a flowing soap film, Phys. Fluids 15 (7) (2003) 1954โ€“1960. 66 [33] Y. Kim, C.S. Peskin, Penalty immersed boundary method for an elastic boundary with mass, peprint, 2005. [34] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220โ€“ 252. [35] Z. Yu, A DLM/FD method for fluid/flexible-body interactions, J. Comput. Phys. 207 (2005) 1โ€“ 27. [36] Zhaowu Lin, Andrew Hess, Zhaosheng Yu, Shengqiang Cai, Tong Gao, โ€œA fluidโ€“structure interaction study of soft robotic swimmer using a fictitious domain/active-strain method,โ€ Journal of Computational Physics, Volume 376, 2019, Pages 1138-1155, ISSN 0021-9991. [37] D.J.J. Farnell, T. David, D.C. Barton, Numerical simulations of a filament in a flowing soap film, Int. J. Numer. Methods Fluids 44 (2004) 313โ€“330. [38] Fang-Bao Tian, Haoxiang Luo, Luoding Zhu, James C. Liao, Xi-Yun Lu, An efficient immersed boundary-lattice Boltzmann method for the hydrodynamic interaction of elastic filaments, Journal of Computational Physics, Volume 230, Issue 19, 2011, Pages 7266-7283, ISSN 0021-9991. [39] Shiyi Chen and Gary D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 1998 30:1, 329-364. [40] Cyrus K. Aidun and Jonathan R. Clausen, Lattice-Boltzmann method for complex flows, Annual Review of Fluid Mechanics 2010 42:1, 439-472. [41] O. Malaspinas, N. Fiรฉtier, M. Deville, Lattice Boltzmann method for the simulation of viscoelastic fluid flows, Journal of Non-Newtonian Fluid Mechanics, Volume 165, Issues 23โ€“24, 2010, Pages 1637-1653, ISSN 0377-0257. [42] Lincheng Xu, Fang-Bao Tian, John Young, Joseph C.S. Lai, A novel geometry-adaptive Cartesian grid based immersed boundaryโ€“lattice Boltzmann method for fluidโ€“structure interactions at moderate and high Reynolds numbers, Journal of Computational Physics, Volume 375, 2018, Pages 22-56, ISSN 0021-9991. [43] N. Goyal, J. Derksen, Direct simulations of spherical particles sedimenting in viscoelastic fluids, J. Non-Newton. Fluid Mech. 183 (2012) 1โ€“13. [44] Y.K. Lee, K.H. Ahn, A novel lattice Boltzmann method for the dynamics of rigid particles suspended in a viscoelastic medium, J. Non-Newton. Fluid Mech. 244 (2017) 75โ€“84. [45] Andrew K. Gunstensen, Daniel H. Rothman, Stรฉphane Zaleski, and Gianluigi Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43, 4320 โ€“ Published 1 April 1991. [46] T. Inamuro, T. Ogata, S. Tajima, and N. Konishi. 2004. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J. Comput. Phys. 198, 2 (10 August 2004), 628โ€“644. 67 [47] Misun Min and Taehun Lee, A spectral-element discontinuous Galerkin lattice Boltzmann method for nearly incompressible flows. J. Comput. Phys. (2011) 230, 1 (January, 2011), 245โ€“ 259. [48] Renwei Mei, Wei Shyy, Dazhi Yu, and Li-Shi Luo, Lattice Boltzmann method for 3-D flows with curved boundary, Journal of Computational Physics 161, 680โ€“699 (2000). [49] F.B. Tian, H. Luo, X.Y. Lu, Coupling modes of three filaments in side-by-side arrangement, Phys. Fluids 23 (2011) 111903. [50] F.B. Tian, H. Luo, L. Zhu, X.Y. Lu, Interaction between a flexible filament and a downstream rigid body, Phys. Rev. E 82 (2010) 026301. [51] F.-B. Tian, Role of mass on the stability of flag/flags in uniform flow, Appl. Phys. Lett. 103 (2013) 034101. [52] L. Zhu, A three-dimensional immersed boundary method for non-Newtonian fluids, Theor. Appl. Mech. Lett. 8 (2018) 193โ€“196. [53] Luoding Zhu, Guowei He, Shizhao Wang, Laura Miller, Xing Zhang, Qian You, Shiaofen Fang, An immersed boundary method based on the lattice Boltzmann approach in three dimensions, with application, Computers & Mathematics with Applications, Volume 61, Issue 12, 2011, Pages 3506-3518, ISSN 0898-1221. [54] D.-K. Sun, Y. Wang, A.-P. Dong, B.-D. Sun, A three-dimensional quantitative study on the hydrodynamic focusing of particles with the immersed boundary-lattice Boltzmann method, Int. J. Heat Mass Transf. 94 (2016) 306โ€“315. [55] Michael R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E 54, 5041 โ€“ Published 1 November 1996. [56] Jian Hao, Luoding Zhu, A lattice Boltzmann based implicit immersed boundary method for fluidโ€“structure interaction, Computers & Mathematics with Applications, Volume 59, Issue 1, 2010, Pages 185-193, ISSN 0898-1221. [57] Cheng, Yongguang & Zhang, Hui. (2010). Immersed boundary method and lattice Boltzmann method coupled FSI simulation of mitral leaflet flow. Computers & Fluids - COMPUT FLUIDS. 39. 871-881. 10.1016/j.compfluid.2010.01.003. [58] D.-K. Sun, Z. Bo, Numerical simulation of hydrodynamic focusing of particles in straight channel flows with the immersed boundary-lattice Boltzmann method, Int. J. Heat Mass Transf. 80 (2015) 139โ€“149. [59] T. Krรผger, F. Varnik, D. Raabe, Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method, Comput. Math. Appl. 61 (2011) 3485โ€“3505. [60] J. Ma, Z. Wang, Y. Sui, J. Young, J. Lai and F. Tian, An IB-LBM for FSI Problems Involving Viscoelastic Fluids and Complex Geometries, in: WCCM-ECCOMAS2020. 68 [61] Injae Lee, Haecheon Choi, A discrete-forcing immersed boundary method for the fluidโ€“ structure interaction of an elastic slender body, Journal of Computational Physics, Volume 280, 2015, Pages 529-546, ISSN 0021-9991. [62] M.D. de Tullio, G. Pascazio, A moving-least-squares immersed boundary method for simulating the fluidโ€“structure interaction of elastic bodies with arbitrary thickness, Journal of Computational Physics, Volume 325, 2016, Pages 201-225, ISSN 0021-9991. [63] Z.G. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems, J. Comput. Phys. 195 (2004) 602โ€“628. [64] Z.G. Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulations of particulate flows, J. Comput. Phys. 202 (2005) 20โ€“51. [65] C. Shu, N. Liu, Y.T. Chew, A novel immersed boundary velocity correction-lattice Boltzmann method and its application to simulate flow past a circular cylinder, J. Comput. Phys. 226 (2007) 1607โ€“1622. [66] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications, J. Comput. Phys. 228 (2009) 1963โ€“1979. [67] J. Wu, C. Shu, An improved immersed boundaryโ€“lattice Boltzmann method for simulating three-dimensional incompressible flows, J. Comput. Phys. 229 (2010) 5022โ€“5042. [68] [23] D. Yu, R. Mei, W. Shyy, A multi-block lattice Boltzmann method for viscous fluid flows, Int. J. Numer. Meth. Fluids 39 (2002) 99โ€“120. [69] Y. Peng, C. Shu, Y.T. Chew, X.D. Niu, X.Y. Lu, Application of multi-block approach in the immersed boundary-lattice Boltzmann method for viscous fluidflows, J. Comput. Phys. 218 (2006) 460โ€“478. [70] Y. Sui, Y.T. Chew, P. Roy, H.T. Low, A hybrid immersed-boundary and multi-block lattice Boltzmann method for simulating fluid and moving boundaries interactions, Int. J. Numer. Meth. Fluids 53 (2007) 1727โ€“1754. [71] J. Wu, C. Shu, An improved immersed boundaryโ€“lattice Boltzmann method for simulating three-dimensional incompressible flows, J. Comput. Phys.229 (2010) 5022โ€“5042. [72] Dai L, He G, Zhang X, Zhang X. 2018 Stable formations of self-propelled fish-like swimmers induced by hydrodynamic interactions. J. R. Soc. Interface 15: 20180490. [73] Ashraf I, Godoy-Diana R, Halloy J, Collignon B, and Thiria B. 2016 Synchronization and collective swimming patterns in fish (Hemigrammus bleheri). J. R. Soc. Interface 13, 1โ€“9. [74] Xu, Y.-Q., Jiang, Y.-Q., Wu, J., Sui, Y., & Tian, F.-B. (2018). Benchmark numerical solutions for two-dimensional fluidโ€“structure interaction involving large displacements with the deforming- spatial-domain/stabilized spaceโ€“time and immersed boundaryโ€“lattice Boltzmann methods. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 232(14), 2500โ€“2514. 69