Study of flapping filaments using the immersed boundary-lattice Boltzmann method

Abstract

In this study, flow over a flexible filament under a wide range of parameters is simulated using the immersed boundary-lattice Boltzmann method (IB-LBM). The leading end of the filament is fixed in the flow field and the trailing end is free to flap. To execute the simulation, we combine the IB-LBM and a semi-implicit discrete equation of force on the filament to better satisfy the boundary condition. After some numerical simulations validating the modified method, the motion of flexible filaments is examined with different dimensionless bending coefficients ( $κ^{*}$ ) and Reynolds numbers ( $Re$ ). From the trajectory of the flapping filament, different flapping modes are found. When the parameter is between that of two modes, the anti-resonance phenomenon is observed. Numerical results show that the dimensionless bending coefficient and Reynolds number both affect the flapping motion, but in different ways. The dimensionless bending coefficient mainly affects the mode of the flapping, while the Reynolds number mainly affects the perturbation to the filament motion, which is related to the motivation of this system. Some other motion characteristics, for example, the function of amplitude and perturbation propagation, are also discussed in this work.

Keywords

fluid–structure interaction immersed boundary-lattice Boltzmann method flexible filament numerical simulation

Fluid–structure interaction (FSI) problems can be observed frequently in biology kinematics, practical engineering application and many other fields. In biology kinematics, many motions involve complicated interactions between deformable bodies and their surrounding viscous incompressible fluids, like the typical phenomena observed in fish swimming. In mechanical engineering, the design of pumps, aircrafts, artificial heart valves and many other devices encounter FSI phenomena frequently. Sophisticated understanding of the involved FSI problems plays a key role in successfully designing equipment with excellent performance.

Experiments have been carried out to study FSI phenomena by investigating a flexible filament in a flowing soap film.¹ In the experiments, two motion states of the flexible filament have been captured: a stretched-straight state and a flapping state. Filament length is the only changeable coefficient in the experiments. It is found that when the filament length exceeds a critical value, the filament motion state will change from the stretched-straight state to the flapping state. Meanwhile, a bifurcation phenomenon is captured in the experiments. In the sub-critical bifurcation range, filament motion can be stable in either the stretched-straight state or the flapping state. Up to now, many efforts have been made regarding numerical simulation of flexible filaments. Mainly, the immersed boundary method (IBM) is used to simulate this problem.^2–6 Numerical simulations match well with the experimental results. In 2002, Zhu and Peskin² first studied this phenomenon numerically and captured both the flapping state and the stretched-straight state when the Reynolds number was about 200. Further investigations on drag and lift coefficients, trajectory, flow field and so on have achieved many valuable results. The bending coefficient, filament length and mass of the filament are all found to be important to the filament motion.^5–7 The interaction between multiple filaments at different positions has also been investigated.³ Using an attached flexible filament to control flow over the circular cylinder is a highlighted work in recent numerical simulations. Some other research has also studied the motion of a flexible filament in shear flow.⁸

Since the pioneering work by Qian et al.,⁹ the lattice Boltzmann method (LBM) with the DdQb model has been widely adopted in computational fluid dynamics (CFD).^8–12 This model can deal with the d(1,2,3)-dimensional problem with b discrete velocity vectors. Different from other traditional CFD methods, the LBM solves the Boltzmann equation and usually uses Eulerian lattice nodes and shows advantages in being mesh-free, parallel computing and so on. Resonantly, Xu et al.^10,11 optimized the LBM for computational aeroacoustics and discussed the parameters in this method. More and more works have been done to understand the bifurcation phenomenon by the LBM.¹² The direct-force immersed boundary-lattice Boltzmann method (IB-LBM) proposed by Feng and Michaelides^13,14 turns out to be a convenient method to solve the FSI problem, especially for deformable structures. In this paper, we combine the IB-LBM and a semi-implicit discrete equation of force on filaments to better simulate a flexible filament in a uniform flow. Considering the motion of the filament, the modified IB-LBM achieves higher accuracy and stability in simulation.

The rest of this paper is organized as follows: the second section describes the governing equation and the modified IB-LBM. After some numerical validations in the third section, the numerical results and a detailed discussion are presented in the fourth section. Finally, a brief conclusion is given in the fifth section.

Numerical model and method

The numerical method used in this work is different from others' researches, not only the governing equation of the flexible filament, but also the IB-LBM. In this section, the governing equation of the flexible filament used in this paper and the basic idea of the combined treatment of the IB-LBM and the kinematic equation of the filament will be introduced.

Governing equation of the filament

Different from previous works,^2–4 we take the large deflection into consideration. In the previous work, the bending moment is calculated by the second partial derivative of deflection. This is derived from the small deflection assumption. However, the flexible filament flapping in the uniform flow is not small deflection, especially when the Reynolds number is large and the bending coefficient is low. The governing equation of a massive and inextensible filament is written as

ρ_{l} \frac{\partial^{2} X_{α}}{\partial t^{2}} - \frac{\partial}{\partial s} (T n_{α}) + \frac{\partial}{\partial s} [\frac{\partial (κ φ)}{\partial s} τ_{α}] = F_{α}^{fluid}

(1)

where

X_{α}

is the position vector of the filament, s is the Lagrangian coordinate along the filament length,

ρ_{l}

denotes the linear density difference between the filament and the surrounding fluid,³ T is the tension force along the filament, κ is the bending coefficient, φ is the curvature, which reads

φ = \frac{\partial θ}{\partial s}

, and

n_{α}

and

τ_{α}

are the unit tangential vector and the normal vector, respectively. The filament simulated in this work is assumed to be inextensible, that is,

T > > 1

. In previous works,⁶ the dimensionless stretching coefficient is about seven orders of magnitude higher than the dimensionless bending coefficient, which satisfies the assumption. The inextensibility condition is expressed by

| \frac{\partial X_{α}}{\partial s} | = 1

(2)

The upstream end of the filament is fixed in a uniform flow field. At the fixed end of the filament, a simply supported condition is used

X_{α} = X_{c α}, \begin{matrix}  \end{matrix} φ = 0

(3)

At the free end, the condition reads

φ = 0, \frac{\partial φ}{\partial s} = 0

(4)

which means no bending momentum and shear force at the free end.

The basic idea of combined treatment

The IB-LBM uses the LBM to simulate the flow field and the IBM to couple the motion of the solid boundary and the flow field. In the following sections, the LBM and the IB-LBM will be introduced.

Lattice Boltzmann method

The evolution equation for the LBM is^9,15,16

f_{i} (x_{α} + c_{i α} δ_{t}, t + δ_{t}) - f_{i} (x_{α}, t) = - ω [f_{i} (x_{α}, t) - f_{i}^{(eq)} (x_{α}, t)] + F_{i}

(5)

where

f_{i}

is the density distribution function,

f_{i}^{(eq)}

is the equilibrium function,

c_{i α}

is the ith discretized velocity vector and ω is the relaxation parameter. For the D2Q9 model,⁹ the equilibrium function

f_{i}^{(eq)}

f_{i}^{(eq)} (x_{α}, t) = w_{i} ρ [1 + \frac{c_{i α} v_{α}}{c_{s}^{2}} + \frac{v_{α} v_{β}}{2 c_{s}^{2}} (\frac{c_{i α} c_{i β}}{c_{s}^{2}} - δ_{α β})]

(6)

where the sound speed

c_{s}

and the weight coefficient

w_{i}

are

c_{s} = \frac{1}{\sqrt{3}} \frac{δ_{x}}{δ_{t}}, \begin{matrix} w_{i} = {\begin{matrix} 4 / 499 & i = 0 \\ 1 / 199 & i = 1 \sim 4 \\ 1 / 13636 & i = 5 \sim 8 \end{matrix} \end{matrix}

(7)

The relation between the relaxation parameter ω and kinematic viscosity ν is

ν = c_{s}^{2} (\frac{1}{ω} - \frac{1}{2}) δ_{t}

(8)

The density (ρ) and the momentum of fluid ( $ρ u_{α}$ ) can be calculated by $ρ = \sum_{i} f_{i}$ and $ρ u_{α} = \sum_{i} f_{i} c_{i α}$ , respectively. We adopt the force scheme proposed by He et al.¹⁶

F_{i} = (1 - \frac{ω}{2}) \frac{(c_{i α} - u_{α}) F_{α}}{ρ c_{s}^{2}} f_{i}^{(eq)}

(9)

where

F_{α}

is the external force.

Taking the external force into consideration, the momentum should be estimated as $ρ u_{α} = \sum_{i} f_{i} c_{i α} + \frac{1}{2} F_{α} δ_{t}$ .

IB-LBM

The direct-force model proposed by Feng and Michaelides¹⁴ achieves excellent results in simulating the fluid–particle interaction problem. It can be summarized in the following four steps.

1) Gain the velocity of every point set on the boundary of the particles. In this step, the delta function is used to gain the velocity of particles from the flow field. Feng et al.¹⁴ suggest

\begin{matrix} δ (r) = {\begin{matrix} \frac{1}{4} [1 + cos (\frac{π | r |}{2})] \begin{matrix} | r | < 2 \end{matrix} \\ 0 \begin{matrix}  \end{matrix} | r | \geq 2 \end{matrix} \end{matrix}

(10)

In two-dimensional simulation, the delta function

Δ (x_{0 α})

is defined as

Δ (x_{0 α}) = δ (x - x_{0}) δ (y - y_{0})

, where

(x, y)

represents the Cartesian coordinates of the lattice node and

(x_{i}, y_{i})

represents the Cartesian coordinates of the boundary points of the particles. Thus, the velocity of the unforced boundary points

u (x_{b α})

can be estimated as

u_{α} (x_{b α}) = \int u_{α} (x_{α}) Δ (x_{α}^{b}) d s

(11)

2) Calculate the restoration force. In this step, the restoration force aims to change the velocity of the flow around the boundary points of the particles in order to satisfy the boundary conditions.

3) Apply the restoration force on the boundary points of the particles into the flow field. The restoration force is stored in the Lagrangian boundary points. So, the delta function is used to transfer the Eulerian lattice nodes to the body force

F_{α} (x_{α}) = \int F_{α} (x_{α}^{b}) Δ (x_{α}^{b}) d l

(12)

4) Update the position of the particle points and the flow field with the standard LBM model. Then go back to step 1.

The key steps are steps 1–3. The correct restoration force is important to the accuracy of simulation. Feng and Michaelides¹⁴ suggest estimating the restoration force by following the direct-force term

F_{α} = ρ \frac{u_{s α} - u_{α}}{Δ t} + \partial_{β} (u_{α} u_{β} - μ \partial_{β} u_{α}) + \partial_{α} p

(13)

where

u_{s α}

is the velocity of the boundary points.

Discrete equations

Considering the previous discussed equations (Equations (1) and (13)), the discrete form of the kinematic equation for the lth nodal point is written as

\begin{matrix} ρ_{s} \frac{d^{2} X_{α}}{d t^{2}} = F^{total} = F^{inner} + F_{α}^{fluid} \\ = \frac{\partial}{\partial s} (T n_{α}) - \frac{\partial}{\partial s} (\frac{\partial (κ φ)}{\partial s} τ_{α}) + \frac{ρ}{Δ t} (u_{s α} - u_{α}) \\ = fixl (X_{α}^{N + 1}) - \frac{κ [(φ_{l + 1}^{N} - φ_{l}^{N}) τ_{(l + 1 / 2) α} - (φ_{l}^{N} - φ_{l - 1}^{N}) τ_{(l - 1 / 2) α}]}{Δ s^{2}} \\ + \frac{ρ^{N}}{Δ t} (\frac{X_{(l) α}^{N + 1} - X_{(l) α}^{N - 1}}{2 Δ t} - u_{α}) + vispre (u_{α}^{N}, ρ^{N}) \end{matrix}

(14)

where N means the Nth iteration step.

fixl (X_{α}^{N})

is determined by the inextensibility condition. Noticing that the tension force is along the filament, the effect of this function deforms the filament point along the tangential direction. The curvature

φ_{l}^{N}

is discrete

\begin{matrix} φ_{l}^{N} = \frac{2 \arcsin (\frac{| δ x_{(l + 1 / 2)} δ y_{(l - 1 / 2)} - δ x_{(l - 1 / 2)} δ y_{(l + 1 / 2)} |}{\sqrt{δ x_{(l + 1 / 2)}^{2} + δ y_{(l + 1 / 2)}^{2}} \sqrt{δ x_{(l - 1 / 2)}^{2} + δ y_{(l - 1 / 2)}^{2}}})}{\sqrt{δ x_{(l + 1 / 2)}^{2} + δ y_{(l + 1 / 2)}^{2}} + \sqrt{δ x_{(l - 1 / 2)}^{2} + δ y_{(l - 1 / 2)}^{2}}} \end{matrix}

(15)

where

δ x_{(l + 1 / 2)} = X_{(l + 1) x}^{N} - X_{(l) x}^{N}

and

δ y_{(l + 1 / 2)} = X_{(l + 1) y}^{N} - X_{(l) y}^{N}

To solve the discrete equation, the semi-implicit equation (Equation (14)) can be rewritten as below

\begin{matrix} (\frac{ρ_{s}}{Δ t^{2}} - \frac{ρ^{N}}{2 Δ t^{2}}) X_{(l) α}^{N + 1} - fixl (X_{α}^{N + 1}) \\ = - \frac{κ [(φ_{l + 1}^{N} - φ_{l}^{N}) τ_{(l + 1 / 2) α} - (φ_{l}^{N} - φ_{l - 1}^{N}) τ_{(l - 1 / 2) α}]}{Δ s^{2}} \\ - \frac{ρ^{N}}{Δ t} (\frac{X_{(l) α}^{N - 1}}{2 Δ t} + u_{α}) - ρ_{s} \frac{X_{(l) α}^{N - 1} - 2 X_{(l) α}^{N}}{Δ t^{2}} + vispre (u_{α}^{N}, ρ^{N}) \end{matrix}

(16)

The numerical implementation of the present combined model in each iteration step can be illustrated as follows.

(1) Gain the velocity of the fluid field at each the position of each Lagrangian point, that is, calculate the term $- \frac{ρ^{N}}{Δ t} u_{α}$ .

(2) Compute the previous force term using the position of the previous Lagrangian point, that is, calculate the term $- ρ^{N} \frac{X_{(l) α}^{N - 1}}{2 Δ t^{2}} - ρ_{s} \frac{X_{(l) α}^{N - 1}}{Δ t^{2}}$ .

(3) Compute the bending force term using the curvature, which is calculated with the position of the present Lagrangian point, that is, calculate the term $- \frac{κ [(φ_{l + 1}^{N} - φ_{l}^{N}) τ_{(l + 1 / 2) α} - (φ_{l}^{N} - φ_{l - 1}^{N}) τ_{(l - 1 / 2) α}]}{Δ s^{2}}$ . Compute the viscosity and pressure term: $vispre (u_{α}^{N}, ρ^{N})$ .

(4) Finally, determine the new position of the Lagrangian point. After divided by the coefficient ahead, $fixl (X_{α}^{N + 1})$ , change the position along the tangential direction to keep the distance of every two points the same as that of the initial condition.

Numerical validation

To validate the present method and our code, simulations on the flow over a circular cylinder and flow over a cylinder with a flexible filament attached in a uniform flow are simulated.

The flow over a circular cylinder has been studied extensively through theoretical, experimental and numerical methods. In the simulation, the computational domain is

20 D \times 20 D

and

D = 20 (l . u .)

(lattice unit). Three different Reynolds numbers based on the diameter D and the inflow speed

U_{\infty}

, namely

Re = 20, 40

and 100, are selected. The drag coefficient

C_{d} = F_{d} / (0.5 ρ U_{\infty}^{2} D)

and separation bubble length (

L_{w}

) of each steady case are listed in Table 1 and compared with other results. Table 1 also shows the time-average drag coefficient and the Strouhal number

St = fD / U_{\infty}

from the simulation. From the results shown in the table, it is clear that the current method can perfectly handle simulations of flow with curved boundaries.

Table 1.

Comparison between other numerical methods and the present method. The table demonstrates the drag coefficient ( $C_{d}$ ), the separation bubble length ( $L_{w}$ ) of flows over a cylinder at Re = 20,40 and the time-average drag coefficient, the Strouhal number for the Re = 100 case using different methods

	$Re = 20$		$Re = 40$		$Re = 100$
Ref.	$C_{d}$	$L_{w}$	$C_{d}$	$L_{w}$	Avg. $C_{d}$	$St$
Niu et al.¹⁷	2.144	1.89	1.589	4.52	−	−
Wu and Shu¹⁸	2.091	1.86	1.565	4.62	1.364	0.163
Kang Zhu and Hassan¹⁹	2.057	1.82	1.538	4.50	1.336	0.165
Park et al.²⁰	2.01	−	1.51	−	1.33	0.165
Present	2.022	1.87	1.495	4.67	1.332	0.166

Wu et al.⁵ studied a flapping filament in the wake of a cylinder with its upstream fixed on the cylinder at low Reynolds number in detail. To further validate the present method for simulating elastic bodies, we simulated this flow and compared it with Wu et al.'s work.⁵ In the simulations, the Reynolds number $Re$ is 100 and the dimensionless filament lengths $L^{*} (= L / D)$ are 1 and 2, respectively. The dimensionless bending coefficients $κ^{*} (= κ / (ρ_{f} U_{\infty}^{2} L^{3}))$ are $1 \times 10^{- 2}$ , $1.5 \times 10^{- 2}$ , $2 \times 10^{- 2}$ and $2.5 \times 10^{- 2}$ .

The results shown in Figure 1 manifest the asymmetrical motion of the filament when the dimensionless filament length equals 1. When the dimensionless filament length equals 2 and the bending coefficient is low, the motion of the filament becomes symmetric. The trajectory of the free end emerges into a pattern shaped like a butterfly. These results qualitatively match with the work of Wu et al.⁵ and validate the reliability of our simulation.

Figure 1.

The trajectory of the free end of the filament: (a) $L^{*} = 1$ ; (b) $L^{*} = 2$ .

Numerical simulations and discussion

In this section, a flexible filament in a uniform flow with its upstream end fixed in the flow field is simulated. It is reported in a laboratory experiment¹ that the motion of a filament in a uniform flow has two different states, that is, the flapping state and the stretched-straight state. The only adjustable parameter is restricted to the length of filament; thus, the obtained results merely show the effects of the filament length on the motion state. One of the advantages of using the numerical method is a wide choice of adjustable parameters compared to laboratory experiments. Various different characters of the flexible filament affect its motion. According to dimensional analysis, we focus on two main dimensionless parameters (i.e. the dimensionless bending coefficient $κ^{*}$ and the Reynolds number $Re$ ) that determine the filament motion state. The filament used in our simulation is represented by 50 Lagrangian points and the filament length $L = 100 (l . u .)$ . The computational domain is $400 \times 600 (l . u .)$ . Both stable and flapping motion states are successfully captured in our simulations.

Effect of the bending coefficient

In this section, the effect of the dimensionless bending coefficient on the filament motion is investigated. The dimensionless bending coefficient $κ^{*}$ is defined as mentioned previously. To study the effect of the bending coefficient, the other parameters are fixed as follows: Reynolds number $Re = 1000$ ; dimensionless linear density difference $ρ_{s}^{*} = \frac{ρ_{s}}{ρ_{f} L} = 0.01$ .

The relation between the amplitude of the free end spanwise and the dimensionless bending coefficient is depicted in Figure 2. It is found that the amplitude of filament tends to decrease with the increase of nondimensionless bending coefficient. When the dimensionless bending coefficient exceeds a critical number, the amplitude decreases to zero and the motion state changes from the flapping state to the stretched-straight state. The critical bending coefficient $κ^{*}$ is between $1.65 \times 10^{- 3}$ and $1.8 \times 10^{- 3}$ , when $Re = 1000$ .

Figure 2.

The amplitude of the filament varies with the dimensionless bending coefficient.

Figure 2 also shows the existence of some special dimensionless bending coefficients, for example, $κ^{*} = 1.5 \times 10^{- 4}$ , $κ^{*} = 3.25 \times 10^{- 4}$ and $κ^{*} = 8.25 \times 10^{- 4}$ . The amplitude of the free end at these bending coefficients is much lower than those at nearby ones. These three special points divide the flapping interval into three intervals. Three cases from these three intervals, namely, .., are simulated and their trajectories are shown in Figure 3. With the increment of the bending coefficient, the wave number decreases. Noticing that the structure is an analogous standing wave, the wave numbers are 2.5, 2 and 1.5 for $κ^{*} = 3 \times 10^{- 4}, 6 \times 10^{- 4}$ and $1.2 \times 10^{- 3}$ , respectively. The flapping is the resonance at these bending coefficients. Consequently, the flapping is anti-resonance at the special coefficients where the amplitude of the free end is lower than those of nearby ones. Consequently, the bending coefficient mainly affects the physical properties of the filament, especially the natural frequency of the system.

Figure 3.

Flapping modes of the filament flapping in the flow with different bending coefficients: (a) $κ^{*} = 3 \times 10^{- 4}$ ; (b) $κ^{*} = 6 \times 10^{- 4}$ ; (c) $κ^{*} = 1.2 \times 10^{- 3}$ ( $Re = 1000$ ).

It is an interesting result that the trajectory of the filament is analogous to a standing wave at the fixed end, but is analogous to a traveling wave at the free end. Figure 4 shows the transformed trajectory of the filament. According to Zhang et al.,¹ the trajectory of the filament can be expressed as the following formula

Y (x, t) = f (x) sin (kx - ω t)

(17)

where

f (x)

is a quadratic function. However, the above equation is not suitable for describing the trajectory of the filament obtained by our simulation. According to the previous analysis, we proposed the following formula instead

\begin{matrix} Y (x, t) = f (x) sin (kx - ω t) + g (x) sin (kx + ω t) \\ = [f (x) + g (x)] sin (kx) cos (ω t) \\ - [f (x) - g (x)] cos (kx) sin (ω t) \end{matrix}

(18)

where

f (x)

and

g (x)

are the amplitudes of the incident wave and the reflected wave, respectively. Considering the structure of the trajectory obtained in our simulations,

f (x)

and

g (x)

should satisfy

{lim}_{x \to 0} \frac{g (x)}{f (x)} = 1, \begin{matrix}  \end{matrix} {lim}_{x \to L} \frac{g (x)}{f (x)} = 0

(19)

Figure 4.

The trajectory of the filament transformed with an exponential function. The dimensionless bending coefficient is $3 \times 10^{- 4}$ . The exponential function is $0.0023 L exp (3.6 x / L)$ .

In the case shown in Figure 4, $f (x) + g (x)$ is approximately an exponential function. However, in other cases, the exponential function rule dissatisfies the trajectory; at least, the quadratic function is not the rule of any amplitude function. Considering the relation between $f (x)$ and $g (x)$ shown in Equation (19), it is found that perturbation starts from the free end, propagates to the fixed end, reflects at the fixed end and then propagates back to the free end. The energy of the wave keeps increasing in both wave propagation processes.

Effect of the Reynolds number

The Reynolds number is a typical coefficient in fluid dynamics. It presents the ratio of inertial forces to viscous forces. In the simulations of flow over a flexible filament, the Reynolds number based on the filament length L is another crucial parameter. In this section, the filament length, bending coefficient and inflow velocity are fixed in all cases. The Reynolds number changes with the viscosity of the fluid. To study the effect of the Reynolds number, two typical cases, $κ^{*} = 3 \times 10^{- 4}$ and $κ^{*} = 6 \times 10^{- 4}$ , are simulated. The relation between the amplitude of the free end A and the Reynolds number $Re$ is shown in Figure 5.

Figure 5.

The amplitude of the free end of the filament varies with the Reynolds number at a fixed dimensionless bending coefficient: (a) $κ^{*} = 3 \times 10^{- 4}$ ; (b) $κ^{*} = 6 \times 10^{- 4}$ .

When the Reynolds number is low, the flow is stable and the perturbation is small. With the increment of the Reynolds number, the perturbation increases and the filament starts flapping. It is obtained that the critical Reynolds number is a little larger than 375 for $κ^{*} = 3 \times 10^{- 4}$ and a little larger than 500 for $κ^{*} = 6 \times 10^{- 4}$ . When the Reynolds number increases after the critical number, the amplitude of the free end increases rapidly, but rarely changes when the Reynolds number further increases. The amplitude even decreases gradually with the increase of the Reynolds number when $Re > 2000$ for $κ^{*} = 6 \times 10^{- 4}$ . There is no anti-resonance phenomenon when $κ^{*} = 6 \times 10^{- 4}$ . When $κ^{*} = 3 \times 10^{- 4}$ , the anti-resonance phenomenon is observed only once. The mode changes at the anti-resonance point, and the wave number changes from 2 to 2.5. From Figure 6, it is found that though wave numbers of the filament trajectory are the same, but the amplitude functions $f (x)$ and $g (x)$ are different. The Reynolds number is high, and the ratio of $g (x)$ to $f (x)$ is high. This means that when the Reynolds number is high, the wave obtained more energy in the incident process than when the Reynolds number is low. To further study the motivation of the system, vortex structures under different Reynolds numbers are depicted in Figure 7. It can be seen from Figure 7 that the vortex sheds when the Reynolds number is beyond the critical number. When the Reynolds number is low, the structure of the shedding vortex resembles a Kármán vortex street. With the increment of the Reynolds number, the Strouhal number ( $St$ ) is around 0.2 and rarely changes in the interval of the Reynolds number in our simulations. The magnitude of the vortex increases as the Reynolds number increases. Considering the analysis of Figure 6, the Reynolds number mainly determines the motivation of the system, especially the magnitude of motivation.

Figure 6.

The trajectory of the flapping filament with the same dimensionless bending coefficient ( $κ^{*} = 6 \times 10^{- 4}$ ) at different Reynolds numbers.

Figure 7.

The vortex structure of the flapping filament with the same dimensionless bending coefficient ( $κ^{*} = 6 \times 10^{- 4}$ ) at different Reynolds numbers: (a) $Re = 500$ ; (b) $Re = 1000$ ; (c) $Re = 2000$ ; (d) $Re = 5000$ .

Conclusions

In this paper, a combined method is used to simulate the flexible filament in a uniform flow. The motion state of the filament is investigated by varying two dimensionless coefficients: the dimensionless bending coefficient ( $κ^{*}$ ) and the Reynolds number ( $Re$ ). Simulations are performed in a wide range of these two parameters to better understand the physical mechanism of the flow. Based on the numerical results, it is found that these two dimensionless coefficients both affect the flapping of the filament but in different ways.

For the bending coefficient, it determines the physical properties of the filament and it mainly affects the inherent parameters, such as the natural frequency of this vibration system. When increasing the bending coefficient and the Reynolds number is fixed, the wave number on the filament decreases. The dimensionless bending coefficient also affects the amplitude of the free end. When it increases, the amplitude decreases tendentiously and, finally, the motion state turns from the flapping state to the stretched-straight state.

For the Reynolds number, it determines the properties of flow and it mainly affects the motivation of the flapping system. When the Reynolds number is lower than the critical number ${Re}_{c}$ , the magnitude of motivation is not large enough to active the flapping state and the motion is in the stretched-straight state. When the Reynolds number exceeds the critical number, the filament motion changes from the stretched-straight state to the flapping state. The critical number is related to the bending coefficient of the filament. It can be expressed as a function of $κ^{*}$ . There is believed to be a positive correlation between the critical Reynolds number $Re$ and the dimensionless bending coefficient $κ^{*}$ . Different from the bending coefficient, the amplitude rapidly increases after the critical number and rarely changes when the Reynolds number further increases. The ratio of the amplitude function of the incident wave to that of the reflected wave increases with the increment of the Reynolds number.

Moreover, the analyses of trajectories of the flapping filaments show that the wave can be decomposed into an incident wave and a reflected wave, perturbation starts from the free end, propagates to the fixed end and then reflects and, finally, propagates back to the free end. The energy of the wave keeps increasing during the entire propagation process.

Footnotes

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (No. 11872337 and No. 11502237), the Natural Science Foundation of Zhejiang Province (No. LY18A020010) and the Project of Ten Thousands Talents (No. W01060069).

Appendix

References

Zhang

Childress

Libchaber

et al.

Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind. Nature 2000; 408: 835–839.

Zhu

Peskin

. Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method. J Comput Phys 2002; 179: 452–468.

Huang

Shin

Sung

. Simulation of flexible filaments in a uniform flow by the immersed boundary method. J Comput Phys 2007; 226: 2206–2228.

Tian

Luo

Zhu

. Interaction between a flexible filament and a downstream rigid body. Phys Rev E 2010; 82: 026301.

Qiu

Shu

et al.

Flow control of a circular cylinder by using an attached flexible filament. Phys Fluid 2014; 26: 103601.

Yuan

Niu

Shu

et al.

A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating a flexible filament in an incompressible flow. Comput Math Appl 2014; 67: 1039–1056.

Connell

Yue

. Flapping dynamics of a flag in a uniform stream. J Fluid Mech 2014; 67: 1039–1056.

. A dlm/fd method for fluid/flexible-body interactions. J Comput Phys 2005; 207: 1–27.

Qian

d'Humières

Lallemand

. Lattice BGK models for Navier-Stokes equation. Europhys Lett 1992; 17: 479.

10.

Sagaut

. Optimal low-dispersion low-dissipation LBM schemes for computational aeroacoustics. J Comput Phys 2011; 230: 5353–5382.

11.

Malaspinas

Sagaut

. Sensitivity analysis and determination of free relaxation parameters for the weakly-compressible mrtlbm schemes. J Comput Phys 2012; 231: 7335–7367.

12.

Shu

et al.

Study of multiple steady solutions for the 2D natural convection in a concentric horizontal annulus with a constant 265 heat flux wall using immersed boundary-lattice Boltzmann method. Int J Heat Mass Transf 2015; 81: 591–601.

13.

Feng

Michaelides

. The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems. J Comput Phys 2004; 195: 602–628.

14.

Feng

Michaelides

. Proteus: a direct forcing method in the simulations of particulate flows. J Comput Phys 2005; 202: 20–51.

15.

Bhatnagar

Gross

Krook

. A model for collision processes in gases. I: Small amplitude processes in charged and neutral one-component systems. Phys Rev 1954; 94: 511–525.

16.

Shan

Doolen

. Discrete Boltzmann equation model for nonideal gases. Phys Rev E 1998; 57: R13–R16.

17.

Niu

Shu

Chew

et al.

A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows. Phys Lett A 2006; 354: 173–182.

18.

Shu

. Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications. J Comput Phys 2009; 228: 1963–1979.

19.

Kang Zhu

Hassan

. A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries. Int J Numer Meth Fluid 2011; 66: 1132–1158.

20.

Park

Kwon

Choi

. Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160. KSME Int J 1998; 12: 1200–1205.