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.

## I. INTRODUCTION

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.

## II. SOLVERS FOR FLUID AND STRUCTURE

### A. Lattice Boltzmann method

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

where *f* is the density distribution function, Ω is the collision operator, ** x** is the lattice position,

*e*_{α}is the discrete velocity, subscript $\alpha =0,1,\u2026,Q\u22121$ 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,t\u2208RQ$ is equivalent to

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

where ** M** is the transformation matrix, $S\u0302$ is the relaxation matrix, $mx,t$ is the moment, and $meqx,t$ is the equilibrium moment. The transformation matrix $M\u2208RQ\xd7Q$ linearly transforms the density distribution function

*f*to the moment as $m=Mf$ and as $meq=Mfeq$, where

*f*

^{eq}is the equilibrium distribution function typically determined by the density

*ρ*

_{f}and the velocity

*u*_{f}of the fluid. The equilibrium distribution function

*f*

^{eq}is given as

where *w*_{α} is the weight coefficient and *c*_{s} 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 *c*_{s} are associated with the lattice discretizing models (Qian, 1993). In this study, the D3Q27 model is employed, and its discrete velocity vector *e*_{α} is

and Table I lists the weight coefficient *w*_{α}.

. | w_{α}
. | ||
---|---|---|---|

8/27 | $\alpha =0$ | ||

2/27 | $\alpha =1,\u2026,6$ | ||

1/54 | $\alpha =7,\u2026,18$ | ||

1/216 | $\alpha =19,\u2026,26$ |

. | w_{α}
. | ||
---|---|---|---|

8/27 | $\alpha =0$ | ||

2/27 | $\alpha =1,\u2026,6$ | ||

1/54 | $\alpha =7,\u2026,18$ | ||

1/216 | $\alpha =19,\u2026,26$ |

The transformation matrix ** M** of the D3Q27 model is set the same as that in Suga

*et al.*(2015). The collision matrix $S\u0302\u2208RQ\xd7Q$, which is a diagonal matrix, is given as

where *s*_{α} are chosen as *s*_{4} = 1.54, *s*_{5} = *s*_{7}, *s*_{10} = 1.5, *s*_{13} = 1.83, *s*_{16} = 1.4, *s*_{17} = 1.61, *s*_{18} = *s*_{20} = 1.98, and *s*_{23} = *s*_{26} = 1.74. The fluid density *ρ*_{f}, velocity *u*_{f}, and pressure *p* are obtained from the conservation law and the equation of state as follows:

The fluid kinematic viscosity *ν*_{f} is given as

where *τ*_{s} is the hydrodynamic relaxation time equivalent to 1/*s*_{5} and 1/*s*_{7}.

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.

### B. Cosserat rod model

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,t\u2208R3$ and a cross-sectional orientation characterized by a director frame $Rs,t\u2208SO3$, 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

**in the Eulerian frame to a vector $aL$ in the Lagrangian frame via $aL=Ra$ and, conversely, $a=RTaL$. Here, $s\u22080,L\u2208R$ is the center line arc-length coordinate in the current configuration,**

*a**L*is the characteristic length, and $t\u2208R+$ 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**

*κ***as**

*R*where $vecV$ is the 3-vector associated with the skew-symmetric matrix ** V**. Denoting $s\u0302$ 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 $es\u0302,t=ds/ds\u0302$ (hatted symbol indicates the reference configuration, hereafter). As

*d*_{3}is the normal to the cross section and

*d*_{1}and

*d*_{2}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

**, and thus,**

*t*

*d*_{3}≠

**. Here, the vector $\sigma L$ is defined as $\sigma L=Ret\u2212d3$ to determine the axial and shear strains. Note that the conditions**

*t**e*= 1 and

*d*_{3}=

**are equivalent to the Kirchhoff theory for inextensible and unshearable rods (Coleman**

*t**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

where $\tau L$ is the internal moment, $nL$ is the internal force, $B\u2208R3\xd73=diagB1,B2,B3$ is the bending/twisting stiffness matrix, and $S\u2208R3\xd73=diagS1,S2,S3$ is the shearing/stretching stiffness matrix. Here, *B*_{1} = *EI*_{1} is the bending rigidity about *d*_{1}, *B*_{2} = *EI*_{2} is the bending rigidity about *d*_{2}, *B*_{3} = *GI*_{3} is the twisting rigidity about *d*_{3}, *S*_{1} = *α*_{c}*GA* is the shearing rigidity along *d*_{1}, *S*_{2} = *α*_{c}*GA* is the shearing rigidity along *d*_{2}, and *S*_{3} = *EA* is the stretching rigidity along *d*_{3} 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 $I\u2208R3\xd73=diagI1,I2,I3$. The initial strain state is expressed by $\kappa Lt=0$ and $\sigma 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

where $dm=\rho sA\u0302ds\u0302=\rho sAds$ is the infinitesimal mass element, $dJ\u0302=\rho sI\u0302ds\u0302$ is the infinitesimal mass second moment of the area with the rod density *ρ*_{s}, $F=efds\u0302$ is the external force, and $CL=ecLds\u0302$ 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 $rit\u2208R3,i\u22081,n+1$ and a set of director frames $Rit\u2208R3\xd73=d1i,d2i,d3i,i\u22081,n$. Each node is represented by the position vector

*r*_{i}and the mass

*m*

_{i}, while each connecting cylindrical segment is represented by the position vectors at both ends $ri,ri+1$ and the current length

*l*

_{i}. 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.

## III. FLUID–STRUCTURE COUPLING SCHEME

### A. Fluid–structure coupling scheme

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:

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

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

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.

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

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.

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.

### B. Hydrodynamics force evaluation

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

where *F*_{α} is the hydrodynamic force and subscripts *α* and *α*′ represent directions opposite to each other, which are defined as $e\alpha =\u2212e\alpha \u2032$. This indicates that the mass component $f\alpha x,t$ has the velocity $e\alpha \u2212uW$ relative to the boundary and increases a momentum $e\alpha \u2212uWf\alpha x,t$ to the boundary, while the mass component $f\alpha \u2032x,t+\delta t$ has the relative velocity $e\alpha \u2032\u2212uW$ and decreases a momentum $e\alpha \u2032\u2212uWf\alpha \u2032x,t+\delta t$ from the boundary.

### C. Contact detection

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

where *F*_{i} is the external force, *C*_{i} is the external moments, *F*_{f} is the hydrodynamic force calculated by the summation of the GI-MEM forces *∑**F*_{α}, *x*_{b} is the position of the fluid–structure interface lattice nodes, *X*_{G} is the position of the discretized rod centroids, and *N*_{b} 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:

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.

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.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.

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.

### D. Refilling

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

where $f\alpha neq$ is the non-equilibrium distribution function, $\rho \u0304f$ is the averaged fluid density of available neighboring fluid lattice nodes, and *c*_{n} is the extrapolation direction. The equilibrium distribution is calculated in terms of the local wall velocity *u*_{W} and the averaged fluid density $\rho \u0304f$ 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 *c*_{N}. Here, *c*_{N} is a specified discrete velocity along which direction the quantity *n*_{W} · *c*_{N} takes the maximum value, where *n*_{W} is the unit normal vector pointing outwards of the moving boundary from the local wall surface.

### E. Fluid–structure boundary reconstruction

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):

where *t*_{init} is the virtual time step during the re-initialization process and *n*_{b} is the unit normal to the fluid–structure interface and given as

with initial conditions $\varphi x,0=\varphi 0x$.

## IV. RESULTS AND DISCUSSION

### A. Verification of the fluid solver

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 *Re*_{D} = *U*_{in}*D*/*ν*_{f}, where *U*_{in} 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 40*D* × 40*D* 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.

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

where *F*_{D} is the drag force and *F*_{L} 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.

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.

Re
. | Researchers . | C_{D}
. | C_{L}
. |
---|---|---|---|

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 |

Re
. | Researchers . | C_{D}
. | C_{L}
. |
---|---|---|---|

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.

### B. Verification of the structure solver

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 $s\u0302=0$ and subjected to the downward uniform force *F*_{U} is considered, as illustrated in Fig. 7. The theoretical solution of the displacement *y* is analytically expressed as

where $L\u0302$ is the length of the rod, $A\u0302$is the constant cross-sectional area, $I\u0302$ 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 $s\u0302$. See Hutchinson (2001), Levinson (1981), and Timoshenko (1921) for more details.

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.

Case . | Density, ρ_{s} (kg/m^{3})
. | Young’s modulus, E (Pa)
. | Poisson’s ratio, ν[−]
. | Length, $L\u0302$ (mm) . | Diameter, D (mm)
. | Force, F_{U} (N)
. |
---|---|---|---|---|---|---|

1-1 | 1000 | 2.25 × 10^{5} | 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 × 10^{6} | 0.48 | 70 | 3.0 | 0.5 × 10^{−3} |

2-2 | 1.0 × 10^{−3} | |||||

2-3 | 1.5 × 10^{−3} |

Case . | Density, ρ_{s} (kg/m^{3})
. | Young’s modulus, E (Pa)
. | Poisson’s ratio, ν[−]
. | Length, $L\u0302$ (mm) . | Diameter, D (mm)
. | Force, F_{U} (N)
. |
---|---|---|---|---|---|---|

1-1 | 1000 | 2.25 × 10^{5} | 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 × 10^{6} | 0.48 | 70 | 3.0 | 0.5 × 10^{−3} |

2-2 | 1.0 × 10^{−3} | |||||

2-3 | 1.5 × 10^{−3} |

### C. Validation of the LBM–CRM–twFSI scheme

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.

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.

. | Density, . | Kinetic viscosity, . | Inlet velocity, . | Reynolds number . |
---|---|---|---|---|

Case . | ρ_{f} (kg/m^{3})
. | ν_{f} (m^{2}/s)
. | U_{in} (m/s)
. | Re
. |

A | 1.21 | 1.51 × 10^{−5} | 1.1 | 29.1 |

B | 2.1 | 55.6 |

. | Density, . | Kinetic viscosity, . | Inlet velocity, . | Reynolds number . |
---|---|---|---|---|

Case . | ρ_{f} (kg/m^{3})
. | ν_{f} (m^{2}/s)
. | U_{in} (m/s)
. | Re
. |

A | 1.21 | 1.51 × 10^{−5} | 1.1 | 29.1 |

B | 2.1 | 55.6 |

. | Density, . | Young’s modulus, . | Poisson’s ratio, . | Length, . | Diameter, . |
---|---|---|---|---|---|

Case . | ρ_{s} (kg/m^{3})
. | E (Pa)
. | ν[−]
. | $L\u0302$ (mm) . | D (mm)
. |

A | 1000 | 2.25 × 10^{5} | 0.01 | 25 | 0.4 |

B |

. | Density, . | Young’s modulus, . | Poisson’s ratio, . | Length, . | Diameter, . |
---|---|---|---|---|---|

Case . | ρ_{s} (kg/m^{3})
. | E (Pa)
. | ν[−]
. | $L\u0302$ (mm) . | D (mm)
. |

A | 1000 | 2.25 × 10^{5} | 0.01 | 25 | 0.4 |

B |

The Reynolds numbers are based on the inlet velocity *U*_{in}, 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 *L*_{x} × *L*_{y} × *L*_{z} = 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.

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.

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.

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.

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.

Density, . | Kinetic viscosity, . | Inlet velocity, . | Reynolds number, . |
---|---|---|---|

ρ_{f} (kg/m^{3})
. | ν_{f} (m^{2}/s)
. | U_{in} (m/s)
. | Re[−]
. |

1.23 | 1.56 × 10^{−5} | 3.5 | 673 |

Density, . | Kinetic viscosity, . | Inlet velocity, . | Reynolds number, . |
---|---|---|---|

ρ_{f} (kg/m^{3})
. | ν_{f} (m^{2}/s)
. | U_{in} (m/s)
. | Re[−]
. |

1.23 | 1.56 × 10^{−5} | 3.5 | 673 |

Density, . | Young’s modulus, . | Poisson’s ratio, . | Length, . | Diameter, . | Rotational stiffness, . |
---|---|---|---|---|---|

ρ_{s} (kg/m^{3})
. | E (MPa)
. | ν[−]
. | $L\u0302$ (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/m^{3})
. | E (MPa)
. | ν[−]
. | $L\u0302$ (mm) . | D (mm)
. | κ_{θ} (N mm)
. |

950 | 1.78 | 0.48 | 70 | 3.0 | 0.41 |

Reynolds numbers are based on the inlet velocity *U*_{in}, 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 *L*_{x} × *L*_{y} × *L*_{z} = 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.

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.

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.05*D* and 0.17*D*, respectively. Compared with the experimental values of 1.04*D* and 0.3*D*, 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.

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.

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.

## V. CONCLUDING SUMMARY

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.

## ACKNOWLEDGMENTS

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.

## AUTHORS’ CONTRIBUTIONS

All authors contributed equally to this work.

## DATA AVAILABILITY

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