A numerical scheme to simulate three-dimensional two-way fluid–structure interaction (twFSI) problems of flows around a flexible fine structure is developed in this study. The partitioned approach is employed to separately calculate fluid flows and structure motions by the lattice Boltzmann method (LBM) and the geometrically exact Cosserat rod model (CRM), respectively. The fluid–structure interactions are calculated by the simple explicit coupling scheme combined with the contact detection algorithm and the fluid–structure interface reconstruction scheme. The contact detection algorithm utilizing the bounding volume hierarchy is adopted to reduce the computing time of data communication between the fluid and the structure solvers, while the fluid–structure interface reconstruction scheme utilizes the level set method to represent the moving fluid–structure interfaces. The proposed LBM–CRM–twFSI scheme is successfully validated in two experimental benchmarks of a single flexible structure deformation in a wind tunnel. The results confirm that the present scheme accurately calculates the equilibrium state and the time-dependent oscillatory motions of the structures exposed to airflows. The errors of the representative rod position between the experimental and numerical results for both benchmarks are within 5%. These validations confirm the practicability of the presently developed LBM–CRM–twFSI scheme for motions of flexible fine structures in fluid flows.

Two-way fluid–structure interaction (twFSI) problems of the slender flexible structures are one of the essential research topics as they play important roles in biological flows and flows treated in ocean, civil, and textile engineering (Bourguet and Triantafyllou, 2016; Dey et al., 2018; Li et al., 2019; 2020; Mattis et al., 2015; Wu et al., 2014; and Zhu and Peskin, 2003). Since these twFSI problems are extremely complicated due to their non-linear nature and interactions between fluids and structures, conducting accurate experiments is very expensive and time-consuming. Thus, for understanding the underlying physics and designing engineering systems, various numerical frameworks capable of solving twFSI problems of the slender flexible structures have been investigated (Abdi et al., 2019; Tian et al., 2014; Wu et al., 2020; and 2014). These numerical methods are distinguished by the temporal and spatial coupling scheme of the structure motions and the surrounding fluid flows.

For the temporal coupling, two strategies (monolithic and partitioned algorithms) have been developed (Degroote et al., 2009; Ha et al., 2017). Compared with the monolithic algorithm that treats the fluid flow and structure motion by a single combined solver, the partitioned algorithm is employed in most cases of the twFSI problems owing to the simplicity and efficiency of the algorithm. The partitioned approach utilizes individual solvers for fluids and structures with an appropriate coupling scheme. Although data communication between the fluid and structure solvers is required, existing solvers can be adopted without modifications. As a fluid solver of the twFSI scheme, we focus on the lattice Boltzmann method (LBM), which is a computational fluid dynamics (CFD) technology based on the discretized Boltzmann equation. Compared with the conventional continuum-based CFD approaches, the LBM has the following advantages: simple fluid/structure representation of lattice nodes by binarized on/off treatments, easy implementation of no-slip boundary conditions of complex geometries, and easy parallelization using graphics processing units (GPUs) due to the discrete nature and locality of the computational algorithm (Agarwal et al., 2020; Lorenz et al., 2018; and Rinaldi et al., 2012). These features allow the LBM to be employed in simulations of the twFSI problems with complex geometry, flexible structures, or interacting boundaries. Owen et al. (2011) presented a serial numerical framework coupling the LBM and the discrete element method (DEM) for solving the particle suspension problems. Combined with the simple bounce-back method and the immersed boundary method (IBM), the LBM–DEM scheme demonstrated a stable and efficient simulation of multi-particle suspension. Favier et al. (2014) proposed a numerical framework based on the LBM and the IBM to calculate two-dimensional twFSI problems of moving structures. Their scheme was validated by the test cases of impulsively started flat plates, rigid particle sedimentation, and flapping filaments. They concluded that the validation cases showed good agreement with the reference numerical and experimental results in two dimensions. Li et al. (2020) presented the coupling of the LBM and the finite element method (FEM) combined with the IBM to solve twFSI problems. The LBM–FEM scheme was validated by the test cases of flows over a deformable plate and a cantilever flexible plate. They showed that numerical results agreed with the previous experimental and numerical investigations. As stated above, combined with various structure solvers suited for the flow cases, the LBM has been widely utilized as a fluid solver in twFSI problems. As for the structure solver of the twFSI problems of flexible fine structures, we focus on the geometrically exact Cosserat rod model (CRM) (Rubin, 2000), which approximates three-dimensional physical motions of the slender structures by one-dimensional mathematical models without losing physical accuracy. Although the DEM- and FEM-based methods require finer meshes to sufficiently resolve the diameter of slender flexible structures to handle large and complex deformations, the CRM does not require finer meshes by the reduction of the degrees of freedom (DOF) to handle large structural deformations or displacement. By the DOF, what we need is solving the linear and angular momentum balance for 1D discretized nodes and cylindrical segments. This process is much simpler than those required for the 3D FEM-based methods as reported in the literature (Gazzola et al., 2018; Spillmann and Teschner, 2007; and Tschisgale and Fröhlich, 2020). Accordingly, the total computational costs of the CRM become lower to simulate bending, twisting, shearing, and stretching of the structure under a wide range of boundary conditions (Cao et al., 2006; Cao and Tucker, 2008; Jung et al., 2012; and Spillmann and Teschner, 2007). Marheineke and Wegener (2011) presented the one-way fluid–structure interaction (owFSI) scheme composed of the Navier–Stokes equations and the CRM. The air drag model is valid for various Reynolds number ranges, and incident flow directions were formulated and introduced as the external force in the CRM. They demonstrated the numerical simulations of the melt-spinning process of nonwoven fabrics and compared with the experimental data. As a result, the probability distributions of the fibers were in good agreement. Tschisgale and Fröhlich (2020) proposed the twFSI scheme by coupling the Navier–Stokes equations and the CRM combined with a new semi-implicit coupling approach based on the IBM. Their proposed scheme was applied to the flow through an artificial aquatic canopy consisting of 800 flexible structures and demonstrated its stability and efficiency. As stated above, combined with the Navier–Stokes equations, the CRM has been applied in a structure solver for some owFSI and twFSI problems. Although the advantages of the LBM and the CRM are suitable for solving twFSI problems of the slender flexible structures, to the best of the authors’ knowledge, there has been no report in the literature on the twFSI scheme based on coupling the LBM and the CRM.

As for the spatial coupling, either a boundary fitted mesh or a spatially uniform grid has been applied to moving structures. The former approach requires mesh conformation along the fluid–structure interfaces or remeshing process to maintain mesh quality near the moving boundaries (Dukowicz and Kodis, 1987; Margolin and Shashkov, 2003). These mesh deformation or regeneration processes are time-consuming and can lead to errors due to transferring data from the previous to new meshes. Thus, a boundary fitted moving mesh technique is not practical to solve the twFSI problems with complex boundaries or large deformations. Recently, the latter approach utilizing a spatially uniform Eulerian grid to solve the fluid flow and a Lagrangian expression to solve the structure motions has been alternatively employed (Bavo et al., 2016). This technique has advantages of the simple algorithm, no grid generation for the fluid solver independent of complex geometries or large deformations of the fluid–structure interfaces, and high efficiency of the fluid solver owing to the regular grid. While various methods to represent the fluid–structure interface with the spatially uniform grid have been reported, such as volume-of-fluid (VOF) methods (Cerroni et al., 2018; Kassiotis et al., 2010; Patel et al., 2017; and Shen and Chan, 2008), phase field (PF) methods (Mokbel et al., 2018; Sun et al., 2014; and Wick, 2016), or IBM (Lim and Peskin, 2012; Suzuki and Inamuro, 2011; and Zhu et al., 2011), we focus on the level set method to represent the structure shapes. The level set method is a well-known method for interface capturing, which utilizes the level set function indicating normal distance from an interface between two phases. With mathematically described structure shapes, the level set method allows us to represent the fluid–structure interfaces in a spatially uniform Eulerian grid and to impose twFSI boundary conditions to interfaces by a simple explicit coupling scheme (Calderer et al., 2014; Legay et al., 2006).

Therefore, this study develops a numerical framework of the three-dimensional twFSI scheme by coupling the LBM and the CRM (LBM–CRM–twFSI) utilizing a simple partitioned explicit coupling scheme for the temporal coupling and the level set method for the spatial coupling. For demonstrating its practical performance, validations by experimental test cases are also carried out. The organization of this paper is as follows: Secs. II A and II B describe the fluid solver of the LBM with the bounce-back model and the structure solver of the CRM, respectively. Section III A presents the overall description of the fluid–structure coupling procedure for the proposed LBM–CRM–twFSI. Sections III B–III E are the descriptions of the individual principal schemes of the coupling procedure. Sections IV A and IV B give the verifications of the fluid solver of the LBM and the structure solver of the CRM, respectively. Section IV C validates the LBM–CRM–twFSI by comparing with two experimental benchmarks. Finally, Sec. V gives the concluding summary of the article.

Based on the discretized Boltzmann equation, the LBM calculates the time-evolution of the distribution function of the groups of fluid particles that stream and collide at lattice nodes (Frisch et al., 1986; Wolfram, 1986). The streaming process propagates the particles to adjacent lattice nodes, while the collision process redistributes the particles at each lattice node toward a local equilibrium state (Aidun and Clausen, 2010). This evolution of the distribution function is discretized in time and space as

fx+eαδt,t+δt=fx,t+Ωx,t,
(1)

where f is the density distribution function, Ω is the collision operator, x is the lattice position, eα is the discrete velocity, subscript α=0,1,,Q1 denotes the finite velocity component of the distribution function, t is the time, and δt is the time step. Here, the notation for the column vector, e.g., fx,tRQ is equivalent to

fx,t=f0x,t,f1x,t,,fQ1x,tT.
(2)

Equation (1) describes that the fluid particles move with their discrete velocity eα to the neighboring lattice node x + eαδt at the next time step t + δt, followed by the collision process under the control of the collision operator Ω. In this study, the multiple-relaxation-time (MRT) model (D’Humières et al., 2002; Suga et al., 2015) is employed, which is expressed as

Ωx,t=M1Ŝmx,tmeqx,t,
(3)

where M is the transformation matrix, Ŝ is the relaxation matrix, mx,t is the moment, and meqx,t is the equilibrium moment. The transformation matrix MRQ×Q linearly transforms the density distribution function f to the moment as m=Mf and as meq=Mfeq, where feq is the equilibrium distribution function typically determined by the density ρf and the velocity uf of the fluid. The equilibrium distribution function feq is given as

fαeqx,t=wαρf1+eαufcs2+eαuf22cs4ufuf2cs2,
(4)

where wα is the weight coefficient and cs is the sound speed that is equivalent to c/3 for isothermal flow with the particle velocity c = Δ/δt and the lattice spacing Δ. The weight coefficient wα and the sound speed cs are associated with the lattice discretizing models (Qian, 1993). In this study, the D3Q27 model is employed, and its discrete velocity vector eα is

e0e1e2e3e4e5e6e7e8e9e10e11e12e13e14e15e16e17e18e19e20e21e22e23e24e25e26/c=000100010100010001001110110110110101011101011101011101011111111111111111111111111,
(5)

and Table I lists the weight coefficient wα.

TABLE I.

Weight coefficient of the D3Q27 discrete velocity model.

wα
 8/27 α=0  
 2/27 α=1,,6  
 1/54 α=7,,18  
 1/216 α=19,,26  
wα
 8/27 α=0  
 2/27 α=1,,6  
 1/54 α=7,,18  
 1/216 α=19,,26  

The transformation matrix M of the D3Q27 model is set the same as that in Suga et al. (2015). The collision matrix ŜRQ×Q, which is a diagonal matrix, is given as

Ŝ=diag0,0,0,0,s4,s5,s5,s7,s7,s7,s10,s10,s10,s13,s13,s13,×s16,s17,s18,s18,s20,s20,s20,s23,s23,s23,s26,
(6)

where sα are chosen as s4 = 1.54, s5 = s7, s10 = 1.5, s13 = 1.83, s16 = 1.4, s17 = 1.61, s18 = s20 = 1.98, and s23 = s26 = 1.74. The fluid density ρf, velocity uf, and pressure p are obtained from the conservation law and the equation of state as follows:

ρf=αfa,
(7)
ρfuf=αfaeα,
(8)
p=ρfcs2.
(9)

The fluid kinematic viscosity νf is given as

νf=cs2τs12δt,
(10)

where τs is the hydrodynamic relaxation time equivalent to 1/s5 and 1/s7.

Throughout this study, we utilize the linear interpolated bounce-back (IPBB) method as a non-slip boundary condition of the LBM to treat the curved and moving wall boundaries (Bouzidi et al., 2001; Lallemand and Luo, 2003; and Pan et al., 2006). Compared to the conventional half-way bounce-back (HWBB) method, the IPBB method can be applicable for more accurate calculations to present curved boundaries with second order accuracy in the case with mathematically described fluid–structure interfaces (De Rosis et al., 2014; Li et al., 2004; and Obrecht et al., 2013). Combined with the fluid–structure boundary reconstruction (fsBR) scheme described in Sec. III E, this intrinsic feature of the IPBB method provides the advantage in simulating the motion of the fine structures with a smooth circular cross-sectional shape.

Quasi-one-dimensional structures usually referred to as beams, cantilevers, or rods have a much smaller radius than their longitudinal length at any cross sections. In this study, we focus on the latter term “rod,” which has a long slender cylindrical shape. In principle, such characteristics of geometric slenderness provide significant simplification of the mathematical model to present the rod motion without losing physical accuracy. One of the typical rod models is the CRM, which approximates three-dimensional physical motions of the slender structures by one-dimensional mathematical models. Since the CRM can handle a large structural deformation and take account of bending, twisting, shearing, and stretching of the rod under a wide range of boundary conditions, we utilize this model as the structure solver in the present study.

The basic theory of the CRM is illustrated in Fig. 1. The three orthonormal unit vectors e1,e2,e3 form a fixed right-handed basis. A three-dimensional rod in the CRM is described as a space curve corresponding to its center line positions rs,tR3 and a cross-sectional orientation characterized by a director frame Rs,tSO3, which is equivalent to the three orthonormal directional unit vectors R=d1,d2,d3. Note that SO3 indicates the 3D rotation group about the origin of R3, often utilized in mechanics and geometry fields, and the matrix R translates a vector a in the Eulerian frame to a vector aL in the Lagrangian frame via aL=Ra and, conversely, a=RTaL. Here, s0,LR is the center line arc-length coordinate in the current configuration, L is the characteristic length, and tR+ is the time. The CRM theory introduces the Lagrangian frame to handle rotations of rod cross sections relative to each other. Combining the Lagrangian frame with the 3D rotation group at every discretized cross section of the rod, the CRM can capture the complex deformations of the filament (i.e., bend, twist, shear, and stretch) as described in the literature (Cosserat and Cosserat, 1909; Love, 1892). The rod angular velocity ω and the generalized curvature κ are defined in terms of the rotation matrix R as

ω=vecRtTR,
(11)
κ=vecRsTR,
(12)

where vecV is the 3-vector associated with the skew-symmetric matrix V. Denoting ŝ is the center line arc-length coordinate in the reference configuration, the stretching of the rod relative to the reference configuration is expressed by the scalar ratio parameter eŝ,t=ds/dŝ (hatted symbol indicates the reference configuration, hereafter). As d3 is the normal to the cross section and d1 and d2 is a pair of director fields spanning the cross-sectional planes along the center line in the unsheared state, the shearing of the rod is represented by the shifting of three orthonormal directional unit vectors d1,d2,d3 away from the unit tangential vector t, and thus, d3t. Here, the vector σL is defined as σL=Retd3 to determine the axial and shear strains. Note that the conditions e = 1 and d3 = t are equivalent to the Kirchhoff theory for inextensible and unshearable rods (Coleman et al., 1993; Duričković et al., 2013; Dill, 1906; and Meier et al., 2019). By assuming a linear elastic material, the constitutive laws to formulate the internal forces and moments generated by bending/twisting and shearing/stretching are given as

τL=BκLκLt=0,
(13)
nL=SσLσLt=0,
(14)

where τL is the internal moment, nL is the internal force, BR3×3=diagB1,B2,B3 is the bending/twisting stiffness matrix, and SR3×3=diagS1,S2,S3 is the shearing/stretching stiffness matrix. Here, B1 = EI1 is the bending rigidity about d1, B2 = EI2 is the bending rigidity about d2, B3 = GI3 is the twisting rigidity about d3, S1 = αcGA is the shearing rigidity along d1, S2 = αcGA is the shearing rigidity along d2, and S3 = EA is the stretching rigidity along d3 with Young’s modulus E, the shear modulus G, the cross-sectional area A, the constant shear factor αc, and the second moment of area IR3×3=diagI1,I2,I3. The initial strain state is expressed by κLt=0 and σLt=0. From the definitions above, the material properties of the rod are represented by Young’s modulus E and the shear modulus G, while its geometric properties are represented by the cross-sectional area A, the second moment of area I, and the constant shear factor αc. Finally, the formulations of the linear and angular momentum balance are given as

dm2rt2=ŝRTŜσLedŝ+F,
(15)
dĴeωLt=ŝB̂κ̂Le3dŝ+κ̂L×B̂κ̂Le3dŝ+Rt×ŜσLdŝ+dĴωLe×ωL+dĴωLe2et+CL,
(16)

where dm=ρsÂdŝ=ρsAds is the infinitesimal mass element, dĴ=ρsÎdŝ is the infinitesimal mass second moment of the area with the rod density ρs, F=efdŝ is the external force, and CL=ecLdŝ is the external moment with the force and moment line density f and cL. The discretized model of a rod segment in the CRM is illustrated in Fig. 2. The continuous rod is discretized into n+1 nodes of center line positions and n connecting cylindrical segments and thus is represented by a set of nodes ritR3,i1,n+1 and a set of director frames RitR3×3=d1i,d2i,d3i,i1,n. Each node is represented by the position vector ri and the mass mi, while each connecting cylindrical segment is represented by the position vectors at both ends ri,ri+1 and the current length li. Therefore, the displacement and rotation of the discretized rod are determined by the internal and external forces applied to the nodes and internal and external moments to the segments, respectively. See Cao and Tucker (2008), Gazzola et al. (2018), Goicoechea et al. (2019), Spillmann and Teschner (2007), and Zhang et al. (2019) for more details.

FIG. 1.

Schematic representation of the basic theory of the CRM.

FIG. 1.

Schematic representation of the basic theory of the CRM.

Close modal
FIG. 2.

Schematic representation of the discretized rod segment in the CRM.

FIG. 2.

Schematic representation of the discretized rod segment in the CRM.

Close modal

In this study, we utilize the partitioned approach that handles the fluid and the structure solvers separately and then couples them by an appropriate coupling algorithm. The procedures of the LBM–CRM–twFSI coupling scheme are based on a simple weak (explicit) coupling scheme (Farhat and Lesoinne, 2000; Felippa et al., 2001; Fernández, 2011; and Park et al., 1977), which is summarized as follows:

  1. A three-dimensional flow field including the structures is computed by the LBM.

  2. The hydrodynamic forces are calculated by the Galilean invariant momentum exchange method (GI-MEM).

  3. The hydrodynamic force information is transferred to the structure solver to calculate the external forces and moments. They are subsequently applied to the corresponding nodes and segments of the discretized rods via the contact detection algorithm.

  4. A three-dimensional motion of the rods is computed by the CRM.

  5. The geometry and velocity information of the rods are transferred back to the fluid solver. They are subsequently assigned to the corresponding lattice nodes to update the structure information in the LBM via the contact detection algorithm.

  6. The refilling algorithm and the fsBR scheme are applied to reconstruct the flow field and the fluid–structure interface, respectively.

The above procedures are repeated until the simulation reaches the pre-defined computational time. Figure 3 shows the schematic illustration of the LBM–CRM–twFSI coupling scheme. The details of the principal schemes in coupling procedures are described in Secs. III B–III E.

FIG. 3.

Diagram of the LBM–CRM–twFSI coupling procedures.

FIG. 3.

Diagram of the LBM–CRM–twFSI coupling procedures.

Close modal

Evaluating hydrodynamic forces is the important process in the twFSI problems. In this study, the GI-MEM is utilized to compute the hydrodynamic forces, which is simple, accurate, and independent of boundary interface geometries (Tan et al., 2018; Tao et al., 2016; Wen et al., 2015; and 2014). According to the theorem of momentum, the hydrodynamic force is represented by the momentum-exchange value during the bounce-back process on the fluid–structure interface. Thus, the GI-MEM is formulated as

Fαx,t=eαuWfαx,teαuWfαx,t+δt,
(17)

where Fα is the hydrodynamic force and subscripts α and α′ represent directions opposite to each other, which are defined as eα=eα. This indicates that the mass component fαx,t has the velocity eαuW relative to the boundary and increases a momentum eαuWfαx,t to the boundary, while the mass component fαx,t+δt has the relative velocity eαuW and decreases a momentum eαuWfαx,t+δt from the boundary.

To compute three-dimensional motions of the rods, the external forces and moments applied to nodes and segments of the discretized rods are given as

Fi=k=1NbFfxb,t,
(18)
Ci=k=1NbxbXG×Ffxb,t,
(19)

where Fi is the external force, Ci is the external moments, Ff is the hydrodynamic force calculated by the summation of the GI-MEM forces Fα, xb is the position of the fluid–structure interface lattice nodes, XG is the position of the discretized rod centroids, and Nb is the number of fluid–structure interface lattice nodes. The contact detection algorithm (Schneider and Eberly, 2003; van den Bergen, 1997) is utilized to accelerate connecting the discretized rod segment and the corresponding fluid–structure interface lattice nodes, which is summarized as follows:

  1. Assigning a Morton code for each lattice node and then sorting the Morton codes along the Z-order curve, which is the simple space-filling curve.

  2. The bounding volume hierarchy (BVH) is built using the parallel BVH generation algorithm proposed by Karras (2012). The BVH with n leaf nodes corresponding to the lattice nodes contains n − 1 internal nodes. The leaf nodes represent a set of n keys expressed as bit strings and each internal node partitions its keys according to their first differing bit and thus corresponds to the longest common prefix shared by the keys in its respective sub-tree.

  3. Inclusion of each lattice node in the discretized rod segment is detected in two steps.

    • Step 1: Generating axis-aligned bounding boxes (AABBs) that cover the discretized rod segments. Subsequently, contact detection between each AABB of the discretized rod segment and the AABBs corresponding to the internal and leaf nodes is recursively explored according to the generated BVH (Gao et al., 2014; Zhang and Kim, 2007). In this step, non-contact lattice nodes are roughly excluded.

    • Step 2: The remaining leaf nodes undergo further detailed contact detection. Here, we consider a lattice node as the centroid of the oriented bounding box (OBB) and a discretized rod segment as the oriented cylinder (OC). In this step, contact detection between a vertical line in the OC axial direction from each lattice node and the OC axis is first explored. Then, inclusion of the lattice node in the discretized rod segment is detected if the vertical line length is shorter than the OC radius (Ketchel and Larochelle, 2006). Note that the lattice spacing is added to the OC axis length and radius above to include the fluid–structure interface lattice nodes. The two-dimensional schematic illustrations of Steps 1 and 2 are shown in Figs. 4(a) and 4(b), respectively.

FIG. 4.

Illustrations of the two-step contact detection between the fluid–structure interface lattice nodes and the discretized rod segment: Contact detection algorithm of (a) AABB–AABB and (b) OBB–OC.

FIG. 4.

Illustrations of the two-step contact detection between the fluid–structure interface lattice nodes and the discretized rod segment: Contact detection algorithm of (a) AABB–AABB and (b) OBB–OC.

Close modal

The same algorithm is utilized for assigning the shape and velocity information of the rods to the corresponding lattice nodes to update the structure information in the LBM.

Due to the rod motion, some lattice nodes initially covered by the rod are converted from the solid-phase to the fluid-phase. The distribution functions at these flesh nodes are unknown and need to be initialized by a refilling algorithm. While a variety of refilling algorithms have been developed in previous studies (Caiazzo, 2008; He et al., 2020; Krithivasan et al., 2014; Peng et al., 2016; and Tao et al., 2016), we utilize the method proposed by Gao et al. (2013). This method reconstructs flesh nodes by the equilibrium distribution with non-equilibrium corrections, which is formulated as

fαx,t=fαeqx;ρ̄f,uW+fαneqx+cnδt,t,
(20)

where fαneq is the non-equilibrium distribution function, ρ̄f is the averaged fluid density of available neighboring fluid lattice nodes, and cn is the extrapolation direction. The equilibrium distribution is calculated in terms of the local wall velocity uW and the averaged fluid density ρ̄f of available neighboring fluid lattice nodes. The non-equilibrium correction is defined by copying the non-equilibrium distribution function from the neighboring fluid lattice nodes along the extrapolation direction cN. Here, cN is a specified discrete velocity along which direction the quantity nW · cN takes the maximum value, where nW is the unit normal vector pointing outwards of the moving boundary from the local wall surface.

The fluid and structure solvers are required to be spatially coupled at the fluid–structure interface in the twFSI problems. In this study, we utilize the level set functions to represent the rod shapes. Since the rod shape information received from the CRM is assigned to the regularly structured lattice nodes via the contact detection algorithm, the transferred rod shape in the LBM is represented by voxels and its surface is in a zig-zag manner. Hence, to recover the original round rod shape, the fsBR scheme developed by Suga et al. (2009) is employed. This method successfully recovered original round fiber shapes of nonwoven fabrics from zig-zag patterned computed tomography data by utilizing the level set functions (Ando et al., 2020). In this study, the level set function ϕ on an Eulerian grid is ϕ = 0, ϕ > 0, and ϕ < 0 to represent the fluid–structure interface, the structure, and fluid regions, respectively. Then, the fsBR scheme applies the level set re-initialization process to reconstruct the original smooth fluid–structure interface. The re-initialization is performed by solving the following partial differential equation (Sussman et al., 1996):

ϕtinit+nbϕ=sgnϕ0,
(21)

where tinit is the virtual time step during the re-initialization process and nb is the unit normal to the fluid–structure interface and given as

nb=sgnϕ0ϕϕ,
(22)

with initial conditions ϕx,0=ϕ0x.

In this section, the LBM combined with the GI-MEM and the fsBR scheme is first verified in terms of the lattice resolution by calculating the flow past a rigid circular cylinder. The flow behavior depends on the Reynolds number ReD = UinD/νf, where Uin is the inlet velocity, D is the cylinder diameter, and νf is the fluid kinetic viscosity. The present simulations are performed at Re = 40 and 200, and the results are compared with the literature values. The computational domain is 40D × 40D with the cylinder centered at 20D,16D to minimize the effect of the free-outflow boundary on the numerical results. The lattice resolutions are set to be D/Δ = 6, 10, 20, and 40. The cylinder shape data required for the IPBB method are obtained by the fsBR scheme. Figures 5(a) and 5(b) show the representative velocity contour maps for the cases of Re = 40 and 200, respectively. As shown in Fig. 5, the wake region is observed for the case of Re = 40, while the vortex shedding is observed for the case of Re = 200.

FIG. 5.

Representative velocity contour maps for the cases of (a) Re = 40 and (b) Re = 200, respectively. The velocity magnitude is normalized by Uin.

FIG. 5.

Representative velocity contour maps for the cases of (a) Re = 40 and (b) Re = 200, respectively. The velocity magnitude is normalized by Uin.

Close modal

The drag and the lift coefficients are calculated in terms of the hydrodynamic forces obtained by the GI-MEM, which are given as

CD=FD12ρfUin2D,
(23)
CL=FL12ρfUin2D,
(24)

where FD is the drag force and FL is the lift force. Figures 6(a) and 6(b) show the drag coefficient against the lattice resolution for the cases of Re = 40 and 200, respectively. As shown in Fig. 6, the hydrodynamic force asymptotically approaches the limiting value with the lattice resolution of finer than 20 for the cylinder diameter.

FIG. 6.

Variation of the drag coefficient against the lattice resolution at (a) Re = 40 and (b) Re = 200, respectively. The error bar shows the oscillation amplitude of the drag coefficient generated by the vortex shedding.

FIG. 6.

Variation of the drag coefficient against the lattice resolution at (a) Re = 40 and (b) Re = 200, respectively. The error bar shows the oscillation amplitude of the drag coefficient generated by the vortex shedding.

Close modal

The drag and the lift coefficients calculated with the lattice resolution of D/Δ = 20 are compared with the data in the literature (Dennis and Chang, 1970; Liu et al., 1998; Lv et al., 2006; Niu et al., 2006; 2003; and Tai and Zhao, 2003) as shown in Table II.

TABLE II.

Comparison of the drag and the lift coefficients.

ReResearchersCDCL
40 Dennis and Chang (1970)  1.522 ⋯ 
Niu et al. (2003)  1.574 ⋯ 
Niu et al. (2006)  1.589 ⋯ 
Present 1.557 ⋯ 
200 Lv et al. (2006)  1.38 ± 0.046 ±0.65 
Tai and Zhao (2003)  1.31 ± 0.041 ±0.64 
Liu et al. (1998)  1.31 ± 0.049 ±0.69 
Present 1.35 ± 0.043 ±0.66 
ReResearchersCDCL
40 Dennis and Chang (1970)  1.522 ⋯ 
Niu et al. (2003)  1.574 ⋯ 
Niu et al. (2006)  1.589 ⋯ 
Present 1.557 ⋯ 
200 Lv et al. (2006)  1.38 ± 0.046 ±0.65 
Tai and Zhao (2003)  1.31 ± 0.041 ±0.64 
Liu et al. (1998)  1.31 ± 0.049 ±0.69 
Present 1.35 ± 0.043 ±0.66 

The present results are in good agreement with the literature values, and therefore, we confirm that the LBM combined with the GI-MEM and the fsBR scheme is correctly implemented and the lattice resolution of D/Δ = 20 is accurate enough to evaluate the hydrodynamic force on the cylindrical structure by the present method.

In this section, the CRM is verified against the Timoshenko cantilever beam (TCB) theory. This model takes account of bending elasticity, rotary inertia, and shear deformations. A cantilever beam with circular cross sections clamped at one end ŝ=0 and subjected to the downward uniform force FU is considered, as illustrated in Fig. 7. The theoretical solution of the displacement y is analytically expressed as

y=FUŝ24EÎŝ3+4L̂ŝ26L̂2ŝFUŝ2αcÂG2L̂ŝ,
(25)

where L̂ is the length of the rod, Âis the constant cross-sectional area, Î is the second moment of area, E is the Young’s modulus, G is the shear modulus, and αc = 4/3 is the Timoshenko shear factor for the circular cross section. Considering the assumption that the TCB theory is for small deformations, the horizontal coordinate x is approximated by the arc-length ŝ. See Hutchinson (2001), Levinson (1981), and Timoshenko (1921) for more details.

FIG. 7.

Schematic representation of the Timoshenko cantilever beam.

FIG. 7.

Schematic representation of the Timoshenko cantilever beam.

Close modal

Figures 8(a) and 8(b) show comparison of the cantilever beam displacements calculated by the CRM and the theoretical solution for the two test cases. Computational conditions are shown in Table III. Here, the material and geometrical properties of the rods in the two test cases are the same as those in the literature (Jung et al., 2018; Silva-Leon et al., 2018) to which we refer for the validation of the LBM–CRM–twFSI later in Sec. IV C. The number of rod discretization nodes is set to be 30. The time step Δt is defined in terms of the spatial discretization as Δt = 10−2Δl, where Δl is the length of the discretized rod segment, as proposed by Gazzola et al. (2018). As shown in Fig. 8, the present results are in good agreement with the theoretical ones, and therefore, we confirm that the CRM is correctly implemented.

FIG. 8.

Comparison between the CRM numerical displacements (black dotted lines) and Timoshenko analytical displacements (colored solid lines) of the cantilever beam for (a) case 1 and (b) case 2, respectively.

FIG. 8.

Comparison between the CRM numerical displacements (black dotted lines) and Timoshenko analytical displacements (colored solid lines) of the cantilever beam for (a) case 1 and (b) case 2, respectively.

Close modal
TABLE III.

Computational conditions for the Timoshenko cantilever beam.

CaseDensity, ρs (kg/m3)Young’s modulus, E (Pa)Poisson’s ratio, ν[−]Length, L̂ (mm)Diameter, D (mm)Force, FU (N)
1-1 1000 2.25 × 105 0.01 25 0.4 0.5 × 10−6 
1-2 1.0 × 10−6 
1-3 1.5 × 10−6 
2-1 950 1.78 × 106 0.48 70 3.0 0.5 × 10−3 
2-2 1.0 × 10−3 
2-3 1.5 × 10−3 
CaseDensity, ρs (kg/m3)Young’s modulus, E (Pa)Poisson’s ratio, ν[−]Length, L̂ (mm)Diameter, D (mm)Force, FU (N)
1-1 1000 2.25 × 105 0.01 25 0.4 0.5 × 10−6 
1-2 1.0 × 10−6 
1-3 1.5 × 10−6 
2-1 950 1.78 × 106 0.48 70 3.0 0.5 × 10−3 
2-2 1.0 × 10−3 
2-3 1.5 × 10−3 

The proposed LBM–CRM–twFSI scheme is validated in two experimental benchmarks of a single flexible rod deformation in a wind tunnel. The first benchmark is the wind tunnel experiment of flexible rods exposed to steady uniform airflows reported by Silva-Leon et al. (2018). The experimental settings and the computational domain are illustrated in Fig. 9.

FIG. 9.

Schematic representation of the experimental setting of Silva-Leon et al. (2018).

FIG. 9.

Schematic representation of the experimental setting of Silva-Leon et al. (2018).

Close modal

The experiments were performed in a wind tunnel of square cross section with a height and width of 350 mm. The upper end of the silicone rubber flexible rod was introduced into a straight stainless steel support tube for about 5 mm and then fixed with glue cement to enforce a cantilever boundary condition. Under the influence of steady uniform airflows and gravity, the flexible rod was deflected into a static equilibrium inclined position depending on the wind velocity and the material and geometric properties of the rod. The equilibrium rod positions and the hydrodynamic forces generated by airflows were observed by the standard image processing of the recorded videos and a non-contact optical technique, respectively. Two experimental conditions in the literature (named cases A and B in this study) are chosen to be compared with our numerical results. The fluid and structure properties for both cases are listed in Tables IV and V.

TABLE IV.

Fluid properties (air).

Density,Kinetic viscosity,Inlet velocity,Reynolds number
Caseρf (kg/m3)νf (m2/s)Uin (m/s)Re
1.21 1.51 × 10−5 1.1 29.1 
2.1 55.6 
Density,Kinetic viscosity,Inlet velocity,Reynolds number
Caseρf (kg/m3)νf (m2/s)Uin (m/s)Re
1.21 1.51 × 10−5 1.1 29.1 
2.1 55.6 
TABLE V.

Structure properties (silicone rubber).

Density,Young’s modulus, Poisson’s ratio, Length, Diameter,
Caseρs (kg/m3)E (Pa)ν[−] L̂ (mm)D (mm)
1000 2.25 × 105 0.01 25 0.4 
Density,Young’s modulus, Poisson’s ratio, Length, Diameter,
Caseρs (kg/m3)E (Pa)ν[−] L̂ (mm)D (mm)
1000 2.25 × 105 0.01 25 0.4 

The Reynolds numbers are based on the inlet velocity Uin, the rod diameter D, and the kinetic viscosity of air νf. We reproduce the above experimental conditions by the LBM–CRM–twFSI scheme. The detailed computational configuration is described in Fig. 10. To save the computational cost, the computational domain is set to be a part of the wind tunnel as Lx × Ly × Lz = 100 × 650 × 740 lattice nodes with the lattice resolution of 0.04 mm/lattice spacing, which is equivalent to ten lattice nodes in the rod diameter. Note that the calculation error of the hydrodynamic force for this lattice resolution is estimated to be about 2.8% referring to the verification results shown in Fig. 6(a) in Sec. IV A. The upper end of the rod is positioned at 50,650,50 with the free-rotation boundary condition. As stated in the literature, the effect of the rotational stiffness at the fixed point is ignored in the simulation model due to the high flexibility of the rod (Silva-Leon et al., 2018). The initial position of the rod is inclined to shorten the time to approach the equilibrium state. The constant inlet velocity is applied, while the free-outflow condition is applied to the outlet. The slip flow condition is applied to the four lateral boundaries. A constant time step is 2.18 × 10−6 s for case A and 1.14 × 10−6 s for case B. The gravitational acceleration g acts opposite to the Y-direction. The rod is discretized by 30 nodes, and the shear factor is set to be αc = 4/3.

FIG. 10.

Schematic representation of the computational domain at the initial state for the cases of Silva-Leon et al. (2018).

FIG. 10.

Schematic representation of the computational domain at the initial state for the cases of Silva-Leon et al. (2018).

Close modal

The rod shape and the fluid velocities around the rod at the time steps (TSs) of 0, 38 000, and 570 000 for case A and 0, 76 000, and 608 000 for case B are shown in Figs. 11(a) and 11(b), respectively. Depending on the hydrodynamic and gravitational force balance, the rod shows inclined oscillatory motion at the transition phase and then gradually approaches the equilibrium state with the decrement of the oscillation amplitude.

FIG. 11.

Time evolutions of the rod shape and the fluid velocity profile for (a) case A (Re = 29.1) and (b) case B (Re = 55.6). The velocity magnitude is normalized by Uin.

FIG. 11.

Time evolutions of the rod shape and the fluid velocity profile for (a) case A (Re = 29.1) and (b) case B (Re = 55.6). The velocity magnitude is normalized by Uin.

Close modal

Figure 12 shows the equilibrium rod position obtained by the experiments and simulations, respectively. The equilibrium inclined angles of the rod obtained from the experiments and simulations are 26.5° and 25.9° for case A and 49.6° and 48.4° for case B, respectively. The hydrodynamic forces acting on the rod are also compared between the experimental and numerical results, resulting in 10.2 µN and 11.3 µN for case A and 19.0 µN and 19.3 µN for case B, respectively. Compared with the experimental results, the hydrodynamic forces though are larger in the simulation, the rod is slightly less deflected. This conflict is considered to be from the measurement uncertainties of the material properties provided in the literature, such as Young’s modulus and the Poisson ratio. Another reason can be related to the isotropic linear elastic constitutive laws, the free-rotation boundary condition applied in this simulation, or the lattice resolution. Despite these uncertainties, the present results are in good agreement with the results in the literature (within a 2.5% error), and thus, we confirm that the proposed method can calculate the rod motion in steady uniform fluid flows with reasonable accuracy.

FIG. 12.

Comparison of the equilibrium rod positions between experimental results (red line) and numerical results (white line). The equilibrium inclined angles of the rod obtained from the experiments and simulations are 26.5° and 25.9° for case A and 49.6° and 48.4° for case B, respectively.

FIG. 12.

Comparison of the equilibrium rod positions between experimental results (red line) and numerical results (white line). The equilibrium inclined angles of the rod obtained from the experiments and simulations are 26.5° and 25.9° for case A and 49.6° and 48.4° for case B, respectively.

Close modal

The second benchmark is the wind tunnel experiment of the flexible rod exposed to the air-flow at the Reynolds number of 673 reported by Jung et al. (2018). The experimental settings and the computational domain are illustrated in Fig. 13.

FIG. 13.

Schematic representation of the experimental setting of Jung et al. (2018).

FIG. 13.

Schematic representation of the experimental setting of Jung et al. (2018).

Close modal

The experiments were performed in a closed-type subsonic wind tunnel of square cross section with a height and width of 150 mm. The lower end of the polydimethylsiloxane (PDMS) flexible rod was mounted on a polylactic acid holder with holes at the center of the tunnel bottom. Under the influence of airflows and gravity, the flexible rod showed oscillatory motions of two DOF along the flow and cross-flow directions at the equilibrium state. The oscillatory motions of the rod and the fluid velocity fluctuation at the downstream of the rod were captured by a high-speed camera and a time-resolved particle image velocimetry (PIV) technique, respectively, and then analyzed by a fast Fourier transform (FFT) to quantify the oscillatory behavior. The experimental conditions of the fluid and structure properties are listed in Tables VI and VII.

TABLE VI.

Fluid properties (air).

Density,Kinetic viscosity,Inlet velocity, Reynolds number,
ρf (kg/m3)νf (m2/s)Uin (m/s)Re[−]
1.23 1.56 × 10−5 3.5 673 
Density,Kinetic viscosity,Inlet velocity, Reynolds number,
ρf (kg/m3)νf (m2/s)Uin (m/s)Re[−]
1.23 1.56 × 10−5 3.5 673 
TABLE VII.

Structure properties (PDMS).

Density,Young’s modulus,Poisson’s ratio,Length, Diameter, Rotational stiffness,
ρs (kg/m3) E (MPa)ν[−]L̂ (mm)D (mm)κθ (N mm)
950 1.78 0.48 70 3.0 0.41 
Density,Young’s modulus,Poisson’s ratio,Length, Diameter, Rotational stiffness,
ρs (kg/m3) E (MPa)ν[−]L̂ (mm)D (mm)κθ (N mm)
950 1.78 0.48 70 3.0 0.41 

Reynolds numbers are based on the inlet velocity Uin, the rod diameter D, and the kinetic viscosity of air νf. We reproduce the above experimental conditions by the LBM–CRM–twFSI scheme. The detailed computational configuration is described in Fig. 14. The computational domain is set to be a part of the wind tunnel as Lx × Ly × Lz = 240 × 560 × 680 lattice nodes with the lattice resolution of 0.15 mm/lattice spacing, which is equivalent to 20 lattice nodes in the rod diameter. The lower end of the rod is positioned at 120,0,100 with the rotational spring boundary condition in terms of the rotational stiffness (Spoon and Grant, 2011). The rotational stiffness of the rod κθ measured by the bending test is provided in the literature (Jung et al., 2019). The constant inlet velocity is applied, while the free-outflow condition is applied to the outlet. The no-slip wall condition is applied to the bottom boundary and the slip flow condition is applied to the remaining top and both side boundaries. A constant time step is 2.57 × 10−6 s. The gravitational acceleration g acts opposite to the Y-direction. The rod is discretized by 30 nodes and the shear factor is set to be αc = 4/3. The grid convergence study is performed for this configuration at the Reynolds number of 673 with the rigid stationary cylinder. The lattice resolutions are set to be D/Δ = 6, 10, and 20. The cylinder shape data required for the IPBB method are obtained by the fsBR scheme. Figure 15(a) shows the drag coefficient against the lattice resolution. The hydrodynamic force asymptotically approaches the limiting value with the lattice resolution of finer than 10 for the cylinder diameter. Figure 15(b) shows the comparison of the numerical and experimental drag coefficients (Williamson, 1996). The shown numerical results are obtained with the lattice resolution of D/Δ = 20. The predicted drag coefficient at the Reynolds number of 673 agrees well with the experimental one, indicating that the applied lattice configuration for the second benchmark of D/Δ = 20 is sufficiently fine to evaluate the hydrodynamic force on the cylindrical structure.

FIG. 14.

Schematic representation of the computational domain at the initial state for the case of Jung et al. (2018).

FIG. 14.

Schematic representation of the computational domain at the initial state for the case of Jung et al. (2018).

Close modal
FIG. 15.

(a) Variation of the drag coefficient for the Reynolds number of 673 against the lattice resolution. (b) Comparison of the numerical and experimental drag coefficients at the Reynolds number of 40, 200, and 673. The numerical results are obtained with the lattice resolution of 20 for the cylinder diameter.

FIG. 15.

(a) Variation of the drag coefficient for the Reynolds number of 673 against the lattice resolution. (b) Comparison of the numerical and experimental drag coefficients at the Reynolds number of 40, 200, and 673. The numerical results are obtained with the lattice resolution of 20 for the cylinder diameter.

Close modal

The time-dependent oscillatory motions of the rod and the fluid velocities around the rod in a half period are shown in Fig. 16. The rod responds by a large deformation to the hydrodynamics forces, and the periodic oscillatory motions along the flow and cross-flow directions are observed at the equilibrium state. The three-dimensional wake generated at the downstream of the rod is also observed. A large scale recirculating flow is formed and moves along in the flow direction. The rod position with the minimal bending is defined as the reference state with the time t/T = 0, where T is the oscillation period. The deformed rod shape at the reference time t/T = 0 is compared between the experimental and numerical results, which is shown in Fig. 17. The present result shows good agreement with the experiment.

FIG. 16.

Time-dependent oscillatory motions of the rod and the fluid velocity profile in a half period. The velocity magnitude is normalized by Uin.

FIG. 16.

Time-dependent oscillatory motions of the rod and the fluid velocity profile in a half period. The velocity magnitude is normalized by Uin.

Close modal
FIG. 17.

Comparison of the rod shape between the experimental (red dotted line) and numerical (gray region) results at the reference time t/T = 0.

FIG. 17.

Comparison of the rod shape between the experimental (red dotted line) and numerical (gray region) results at the reference time t/T = 0.

Close modal

The representative oscillatory motions of the rod tip position, which is equivalent to the discretized upper-end node position, normalized by the rod diameter D are shown in Fig. 18(a). The maximum differences of the oscillation amplitude for the streamwise flow and cross-flow directions are 1.05D and 0.17D, respectively. Compared with the experimental values of 1.04D and 0.3D, the present results show good agreement in the streamwise flow direction, while the oscillation amplitude in the cross-flow direction is underestimated. We consider that one of the reasons for this error comes from the measurement uncertainties of the rotational stiffness of the PDMS rod provided in the literature. The rotational stiffness was measured by the bending test of a point force applied at the free end of the horizontally clamped rod. The linear fitting was employed to the experimental deformation data obtained with varying forces, which resulted in the constant rotational stiffness independent of rotation angles. The normalized power spectral density (PSD) distributions for the oscillatory motions of the rod tip position obtained by the FFT analysis are represented in Fig. 18(b). The PSD distributions show dominant frequencies at 2.17 and 3.62 Hz for the streamwise flow and cross-flow directions, respectively. Compared with the experimental values of 2.06 and 3.09 Hz, the present results show good agreement (with a 5% error) for the streamwise flow direction, while the error due to the same measurement uncertainties stated above is observed for the cross-flow direction. As the numerical oscillation amplitude along the cross-flow direction is smaller than the experimental value, the dominant frequency exhibits a higher value in the numerical results.

FIG. 18.

(a) Oscillatory motions of the rod tip position normalized by the rod diameter D. The variation of Cz and Cx represents the motions along the streamwise flow and cross-flow directions, respectively. The left and the right axes correspond to the variation of Cz and Cx, respectively. (b) Normalized power spectral density (PSD) distributions of the oscillatory motions of the rod tip.

FIG. 18.

(a) Oscillatory motions of the rod tip position normalized by the rod diameter D. The variation of Cz and Cx represents the motions along the streamwise flow and cross-flow directions, respectively. The left and the right axes correspond to the variation of Cz and Cx, respectively. (b) Normalized power spectral density (PSD) distributions of the oscillatory motions of the rod tip.

Close modal

The time-dependent total forces acting on the rod in the cross-flow and the streamwise flow directions are shown in Figs. 19(a) and 19(b), respectively. By comparing the oscillation behavior of the tip positions with that of the total forces, we realize that the forces act approximately in the opposite directions to those of the tip motions. This implies that the oscillatory motions of the rod are induced by the hydrodynamic forces and the rotational stiffness of the rod. The rod is first bent by the hydrodynamic forces, and then the rotational spring torques corresponding to the rod bending angle are generated to act as the damping force against the rod motion.

FIG. 19.

The time-dependent total forces acting on the rod in (a) the cross-flow direction and (b) the streamwise flow direction. The forces are smoothed by the simple moving average filter.

FIG. 19.

The time-dependent total forces acting on the rod in (a) the cross-flow direction and (b) the streamwise flow direction. The forces are smoothed by the simple moving average filter.

Close modal

The normalized PSD distributions of the streamwise fluid velocity component at the downstream of the rod are shown in Fig. 20. The fluid velocities are obtained in the far wake region (FWR) located at x/D = 1, y/D = −1, and z/D = 3 and in the near wake regions (NWRs) located at x/D = 0, y/D = −1, and z/D = 0.8, 1.1, and 1.4, which correspond to Figs. 20(a) and 20(b), respectively. The PSD distributions of the FWR show two peaks at 2.17 and 208 Hz. Compared with the experimental results of 3.09 and 213 Hz, the higher frequency agrees well with the literature (within a 2.5% error), while each lower frequency is influenced by the oscillatory motions of the rod in two different directions. The numerical value of 2.17 Hz is equivalent to the numerical dominant frequency of the oscillatory motion of the rod tip in the flow direction, indicating that the FWR streamwise velocity is influenced by the oscillatory rod motion in the flow direction. On the other hand, the experimental value of 3.09 Hz is equivalent to the experimental dominant frequency of the oscillatory motion of the rod tip in the cross-flow direction. Thus, they indicate that the FWR streamwise velocity is related to the oscillatory rod motion in the cross-flow direction. Although further discussions are required, we believe that the FWR streamwise velocity is reasonably influenced by the oscillatory rod motion in the flow direction due to the larger oscillation amplitude compared with the cross-flow direction. For the NWRs, although the dominant frequency peak is not observed, many small-scale vortices generated at the NWRs as shown in Fig. 16 are observed. The present results represent the same tendency as the experimental results. Despite the difficulties in precisely determining the material properties of the rods and the fluid velocity profiles by the experimental measurements, we confirm that the comparison between experiments and simulations is overall in good agreement.

FIG. 20.

Normalized power spectral density (PSD) distributions of the streamwise velocity component extracted at the downstream location of (a) the far wake region (FWR) and (b) the near wake regions (NWRs). The FWR is located at x/D = 1, y/D = −1, and z/D = 3, while the NWR1, 2, and 3 are located at x/D = 0, y/D = −1, and z/D = 0.8, 1.1, and 1.4, respectively.

FIG. 20.

Normalized power spectral density (PSD) distributions of the streamwise velocity component extracted at the downstream location of (a) the far wake region (FWR) and (b) the near wake regions (NWRs). The FWR is located at x/D = 1, y/D = −1, and z/D = 3, while the NWR1, 2, and 3 are located at x/D = 0, y/D = −1, and z/D = 0.8, 1.1, and 1.4, respectively.

Close modal

In this study, a numerical scheme to simulate three-dimensional twFSI problems of flows around a flexible fine rod is developed. The D3Q27 MRT-LBM and the one-dimensional geometrically exact CRM are employed to compute the fluid flow and the rod motion, respectively. The partitioned approach is employed to handle the fluid and the structure solvers separately, and the fluid–structure interactions are calculated by the simple explicit coupling scheme combined with the contact detection algorithm and the fsBR scheme. The contact detection algorithm utilizing the BVH is adopted to reduce the computing time of exchanging information between the fluid solver and the structure solver, while the fsBR scheme utilizes the level set method to reconstruct the fluid–structure interface. The fluid solver and the structure solver are implemented by the in-house MPI-based multi-GPU code written in CUDA-Fortran programming language and the multi-thread code written in C++ programming language, respectively. The proposed LBM–CRM–twFSI scheme stably and efficiently calculates rod motions in fluid flows owing to its simple algorithm.

The fluid and the structure solvers are first verified individually. For the fluid solver, flows around the circular cylinder at the Reynolds numbers of 40 and 200 are computed by the LBM combined with the fsBR scheme, and then the drag and the lift coefficients obtained by the GI-MEM are compared with the literature values. The results indicate that the LBM combined with the fsBR scheme and the GI-MEM accurately evaluates the hydrodynamic force on the cylindrical structure with the lattice resolution of 20 in the cylinder diameter. For the structure solver, a cantilever beam with circular cross sections horizontally clamped at one end and subjected to the downward uniform force is computed by the CRM, and then the deformation is compared with the analytical solution of the Timoshenko cantilever beam theory. The results show good agreement that the CRM well represents the rod deformation. Finally, the proposed LBM–CRM–twFSI scheme is successfully validated in two experimental benchmarks of a single flexible rod deformation in a wind tunnel. The results confirm that the present scheme accurately calculates the equilibrium state of the rod exposed to uniform airflows at the Reynolds numbers of 29.1 and 55.6. The time-dependent oscillatory motions of the rod and the fluid velocity fluctuation at the downstream of the rod at the Reynolds number of 673 are also reasonably captured by the present scheme. We confirm that the errors of the representative streamwise rod position between the experimental and numerical results for both benchmarks are within 5%. These validations confirm the practicability of the presently developed LBM–CRM–twFSI scheme for motions of flexible fine rod in fluid flows.

The technical support by Dr. Y. Kuwata and Mr. M. Sugimoto is appreciated.

The authors declare that there is no conflict of interest with respect to the research, authorship, and/or publication of this article.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

All authors contributed equally to this work.

The data that support the findings of this study are available within the article.

1.
Abdi
,
R.
,
Rezazadeh
,
N.
, and
Abdi
,
M.
, “
Investigation of passive oscillations of flexible splitter plates attached to a circular cylinder
,”
J. Fluids Struct.
84
,
302
317
(
2019
).
2.
Agarwal
,
A.
,
Gupta
,
S.
, and
Prakash
,
A.
, “
A comparative study of bounce-back and immersed boundary method in LBM for turbulent flow simulation
,”
Mater. Today Proc.
28
,
2387
2392
(
2020
).
3.
Aidun
,
C. K.
and
Clausen
,
J. R.
, “
Lattice-Boltzmann method for complex flows
,”
Annu. Rev. Fluid Mech.
42
,
439
472
(
2010
).
4.
Ando
,
S.
,
Kaneda
,
M.
, and
Suga
,
K.
, “
Permeability prediction of fibrous porous media by the lattice Boltzmann method with a fluid-structure boundary reconstruction scheme
,”
J. Ind. Text.
152808372097891
(
2020
)
5.
Bavo
,
A. M.
,
Rocatello
,
G.
,
Iannaccone
,
F.
,
Degroote
,
J.
,
Vierendeels
,
J.
, and
Segers
,
P.
, “
Fluid-structure interaction simulation of prosthetic aortic valves: Comparison between immersed boundary and arbitrary Lagrangian-Eulerian techniques for the mesh representation
,”
PLoS One
11
,
e0154517
(
2016
).
6.
Bourguet
,
R.
and
Triantafyllou
,
M. S.
, “
The onset of vortex-induced vibrations of a flexible cylinder at large inclination angle
,”
J. Fluid Mech.
809
,
111
134
(
2016
).
7.
Bouzidi
,
M. H.
,
Firdaouss
,
M.
, and
Lallemand
,
P.
, “
Momentum transfer of a Boltzmann-lattice fluid with boundaries
,”
Phys. Fluids
13
,
3452
3459
(
2001
).
8.
Caiazzo
,
A.
, “
Analysis of lattice Boltzmann nodes initialisation in moving boundary problems
,”
Prog. Comput. Fluid Dyn.
8
,
3
10
(
2008
).
9.
Calderer
,
A.
,
Kang
,
S.
, and
Sotiropoulos
,
F.
, “
Level set immersed boundary method for coupled simulation of air/water interaction with complex floating structures
,”
J. Comput. Phys.
277
,
201
227
(
2014
).
10.
Cao
,
D. Q.
,
Liu
,
D.
, and
Wang
,
C. H.-T.
, “
Three-dimensional nonlinear dynamics of slender structures: Cosserat rod element approach
,”
Int. J. Solids Struct.
43
,
760
783
(
2006
).
11.
Cao
,
D. Q.
and
Tucker
,
R. W.
, “
Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation
,”
Int. J. Solids Struct.
45
,
460
477
(
2008
).
12.
Cerroni
,
D.
,
Da Vià
,
R.
, and
Manservisi
,
S.
, “
A projection method for coupling two-phase VOF and fluid structure interaction simulations
,”
J. Comput. Phys.
354
,
646
671
(
2018
).
13.
Coleman
,
B. D.
,
Dill
,
E. H.
,
Lembo
,
M.
,
Lu
,
Z.
, and
Tobias
,
I.
, “
On the dynamics of rods in the theory of Kirchhoff and Clebsch
,”
Arch. Ration. Mech. Anal.
121
,
339
359
(
1993
).
14.
Cosserat
,
E.
and
Cosserat
,
F.
,
Theory of Deformable Bodies
(
Cornell University Library
,
1909
), pp.
1
246
.
15.
D’Humières
,
D.
,
Ginzburg
,
I.
,
Krafczyk
,
M.
,
Lallemand
,
P.
, and
Luo
,
L. S.
, “
Multiple-relaxation-time lattice Boltzmann models in three dimensions
,”
Philos. Trans. R. Soc., A
360
,
437
451
(
2002
).
16.
De Rosis
,
A.
,
Ubertini
,
S.
, and
Ubertini
,
F.
, “
A comparison between the interpolated bounce-back scheme and the immersed boundary method to treat solid boundary conditions for laminar flows in the lattice Boltzmann framework
,”
J. Sci. Comput.
61
,
477
489
(
2014
).
17.
Degroote
,
J.
,
Bathe
,
K. J.
, and
Vierendeels
,
J.
, “
Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction
,”
Comput. Struct.
87
,
793
801
(
2009
).
18.
Dennis
,
S. C. R.
and
Chang
,
G.-Z.
, “
Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100
,”
J. Fluid Mech.
42
,
471
489
(
1970
).
19.
Dey
,
A. A.
,
Modarres-Sadeghi
,
Y.
, and
Rothstein
,
J. P.
, “
Viscoelastic fluid-structure interactions between a flexible cylinder and wormlike micelle solution
,”
Phys. Rev. Fluids
3
,
063301
(
2018
).
20.
Dill
,
E. H.
, “
Kirchhoff’s theory of rods
,”
Arch. Hist. Exact Sci.
44
,
1
23
(
1906
).
21.
Dukowicz
,
J. K.
and
Kodis
,
J. W.
, “
Accurate conservative remapping (rezoning) for arbitrary Lagrangian-Eulerian computations
,”
SIAM J. Sci. Stat. Comput.
8
,
305
321
(
1987
).
22.
Duričković
,
B.
,
Goriely
,
A.
, and
Maddocks
,
J. H.
, “
Twist and stretch of helices explained via the Kirchhoff-Love rod model of elastic filaments
,”
Phys. Rev. Lett.
111
,
108103
(
2013
).
23.
Farhat
,
C.
and
Lesoinne
,
M.
, “
Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems
,”
Comput. Methods Appl. Mech. Eng.
182
,
499
515
(
2000
).
24.
Favier
,
J.
,
Revell
,
A.
, and
Pinelli
,
A.
, “
A lattice Boltzmann–immersed boundary method to simulate the fluid interaction with moving and slender flexible objects
,”
J. Comput. Phys.
261
,
145
161
(
2014
).
25.
Felippa
,
C. A.
,
Park
,
K. C.
, and
Farhat
,
C.
, “
Partitioned analysis of coupled mechanical systems
,”
Comput. Methods Appl. Mech. Eng.
190
,
3247
3270
(
2001
).
26.
Fernández
,
M. A.
, “
Coupling schemes for incompressible fluid-structure interaction: Implicit, semi-implicit and explicit
,”
SeMA J.
55
,
59
108
(
2011
).
27.
Frisch
,
U.
,
Hasslacher
,
B.
, and
Pomeau
,
Y.
, “
Lattice-gas automata for the Navier-Stokes equation
,”
Phys. Rev. Lett.
56
,
1505
1508
(
1986
).
28.
Gao
,
B.
,
Hu
,
K.
, and
Guo
,
S.
, “
Collision detection algorithm based on AABB for minimally invasive surgery
,” in
Proceedings of 2014 IEEE International Conference on Mechatronics and Automation
(
IEEE
,
2014
), pp.
315
320
.
29.
Gao
,
H.
,
Li
,
H.
, and
Wang
,
L.-P.
, “
Lattice Boltzmann simulation of turbulent flow laden with finite-size particles
,”
Comput. Math. Appl.
65
,
194
210
(
2013
).
30.
Gazzola
,
M.
,
Dudte
,
L. H.
,
McCormick
,
A. G.
, and
Mahadevan
,
L.
, “
Forward and inverse problems in the mechanics of soft filaments
,”
R. Soc. Open Sci.
5
,
171628
(
2018
).
31.
Goicoechea
,
H. E.
,
Buezas
,
F. S.
, and
Rosales
,
M. B.
, “
A non-linear Cosserat rod model for drill-string dynamics in arbitrary borehole geometries with contact and friction
,”
Int. J. Mech. Sci.
157-158
,
98
110
(
2019
).
32.
Ha
,
S. T.
,
Ngo
,
L. C.
,
Saeed
,
M.
,
Jeon
,
B. J.
, and
Choi
,
H.
, “
A comparative study between partitioned and monolithic methods for the problems with 3D fluid-structure interaction of blood vessels
,”
J. Mech. Sci. Technol.
31
,
281
287
(
2017
).
33.
He
,
Q.
,
Li
,
Y.
,
Huang
,
W.
,
Hu
,
Y.
,
Li
,
D.
, and
Wang
,
Y.
, “
Lattice Boltzmann model for dense suspended particles based on improved bounce-back method
,”
Comput. Math. Appl.
80
,
552
567
(
2020
).
34.
Hutchinson
,
J. R.
, “
Shear coefficients for Timoshenko beam theory
,”
J. Appl. Mech.
68
,
87
92
(
2001
).
35.
Jung
,
P.
,
Leyendecker
,
S.
,
Linn
,
J.
, and
Ortiz
,
M.
, “
A discrete mechanics approach to the Cosserat rod theory—Part 1: Static equilibria
,”
Int. J. Numer. Methods Eng.
85
,
31
60
(
2012
).
36.
Jung
,
S. Y.
,
Kim
,
J. J.
, and
Lee
,
S. J.
, “
Effect of material stiffness on the motion and flow near the free end of a finite cylinder surface
,”
Exp. Therm. Fluid Sci.
102
,
548
558
(
2019
).
37.
Jung
,
S. Y.
,
Kim
,
J. J.
,
Park
,
H. W.
, and
Lee
,
S. J.
, “
Comparison of flow structures behind rigid and flexible finite cylinders
,”
Int. J. Mech. Sci.
142-143
,
480
490
(
2018
).
38.
Karras
,
T.
, “
Maximizing parallelism in the construction of BVHs, octrees, and k-d trees
,” in
Proceedings of the High-Performance Graphics 2012, HPG 2012—ACM SIGGRAPH/Eurographics Symposium
(
HPG
,
2012
), pp.
33
37
.
39.
Kassiotis
,
C.
,
Ibrahimbegovic
,
A.
, and
Matthies
,
H.
, “
Partitioned solution to fluid structure interaction problem in application to free-surface flows
,”
Eur. J. Mech. - B/Fluids
29
,
510
521
(
2010
).
40.
Ketchel
,
J.
and
Larochelle
,
P.
, “
Collision detection of cylindrical rigid bodies for motion planning
,” in
Proceedings of the 2006 IEEE International Conference on Robotics and Automation
(
IEEE
,
2006
), pp.
1530
1535
.
41.
Krithivasan
,
S.
,
Wahal
,
S.
, and
Ansumali
,
S.
, “
Diffused bounce-back condition and refill algorithm for the lattice Boltzmann method
,”
Phys. Rev. E
89
,
033313
(
2014
).
42.
Lallemand
,
P.
and
Luo
,
L.-S.
, “
Lattice Boltzmann method for moving boundaries
,”
J. Comput. Phys.
184
,
406
421
(
2003
).
43.
Legay
,
A.
,
Chessa
,
J.
, and
Belytschko
,
T.
, “
An Eulerian-Lagrangian method for fluid-structure interaction based on level sets
,”
Comput. Methods Appl. Mech. Eng.
195
,
2070
2087
(
2006
).
44.
Levinson
,
M.
, “
A new rectangular beam theory
,”
J. Sound Vib.
74
,
81
87
(
1981
).
45.
Li
,
G.
,
Staszel
,
C.
,
Yarin
,
A. L.
, and
Pourdeyhimi
,
B.
, “
Hydroentanglement of polymer nonwovens. 1: Experimental and theoretical/numerical framework
,”
Polymer
164
,
191
204
(
2019
).
46.
Li
,
W.
,
Wang
,
W.-Q.
,
Yan
,
Y.
, and
Yu
,
Z.-F.
, “
A strong-coupled method combined finite element method and lattice Boltzmann method via an implicit immersed boundary scheme for fluid structure interaction
,”
Ocean Eng.
214
,
107779
(
2020
).
47.
Li
,
Y.
,
Shock
,
R.
,
Zhang
,
R.
, and
Chen
,
H.
, “
Numerical study of flow past an impulsively started cylinder by the lattice-Boltzmann method
,”
J. Fluid Mech.
519
,
273
300
(
2004
).
48.
Lim
,
S.
and
Peskin
,
C. S.
, “
Fluid-mechanical interaction of flexible bacterial flagella by the immersed boundary method
,”
Phys. Rev. E
85
,
036307
(
2012
).
49.
Liu
,
C.
,
Zheng
,
X.
, and
Sung
,
C. H.
, “
Preconditioned multigrid methods for unsteady incompressible flows
,”
J. Comput. Phys.
139
,
35
57
(
1998
).
50.
Lorenz
,
E.
,
Sivadasan
,
V.
,
Bonn
,
D.
, and
Hoekstra
,
A. G.
, “
Combined lattice–Boltzmann and rigid-body method for simulations of shear-thickening dense suspensions of hard particles
,”
Comput. Fluids
172
,
474
482
(
2018
).
51.
Love
,
A.
,
A Treatise on the Mathematical Theory of Elasticity
(
Cambridge University Press
,
1892
), Vol. 1, pp.
1
368
.
52.
Lv
,
X.
,
Zhao
,
Y.
,
Huang
,
X. Y.
,
Xia
,
G. H.
, and
Wang
,
Z. J.
, “
An efficient parallel/unstructured-multigrid preconditioned implicit method for simulating 3D unsteady compressible flows with moving objects
,”
J. Comput. Phys.
215
,
661
690
(
2006
).
53.
Margolin
,
L. G.
and
Shashkov
,
M.
, “
Second-order sign-preserving conservative interpolation (remapping) on general grids
,”
J. Comput. Phys.
184
,
266
298
(
2003
).
54.
Marheineke
,
N.
and
Wegener
,
R.
, “
Modeling and application of a stochastic drag for fibers in turbulent flows
,”
Int. J. Multiphase Flow
37
,
136
148
(
2011
).
55.
Mattis
,
S. A.
,
Dawson
,
C. N.
,
Kees
,
C. E.
, and
Farthing
,
M. W.
, “
An immersed structure approach for fluid-vegetation interaction
,”
Adv. Water Resour.
80
,
1
16
(
2015
).
56.
Meier
,
C.
,
Popp
,
A.
, and
Wall
,
W. A.
, “
Geometrically exact finite element formulations for slender beams: Kirchhoff–Love theory versus Simo–Reissner theory
,”
Arch. Comput. Methods Eng.
26
,
163
(
2019
).
57.
Mokbel
,
D.
,
Abels
,
H.
, and
Aland
,
S.
, “
A phase-field model for fluid–structure interaction
,”
J. Comput. Phys.
372
,
823
840
(
2018
).
58.
Niu
,
X. D.
,
Chew
,
Y. T.
, and
Shu
,
C.
, “
Simulation of flows around an impulsively started circular cylinder by Taylor series expansion- and least squares-based lattice Boltzmann method
,”
J. Comput. Phys.
188
,
176
193
(
2003
).
59.
Niu
,
X. D.
,
Shu
,
C.
,
Chew
,
Y. T.
, and
Peng
,
Y.
, “
A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows
,”
Phys. Lett. A
354
,
173
182
(
2006
).
60.
Obrecht
,
C.
,
Kuznik
,
F.
,
Tourancheau
,
B.
, and
Roux
,
J.-J.
, “
Efficient GPU implementation of the linearly interpolated bounce-back boundary condition
,”
Comput. Math. Appl.
65
,
936
944
(
2013
).
61.
Owen
,
D. R. J.
,
Leonardi
,
C. R.
, and
Feng
,
Y. T.
, “
An efficient framework for fluid–structure interaction using the lattice Boltzmann method and immersed moving boundaries
,”
Int. J. Numer. Methods Eng.
87
,
66
95
(
2011
).
62.
Pan
,
C.
,
Luo
,
L.-S.
, and
Miller
,
C. T.
, “
An evaluation of lattice Boltzmann schemes for porous medium flow simulation
,”
Comput. Fluids
35
,
898
909
(
2006
).
63.
Park
,
K. C.
,
Felippa
,
C. A.
, and
DeRuntz
,
J. A.
, “
Stabilization of staggered solution procedures for fluid-structure interaction analysis
,” in
Computational Methods for Fluid-Structure Interaction Problems
(
American Society of Mechanical Engineers, Applied Mechanics Division, AMD
,
1977
), pp.
95
124
.
64.
Patel
,
H. V.
,
Das
,
S.
,
Kuipers
,
J. A. M.
,
Padding
,
J. T.
, and
Peters
,
E. A. J. F.
, “
A coupled volume of fluid and immersed boundary method for simulating 3D multiphase flows with contact line dynamics in complex geometries
,”
Chem. Eng. Sci.
166
,
28
41
(
2017
).
65.
Peng
,
C.
,
Teng
,
Y.
,
Hwang
,
B.
,
Guo
,
Z.
, and
Wang
,
L.-P.
, “
Implementation issues and benchmarking of lattice Boltzmann method for moving rigid particle simulations in a viscous flow
,”
Comput. Math. Appl.
72
,
349
374
(
2016
).
66.
Qian
,
Y. H.
, “
Simulating thermohydrodynamics with lattice BGK models
,”
J. Sci. Comput.
8
,
231
242
(
1993
).
67.
Rinaldi
,
P. R.
,
Dari
, E. A.
,
Vénere
,
M. J.
, and
Clausse
,
A.
, “
A lattice-Boltzmann solver for 3D fluid simulation on GPU
,”
Simul. Model. Pract. Theory
25
,
163
171
(
2012
).
68.
Rubin
,
M. B.
, “
Cosserat theories: Shells, rods and points
,”
Appl. Mech. Rev.
55
,
B109
(
2000
).
69.
Schneider
,
P.
and
Eberly
,
D.
,
Geometric Tools for Computer Graphics
(
Morgan Kaufmann Publishers
,
2003
).
70.
Shen
,
L.
and
Chan
,
E.-S.
, “
Numerical simulation of fluid-structure interaction using a combined volume of fluid and immersed boundary method
,”
Ocean Eng.
35
,
939
952
(
2008
).
71.
Silva-Leon
,
J.
,
Cioncolini
,
A.
, and
Filippone
,
A.
, “
Determination of the normal fluid load on inclined cylinders from optical measurements of the reconfiguration of flexible filaments in flow
,”
J. Fluids Struct.
76
,
488
505
(
2018
).
72.
Spillmann
,
J.
and
Teschner
,
M.
, “
CORDE: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects
,” in
Eurographics/ACM SIGGRAPH Symposium on Computer Animation
(
SCA
,
2007
), pp.
1
10
.
73.
Spoon
,
C.
and
Grant
,
W.
, “
Biomechanics of hair cell kinocilia: Experimental measurement of kinocilium shaft stiffness and base rotational stiffness with Euler–Bernoulli and Timoshenko beam analysis
,”
J. Exp. Biol.
214
,
862
870
(
2011
).
74.
Suga
,
K.
,
Kuwata
,
Y.
,
Takashima
,
K.
, and
Chikasue
,
R.
, “
A D3Q27 multiple-relaxation-time lattice Boltzmann method for turbulent flows
,”
Comput. Math. Appl.
69
,
518
529
(
2015
).
75.
Suga
,
K.
,
Tanaka
,
T.
,
Nishio
,
Y.
, and
Murata
,
M.
, “
A boundary reconstruction scheme for lattice Boltzmann flow simulation in porous media
,”
Prog. Comput. Fluid Dyn.
9
,
201
207
(
2009
).
76.
Sun
,
P.
,
Xu
,
J.
, and
Zhang
,
L.
, “
Full Eulerian finite element method of a phase field model for fluid-structure interaction problem
,”
Comput. Fluids
90
,
1
8
(
2014
).
77.
Sussman
,
M.
,
Almgren
,
A. S.
,
Bell
,
J. B.
,
Colella
,
P.
,
Howell
,
L. H.
, and
Welcome
,
M.
, “
An adaptive level set approach for incompressible two-phase flows
,”
Am. Soc. Mech. Eng. Fluids Eng. Div. FED
238
,
355
360
(
1996
).
78.
Suzuki
,
K.
and
Inamuro
,
T.
, “
Effect of internal mass in the simulation of a moving body by the immersed boundary method
,”
Comput. Fluids
49
,
173
187
(
2011
).
79.
Tai
,
C. H.
and
Zhao
,
Y.
, “
Parallel unsteady incompressible viscous flow computations using an unstructured multigrid method
,”
J. Comput. Phys.
192
,
277
311
(
2003
).
80.
Tan
,
W.
,
Wu
,
H.
, and
Zhu
,
G.
, “
Fluid-structure interaction using lattice Boltzmann method: Moving boundary treatment and discussion of compressible effect
,”
Chem. Eng. Sci.
184
,
273
284
(
2018
).
81.
Tao
,
S.
,
Hu
,
J.
, and
Guo
,
Z.
, “
An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows
,”
Comput. Fluids
133
,
1
14
(
2016
).
82.
Tian
,
F.-B.
,
Dai
,
H.
,
Luo
,
H.
,
Doyle
,
J. F.
, and
Rousseau
,
B.
, “
Fluid-structure interaction involving large deformations: 3D simulations and applications to biological systems
,”
J. Comput. Phys.
258
,
451
469
(
2014
).
83.
Timoshenko
,
S. P.
, “
LXVI. On the correction for shear of the differential equation for transverse vibrations of prismatic bars
,”
Philos. Mag.
41
,
744
746
(
1921
).
84.
Tschisgale
,
S.
and
J.
Fröhlich
, “
An immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid
,”
J. Comput. Phys.
423
,
109801
(
2020
).
85.
van den Bergen
,
G.
, “
Efficient collision detection of complex deformable models using AABB trees
,”
J. Graph. Tools
2
,
131
(
1997
).
86.
Wen
,
B.
,
Zhang
,
C.
, and
Fang
,
H.
, “
Hydrodynamic force evaluation by momentum exchange method in lattice Boltzmann simulations
,”
Entropy
17
,
8240
8266
(
2015
).
87.
Wen
,
B.
,
Zhang
,
C.
,
Tu
,
Y.
,
Wang
,
C.
, and
Fang
,
H.
, “
Galilean invariant fluid-solid interfacial dynamics in lattice Boltzmann simulations
,”
J. Comput. Phys.
266
,
161
170
(
2014
).
88.
Wick
,
T.
, “
Coupling fluid–structure interaction with phase-field fracture
,”
J. Comput. Phys.
327
,
67
96
(
2016
).
89.
Williamson
,
C. H. K.
, “
Vortex dynamics in the cylinder wake
,”
Annu. Rev. Fluid Mech.
28
,
477
539
(
1996
).
90.
Wolfram
,
S.
, “
Cellular automaton fluids 1: Basic theory
,”
J. Stat. Phys.
45
,
471
526
(
1986
).
91.
Wu
,
J.
,
Cheng
,
Y.
,
Zhang
,
C.
, and
Diao
,
W.
, “
Simulating vortex induced vibration of an impulsively started flexible filament by an implicit IB-LB coupling scheme
,”
Comput. Math. Appl.
79
,
159
173
(
2020
).
92
Wu
,
J.
,
Qiu
,
Y. L.
,
Shu
,
C.
, and
Zhao
,
N.
, “
Flow control of a circular cylinder by using an attached flexible filament
,”
Phys. Fluids
26
,
103601
(
2014
).
93.
Zhang
,
X.
,
Chan
,
F. K.
,
Parthasarathy
,
T.
, and
Gazzola
,
M.
, “
Modeling and simulation of complex dynamic musculoskeletal architectures
,”
Nat. Commun.
10
,
4825
(
2019
).
94.
Zhang
,
X.
and
Kim
,
Y. J.
, “
Interactive collision detection for deformable models using streaming AABBs
,”
IEEE Trans. Vis. Comput. Graph.
13
,
318
329
(
2007
).
95.
Zhu
,
L.
,
He
,
G.
,
Wang
,
S.
,
Miller
,
L.
,
Zhang
,
X.
,
You
,
Q.
, and
Fang
,
S.
, “
An immersed boundary method based on the lattice Boltzmann approach in three dimensions, with application
,”
Comput. Math. Appl.
61
,
3506
3518
(
2011
).
96.
Zhu
,
L.
and
Peskin
,
C. S.
, “
Interaction of two flapping filaments in a flowing soap film
,”
Phys. Fluids
15
,
1954
1960
(
2003
).