Coupling Motion and Energy Harvesting of Two Side-by-side Flexible Plates in a 3d Uniform Flow

The fluid-structure interaction problems of two side-by-side flexible plates with a finite aspect ratio in a three-dimensional (3D) uniform flow are numerically studied. The plates' motions are entirely passive under the force of surrounding fluid. By changing the aspect ratio and transverse distance, the coupling motions, drag force and energy capture performance are analyzed. The mechanisms underlying the plates' motion and flow characteristics are discussed systematically. The adopted algorithm is verified and validated by the simulation of flow past a square flexible plate. The results show that the plate's passive flapping behavior contains transverse and spanwise deformation, and the flapping amplitude is proportional to the aspect ratio. In the side-by-side configuration, three distinct coupling modes of the plates' motion are identified, including single-plate mode, symmetrical flapping mode and decoupled mode. The plate with a lower aspect ratio may suffer less drag force and capture less bending energy than in the isolated situation. The optimized selection for obtaining higher energy conversion efficiency is the plate flapping in single-plate mode, especially the plate with a higher aspect ratio. The findings of this work provide several new physical insights into the understanding of fish schooling and are expected to inspire the developments of underwater robots or energy harvesters.


Introduction
Interactions between flexible structures and viscous flow are ubiquitous and familiar in nature and engineering applications.The physical mechanisms underlying these interactions are of interest because they provide us novel insights into the fluid dynamics in biological growth or locomotion.In natural creatures, the flexibility of the structure plays an important role to improve its body stability or propulsion performance in surrounding flow.For instance, through its passive shape reconfiguration in the flow, a tree leaf can achieve substantial and beneficial drag reduction [1,2].Similarly, a fish can exploit energy from the surrounding vortex street and achieve cost savings by undulating its body [3][4][5].The famous experiment conducted by Liao et al. revealed that a trout can decrease muscle activity by synchronizing its kinematics to incoming vortices appropriately [6].Further, a dead fish was observed to propel upstream when its flexible body resonated with the wake of a bluff cylinder [7]; this entirely passive locomotion without energy input was also found in a simple articulated fish-like system [8].Regarding the applications, the underlying hydrodynamic principles also inspire researchers to develop advanced equipment involving the interactions of unsteady flows and deformable bodies.Recently, the robotic models of fish swimming have been of great interest in ocean engineering, and a wide diversity of approaches have been taken for the mechanical design of fish-inspired systems [9][10][11], such as a carangiform fish robot [12,13] and a batoid-inspired robot [14][15][16].Besides the active flow control mechanisms used in bionic propulsion, passive flapping or vibrating dynamics have also been applied in developing renewable energy harvesters, which are usually based on compliant materials, such as elastic-mounted cylinders and piezoelectric membranes [17][18][19].Although these interaction phenomena are common, they still puzzle researchers due to their complicated dynamics.
To obtain physical insights into the underlying mechanisms and dynamics of the flexible structure-fluid interactions, flexible filaments or plates are usually adopted due to their simplicities and usefulness [20,21].The deformation of a flexible filament or plate in a motion flow is physically determined by the competition between the fluid load on it and internal elasticity.During the past few years, these challenging problems have attracted many researches.Tenada first conducted experiments on a flag in a wind tunnel and observed the oscillation modes of the waving motions [20].The existence of bistability for an isolated filament, including stretched-straight mode and self-sustained flapping mode, was observed in Zhu and Peskin's research [22], where they carried out a series of computer simulations by an immersed boundary method.Zhu and Peskin's research further revealed that the sustained flapping only occurred when the filament mass was included, and the bistability depended on the filament's length and initial condition [22].Moreover, three distinct regimes, including the straight, regular flapping and chaotic flapping regimes, have been observed by changing the mass ratio, bending rigidity or incoming velocity [23][24][25].Additionally, flow past multiple filaments or plates arranged side-by-side is also an important topic for understanding the mechanism of fish schooling or birds flocking.Both in-phase, out-of-phase and erratic flapping interaction modes have been found by experiments and numerical simulations [26][27][28].More recently, it was found that schooling fish can improve locomotion efficiency in side-by-side configuration [29,30].Although most of the flexible bodies physically exist with a finite aspect ratio, only a few 3D simulations of an isolated flexible plate have been reported [31][32][33][34].To our knowledge, the coupling motions and energy conversion performance of two side-by-side flexible plates in a 3D viscous flow have not been modeled and studied numerically yet.
The primary purpose of this work is to correlate the fluid-induced flapping dynamics and energy conversion performance of the flexible plates with their aspect ratios and transverse gap distances in 3D.Towards this end, we apply a hybrid algorithm combined with the lattice Boltzmann method (LBM) [35,36] and immersed boundary method (IBM) [37] to solve the fluid flow and fluid-structure interactions.To solve the passive deformation of the flexible body, a finite element method based on a triangular mesh [38,39] is adopted.To validate this coupled algorithm, the dynamics of flow past a square flexible plate were simulated at low Reynolds number.Then, numerical simulations based on various parameters were carried out, and the coupled flapping modes and energy conversion performance were analyzed systematically.The rest of this paper is organized as follows.In Section 2, we detail the physical model and numerical method, including the algorithm validation.In Section 3, we perform the simulations and discuss the results.Finally, a brief conclusion is drawn in Section 4.

Physical Problem
A 3D model of the flexible plates in a uniform flow is considered.Figure 1 shows the schematic of flow over two flexible plates in side-by-side arrangement.Both of the two plates are equal in geometry and mechanical property.As shown in Figure 1, the two plates with a fixed leading edge and a free trailing edge are immersed in a uniform flow.The flow's incoming velocity, dynamic viscosity and density are denoted as U ∞ , µ and ρ f , respectively.In this work, the aspect ratio is finite, defined as AR = W/L, where W is the plate's width and L is the plate's length.Furthermore, the plate has a density ρ s , thickness h and Young's modulus E. In the initial status, the upper and lower plates keep parallel in the field and have an inclined angle γ 0 (γ 0 =π/10) with the streamwise plane.The points C 1 and B 1 are located on the mid-point of the upper plate's leading edge and trailing edge, respectively.
The line C 1 B 1 is called the upper plate's mid-line, as well as C 2 B 2 to the lower plate.The transverse gap distance G z is defined as the distance from point C 1 to point C 2 , denoted as G z = |C 1 C 2 |.For the purposes of simplification, the non-dimensional parameters L, W and G z will be used to denote the variants L/L, W/L and G z /L, hereafter.Based on the dimensionless parameters U ∞ , ρ f and L, the 3D flow field with flexible plates can be described as forced incompressible Navier-Stokes equations [31,32,40]: where f s is the body force exerted on the fluid by the structure's boundary.Equation ( 3) is the plate's governing equation, which describes its deformation with large displacement [31,32,40].In this equation, X is the position of the Lagrangian node on the plate in global coordinates XYZ; K b and K s are the bending and stretching stiffness of the plate, respectively.F is the force density acted on the plate's boundary by the surrounding flow; δ ij is the Kronecker δ function; and γ ij and ϕ ij are the out-of-plane and in-plane matrices, respectively.The mentioned above dimensionless variants concerned in the present model are given as follows: Meanwhile, the Reynolds number in this work is defined based on the plate's length and incoming flow and can be expressed as: The drag and lift force are also defined in a dimensionless form, as follows: where F D is the force exerted on the plate in the X-axis direction, while F L is in the Z-axis direction.

Numerical Methods
The governing equations of the viscous flow with flexible plates are solved numerically by a hybrid algorithm coupled with the immersed boundary-lattice Boltzmann method and the finite element method.The lattice Boltzmann is adopted to compute the flow field, while the finite element method is applied to compute the structures' motion.To handle the interactions between plates and surrounding fluid, a four-point-type δ function is introduced.

Mathematical Formulation
Due to its efficiency and the relative simplicity of the computation, the lattice Boltzmann equation (LBE) with the BGK model has been widely employed to simulate a complex flow field.This method acts as an alternative to the conventional numerical methods in solving Navier-Stokes equations.The BGK-based LBM can be expressed as [35,36,41]: where τ is the dimensionless relaxation time, τ = µ/c 2 s + 1/2∆t; ∆t is the time increment and t is the lattice time step.Additionally, f i (x, t) is the distributed function for particles at position x with velocity c i .The lattice equilibrium distributed function f eq i and the forcing function f F i are given as [42]: where ω i is the weighted factor, in the present D3Q19 model [36,43], ω 0 = 1/3, ω 1−6 = 1/18, ω 7−18 = 1/36.Additionally, c s is the lattice speed of sound.The mass density ρ and fluid velocity u can be obtained by: With the Dirac delta function in the immersed boundary model, the Lagrangian momentum forcing term is spread to the Eulerian grids, while the Eulerian velocities are interpolated to the Lagrangian meshes; the numerical schemes can be expressed as [36,37]: Based on our careful validations and convergence examinations with different computational sizes and lattice spacings, the computational domain is modeled as a rectangular box, extending from [0, 0, 0] to L x , L y , L z L x = 8, L y = 4, L z = 12 in the X (streamwise), Y (spanwise) and Z (transverse) directions, respectively.The grid spacing has a uniform value of ∆x = ∆y = ∆z = 0.02.Here, both the domain size and grid spacing are scaled by the plate's length L. To simulate the plates' flap in free stream, the Dirichlet boundary condition (u = U ∞ , v = 0, w = 0) is adopted at the inflow (x = 0), and a convective boundary condition (∂u/∂t + U ∞ ∂u/∂x = 0) is applied at the outflow (x = L x ).Meanwhile, the far-field boundary condition and periodic boundary condition are used in the transverse (z = 0 and z = L z ) and spanwise (y = 0 and y = L y ) direction, respectively.The fixed boundaries of the plates (x b = 0) are aligned with the Y-axis, with their mid-points C 1 and C 2 located in the mid-plane of the computational domain in the Eulerian coordinate system (x = L,y = 1/2L y ).Each of the plates is discretized using a three-node triangular element [38].In our simulations, the parameters of ϕ 11 = ϕ 22 = 10 3 , ϕ 12 = 10, γ 11 = γ 12 = γ 22 = 0.0001, M s = 1.0,K b = 0.0001 and K s = 1000 are also taken for the plates by default.Here, the shearing and tension coefficients are taken to be large enough to satisfy the inextensibility conditions for the plates.For the comparison with previous studies, the gravity force is neglected in the present work.All of the computations are implemented at Re = 200 in this work.

Validation
The validation of the present numerical algorithm will be briefly described in this section.Considering the typical case, we simulate a square (AR = 1) flexible plate in a uniform flow.Compared to the twin plate model shown in Figure 1, there is only an isolated flexible plate immersed in the flow field in this case.Meanwhile, the leading edge of the plate is located in the mid-plane of computational domain (z = 1/2L z ). Figure 2 shows the transverse displacement of the mid-point on the trailing edge; the black solid and red dashed line are the present result and given by Huang and Sung [31], respectively.As it shows, the transverse displacement varies periodically with time increasing, meaning that the plate is flapping passively in the transverse direction.Furthermore, the flapping amplitude of the trailing edge's mid-point is about 0.77.The present result is almost coincident with the previous studies [31][32][33].The good agreement indicates that the present numerical algorithm is accurate.To further analyze the motion of the plate, Figure 3 shows the envelop of an isolated flapping flexible plate, including (a) the 3D envelop of the flapping plate and (b) the envelop of the plate's mid-line.In the plate's 3D flapping envelop, the black solid line denotes its mid-line, the segmentation before and after which are displayed in semi-transparent mesh and solid surface style, respectively.As shown in Figure 3a, the plate's passive flapping is around its fixed leading edge.In the process of passive flapping, the plate loses its initial flat status, which is obviously different from the rigid plate's flapping [44,45].Moreover, the motion of the flexible plate apparently contains the deformation streamwise and the deformation spanwise.Due to its flexibility, the plate's motion in the uniform flow is entirely passive, induced by the interaction between surrounding flow and the structure itself.This deformation characteristic is similar to flapping wings or fins in some extent, especially in spanwise flexibility.According to our previous study [46], the deformation of the tuna's caudal fin spanwise can improve its propulsive efficiency.Compared to imposing on the fin the prescribed motion equation spanwise, here the deformation is the result of fluid and inertial force, which is much closer to fish swimming in nature.In Figure 3b, point B is the mid-point on the pate's trailing-edge.The trajectory of point Blooks like the number "8", which is formed by the free flapping of the plate's trailing-edge in the transverse direction with the longitudinal traveling wave.Meanwhile, the trajectory of point B is symmetrical about the plate's local x b -axis, indicating that passive flapping of the plates keeps symmetry in the positive and negative z b -axis direction.Due to the phenomenon of asymmetrical flapping occurring in the side-by-side arrangement, we define Z B,out and Z B,in as the absolute value of the plate's maximum transverse displacement in the positive and negative z b -axis direction, respectively.Additionally, the flapping amplitude A f of an isolated plate is defined as The hydrodynamics and motion of the flexible plate are closely associated with the vortex in the wake.Figure 4 shows the instantaneous vortical structure of flow over an isolated square flexible plate at Re = 200.Here, the vortical structure is identified by the λ 2 criterion [47], and an iso-surface of λ 2 = −0.2 is taken to plot the vorticity contour in 3D.From top to bottom, the three snapshots are chosen at the same instant, but with different projection directions.It can be clearly seen from the figure that a chain of hairpin-like vortical structures is arrayed in order around the flexible plate.In each of the flapping cycles, the vortex formed in the side edges will connect with that formed in the trailing edge when they are moving downstream by the force of inflow.Hence, the connecting vortex actually contains two parts, named the side edge vortex and the trailing edge vortex.Depending on the relative strength of the side edge vortex and trailing edge vortex, the structure behind a plate will show as a double line, hairpin, or ring shaped.Further, it is inferred that the side edge vortex is the main cause of the plate's deformation in the spanwise direction.Compared to the previous research [31,32], the present study also shows good agreement.In this study, the aspect ratio is changed from 0.5 to 2.0.As shown in Figure 5a, both the mean drag coefficient and the maximum lift coefficient increase almost linearly with the increasing aspect ratio.As regards the increasing slope, the maximum lift coefficient is clearly higher than that of the mean drag coefficient.This indicates that the maximum lift coefficient is much more sensitive to the aspect ratio.Despite both of the force coefficients monotonically increasing with the aspect ratio, the slopes are less than 1.0.This means that the forces on the plate are not able to increase with the increasing forced area by an equivalent percentage.According to Equation (6), it seems that the force coefficients depend on the plate's forced area, but actually, the influence of the plate's geometry shape also has to be considered.As regards the value, the maximum lift coefficient is obviously higher than the mean drag coefficient, and the difference between them also increases with the increasing aspect ratio.This indicates that incoming flow has a stronger influence on the plate in the transverse direction than in the streamwise direction.As shown in Figure 5b, both the flapping amplitude and Strouhal number increase monotonically with the increasing aspect ratio.Moreover, their increasing slope tends to decrease, which indicates that the effect of aspect ratio on the plate's flapping amplitude gets weaker as it increases.It is obvious that the two curves are not wholly similar in the increasing trends.Concerning the definition of the Strouhal number, St = f A f /U ∞ , the flapping amplitude A f and flapping frequency f are both important influencing factors.Regarding the value, the Strouhal number falls in the range from 0.16 to 0.27 in the present study.Despite the difference in the Reynolds number, the similar values of the Strouhal number were obtained in the simulations of flexible filaments.It is much more interesting that even for an active flapping rigid hydrofoil propelled in a tank, the higher thrust efficiency was found to fall in the range of 0.2 ∼ 0.4, which is close to the present study [48].Furthermore, most fish and birds have been shown to swim or fly near an optimal Stvalue at 0.3 due to the high propulsion efficiency [49].The similarity in Strouhal number value in animals' propulsion has been found to be strongly related to the formation of the optimal vortex [50].These comparisons indicate that wake instability existing in the passive flapping plate is similar to that of swimming fish or flying birds in nature.

Flow around an Isolated Flexible Plate
The hydrodynamics and flapping motion have close associations with the fluid characteristics.Even if the flow is moving, the relation of fluid-fluid and fluid-structure can be identified to some extent.Figure 6   Figure 7 shows the instantaneous vortical structures of an isolated flexible plate from the top and axonometric view, including (a) AR = 0.5 and (b) AR = 2.0.Both of the cases show a row of a hairpin-like vortical structure arrayed in order behind the plate.The vortical structures are symmetrical about the plate's centerline, corresponding to the plate's spanwise deformation.Due to the vortex shedding from the plate being the combination of the side edge vortex and trailing edge vortex, the difference in strength between them has a significant effect on the final combined shape.To observe this difference more clearly, the vortical structures have been labeled as V1 and V2 in the figure, respectively.For AR = 0.5, the downstream vortical structures are close to a double line shape.This indicates that the plate's trailing edge vortex has been weakened under the low aspect ratio.For AR = 2.0, the downstream vortical structures are close to the O shape.This indicates that both of the vortices shedding from the side edges and trailing edge have been strengthened under a high aspect ratio.

Coupling Motions and Vortical Structures
To understand the interaction dynamics between two flexible plates in a uniform flow, we perform the simulations for side-by-side arrangement.According to the results, the flapping behavior of each plate is asymmetrical about its local x b -axis when the transverse gap distance is small.Meanwhile, the flapping motion of the upper plate is similar to the lower one about the flow domain's center-plane in the Z-axis direction.For simplicity, the dynamics of the upper plate will be taken for the following description due to their symmetry or similarity.It also should be noted that the following relative flapping amplitude and force coefficient have been scaled by the result of an isolated plate under the same flow condition.
Figure 8 plots the relative flapping amplitude versus transverse distance, including (a) at the outer side and (b) at the inner side.As shown in Figure 8a, Z B,out of the plates under a small aspect ratio (AR < 1.0) increases monastically with the increasing transverse distance.For instance, as AR = 0.5, Z B,out increases smoothly and tends to reach the value 1.0.Obviously, Z B,out is less than 1.0 when the two plates are placed nearby, which means that the plates' flapping is weakened; here, the largest weakened extent is about 10%.However, as the transverse distance increases, Z B,out of the plates under a large aspect ratio (AR > 1.0) decreases correspondingly until reaching the lowest point, then starts to increase smoothly.For instance, as AR = 1.5, the value of Z B,out decreases from above the 1.0 level to below the 1.0 level and then starts to increase smoothly.This trend means that the plate's flapping status transits from strengthened to weakened, and the largest varied extent is about 24% and 5%, respectively.Obviously, the wider (higher aspect ratio) plates can achieve better flapping amplification than the longer plates.In other words, the longer plates in the streamwise direction can achieve better flapping weakened performance in the side-by-side configuration.
As shown in Figure 8b, the trends of Z B,in of the plates under a small aspect ratio are similar to the Z B,out .However, the trends of Z B,in of the plates under a large aspect ratio show some differences.Take AR = 1.5 for instance again; as the transverse distance increases, Z B,in of the plate increases until it reaches the peak.The value of Z B,in is less than 1.0 when the two plates are placed nearby.Recall that Z B,out is increasing, and above 1.0 in this regime, it can be inferred that the flapping motions towards the inside are resisted by the advection in the gap domain.The ratio of Z B,out to Z B,in is always larger than 1.0, and this value gets higher as the aspect ratio increases, meaning that the gap flow has a stronger influence on the wider plates.Despite the difference in various trends, the plates share the common point that the flapping asymmetry phenomenon tends to disappear as the gap distance increases.Based on a variety of simulations, three distinct coupling modes of the two side-by-side plates are identified depending on the transverse gap distance.Hereafter, the three coupling modes of the plates' motion are named as single-plate mode, symmetrical flapping mode and decoupled mode, respectively.To clearly introduce the dynamics and interaction mechanisms of these typical motion modes, the plates at AR = 1.5 are taken as examples.According to the computation results, the parameters at G z = 0.1, G z = 3.0 and G z = 4.0 are corresponding to single-plate mode, symmetrical flapping mode and decoupled mode, respectively.The dynamics of these cases are briefly analyzed as follows.
In the single-plate mode, the two plates are flapping in phase with just a row of vortex existing in the wake.Figure 9 shows the time evolutions of the transverse displacements of points B 1 and B 2 at G z = 0.1.As shown in the figure, the two plates flap synchronously with the same phase.Meanwhile, the two plates' flapping motion share the equivalent frequency and amplitude.Despite the two displacement curves not overlapping, the peak of Z B,up is equal to the absolute value of minimum Z B,un .Therefore, the motion of each plate is asymmetry about its local x b -axis, but the two plates' flapping amplitude is symmetrical about the flow domain's transverse center-plane.Figure 10 shows the envelops of the two flapping flexible plates at G z = 0.1, including (a) envelops in 3D and (b) envelops of the plates' mid-lines.In Figure 10a, the upper and lower plate are displayed in a semi-apparent mesh and a solid surface, respectively.As shown in Figure 10a, the deformation of the upper plate is always similar to the lower one.Additionally, the deformation in the spanwise direction still exists in this situation and keeps symmetrical about the flow domain's center-plane in the Y-axis direction.As shown in Figure 10b, each of the trailing edge's mid-point shows an "8" shape trajectory from the side view.The point N 1 denotes the intersection node of the trajectory of point B 1 , so as with point N 2 to point B 2 .It is observed that N 1 and N 2 are symmetrical about the transverse mid-plane; this result is in accord with Figure 9. Figure 11 shows the time evolutions of the vortical structures of two side-by-side flapping plates at G z = 0.1.The instants t 0 ∼ t 3 in the figure are equal to those in Figure 9.As the plates move towards the positive Z-axis direction, a clockwise vortex gradually forms on the upper surface of the upper plate.As the plates start to move towards the negative Z-axis direction, this fully-developed vortex V2 starts to shed from the trailing edge.Concerning the fluid-induced flapping, the strength of vorticity in both sides of the plates is stronger than in the gap flow; this causes the plate's asymmetrical flapping motion to some extent.Although there are two flexible plates in the flow domain, just one vortex sheds in half of the flapping cycle.The vorticity pattern in this situation is similar to an isolated plate in the same flow condition, so it is called single-plate mode.In the symmetrical flapping mode, the two plates are flapping out of phase, and two rows of vortices with opposite rotating directions exist in the wake.Figure 12 shows the time evolutions of the transverse displacements of points B 1 and B 2 at G z = 3.0.As shown in the figure, the two plates flap synchronously, but with the opposite phase.For instance, as the upper plate's displacement is moving upwards its peak, the lower plate's displacement is moving downwards its trough at the same time.Each of the plates' displacements are equal in absolute value, but with the opposite sign.Meanwhile, the two plates' flapping is periodical and sharing equivalent frequency and amplitude.Clearly, the motion of each plate is symmetrical about its local x b -axis, and the two plates' flapping is symmetrical about the flow domain's transverse center-plane.Figure 13 shows the envelops of the upper flapping plate at G z = 3.0, including (a) the axonometric view and (b) the side view.As the plate's flapping is similar to the other, here just the envelops of the upper one are shown in a semi-apparent mesh.It can be observed that the spanwise deformation still exists in this coupling flapping mode and keeps symmetrical about the flow domain's center-plane in the Y-axis direction.As shown in Figure 13b, the trailing edge's mid-point shows an "8" shape trajectory from the side view, as well.Generally, the deformation of the lower plate in this coupling mode is similar to that of an isolated one.Figure 14 shows the time evolutions of the vortical structures of two side-by-side flapping plates at G z = 3.0.The instants t 0 ∼ t 3 in the figure are equal to those in Figure 12.As the upper plate moves downwards (from t 0 to t 2 ), a counter-clockwise vortex (labeled V2) forms on the under surface and gradually separates from the trailing edge.During the same time, the lower plate moves upwards, and a clockwise vortex (labeled V1) forms on the upper surface and gradually separates.It is observed the fully-developed vortex shed from one plate will collide with that from the other one.This collision only occurs in the domain away from the trailing edges, so the influence to the flapping motion's symmetrical characteristic is weak.However, this collision between the opposite vortex will weaken the strength of the separation flow streamwise, then lead to the phenomenon of drag reduction.Considering that the flapping motion and vortex pattern of the plates are symmetrical, we name this coupling status the symmetrical flapping mode.In the decoupled flapping mode, the two plates are flapping in phase, and two rows of vortex with the same rotating direction exist in the wake.Figure 15 shows the time evolutions of the transverse displacements of points B 1 and B 2 at G z = 4.0.As shown in the figure, the two plates flap synchronously with the same phase.Meanwhile, the two plates' flapping motions share equivalent frequency and amplitude.Obviously, the two periodical displacement curves almost overlap.Regarding the value of displacement, the flapping motion of each plate is symmetrical about its local x b -axis.The envelops of the upper flapping plate at G z = 4.0 are also similar to those of an isolated one (see Figure 3).Figure 16 shows the time evolutions of the vortical structures of two side-by-side flapping plates at G z = 4.0.The instants t 0 ∼ t 3 in the figure are equal to those in Figure 15.As the upper plate moves downwards (from t 0 to t 2 ), a counter-clockwise vortex (labeled V2) forms on the under surface and gradually separates from the trailing edge.During the same time, the lower plate moves downwards, as well, and a counter-clockwise vortex (labeled V1) also forms on the under surface and gradually separates from its trailing edge.The two vortexes generated in one-half cycle keep independent without obvious collisions or other interactions.Generally, the flapping motion and vortex pattern of each plate are extremely similar to each other, so it seems like there does not exists any coupling behavior between them.Hence, this coupling status is called the decoupled mode.From the discussions above, it is obvious that the transverse gap distance can evidently change the flapping dynamics and coupling modes of the two side-by-side flexible plates.The interactions between the two plates are decreasing with the increasing gap distance.According to the characteristics of flapping motion and fluid pattern obtained in the present study, three typical wake modes in such plates' configuration are classified.The three distinct modes are categorized as single-plate mode, symmetrical flapping mode and decoupled mode, respectively.To show the features more clearly, Figure 17 plots the three typical vortical patterns in the wake of two side-by-side flexible plates.In this figure, the vortical structures are plotted by the λ 2 criterion, including the top view and side view in each subfigure.As shown in the top views (XY plane), all of the plates' wakes exist in a chain of hairpin-like vortical structures, which are arrayed in order along the streamwise direction.However, the wake patterns are obviously different in the side views (XZ plane).In the single plate mode, the vortex is extremely structured like the result of an isolated plate.In the symmetrical flapping mode and decoupled mode, the two rows of vortical structures are symmetrical and parallel, respectively.

Force and Energy Conversion
The plate's deformation is associated with the fluid force and energy harvesting performance.Figure 18 shows the relative force coefficient versus transverse distance, including (a) the maximum lift coefficient and (b) the mean drag coefficient.As shown in Figure 18a, the trends of lift amplitude are similar to those of the flapping amplitude.This can be explained briefly by the fact that the plates' passive flapping is mainly dependent on the lift fore on the boundary.The higher lift force can induce a larger flapping amplitude.As shown in Figure 18b, the trends of the drag coefficient are similar to those of the flapping amplitude, as well, but for different reasons.This can be explained briefly as follows.The drag force on the plate's boundary is due to the resistance to the incoming flow.This resistance force is proportional to the up-wind area.Hence, the higher amplitude produces a higher projection area vertical to the streamwise, leading to the production of a higher drag force.Regarding the value, the phenomenon of drag reduction is observed when the plates' aspect ratio is less than 1.0 or the higher aspect ratio plates flap under a large transverse distance.Here, we examine the relation between coupling modes and energy capture efficiency.Once the plate's flapping motion is induced by the incoming flow, the deformation is highly dependent on the bending energy captured from the fluid's kinetic energy.The bending energy E b (t) is expressed as: The conversion efficiency from the fluid kinetic energy to the plate's bending energy can be defined as [51]: where E b,max is the maximum bending energy at the largest deformation status; E k is the kinetic energy of the incoming fluid passing through the maximum up-wind area (Z B,out W) during a bending phase; T out is the time from when point B crosses the equilibrium position to when it reaches the peak position.The performance of energy harvesting is shown in Figure 19, including (a) the bending energy and (b) the efficiency of bending energy conversion.In Figure 19a, the red dashed line denotes the result of an isolated plate.It is shown that E b,max /E b0,max of the wider (AR = 1.5) plate decreases monotonically as G z increases, while the result of the slender (AR = 0.5) plate increases monotonically.As G z increases, both of the results tend to approach the value 1.0.Obviously, the wider plate in side-by-side configuration can capture more energy than in the isolated situation, while the slender plate captures less bending energy.For the purpose of harvesting more bending energy, the optimized selection is the higher aspect ratio plate flapping in the single-plate mode.
In Figure 19b, the η for both of the plates under different AR is decreasing monotonically with G z .Compared to the isolated result, it is found that η is amplified when the plates are placed in the side-by-side configuration.In the present study, the maximum η is 0.57, which is about 32% higher than in the isolated situation.For the purpose of improving the efficiency of bending energy conversion, the plate in the side-by-side configuration is better than in the isolated situation.More specifically, the optimized selection for higher bending energy conversion efficiency is the plate flapping in single plate mode, especially the plate with a higher aspect ratio.To predict the actual energy harvesting performance, the conversion of bending energy to other direct useable energy and its interaction with the plate's flapping dynamics should be studied, and this is beyond the scope of the present work.

Conclusions
The fluid-structure interaction problems of finite aspect ratio flexible plates in a uniform flow have been numerically studied.In the framework of the simulations, the lattice Boltzmann method combined with the immersed boundary and finite element method was used to solve the flow field and the plate's passive flapping.To compute the equations efficiently in the 3D flow field, the fluid domain was discretized by unifying Cartesian grids, while the structure domain was discretized by Lagrangian grids.To validate this coupled algorithm, the dynamics of flow past a square flexible plate were simulated at low Reynolds number.The results of the validation show good agreements with previous studies by CFD simulations, indicating the efficiency and accuracy of the present method.Then, numerical simulations based on various parameters were carried out, and the mechanisms governing the motions of flexible plates were analyzed systematically.The results related to the plates' passive flapping dynamics are briefly summarized as follows.
Firstly, the effects of aspect ratio on the fluid-structure interaction dynamics of an isolated plate in a uniform flow have been studied.We have found that the flexible plate can flap periodically with a large amplitude.This passive flapping behavior contains transverse and spanwise deformation.The flexible plate with a larger aspect ratio induced a larger flapping amplitude and Strouhal number.Furthermore, the flexible plate with a larger aspect ratio was subjected to greater drag and lift force.Further, the vortex shed from the side-edge was merged with the one shed from the trailing-edge and formed a hairpin-like structure.The shedding of the side-edge vortex was the main cause for its spanwise deformation.
Secondly, the effects of aspect ratio and transverse gap distance on the hydrodynamics interactions between two side-by-side flexible plates have been studied.With the transverse distance getting higher, we have found three distinct coupling modes of the plates' motion, including single-plate mode, symmetrical flapping mode and decoupled mode.In the single-plate mode, both of the two plates flapped synchronously, and there was a row of vortex shed from the plates' system, acting like the characteristic of a single plate in the same flow.In the symmetrical flapping mode, the two plates flapped out of phase, and there was a row of vortex shed form each plate, but with opposite directions.In the decoupled mode, each of the plates flapped synchronously, and there was a row of vortex with the same direction shed from each one, acting like two isolated plates flapping passively in the same flow field.Compared to the flow past an isolated plate, when the plates' motion was in the single-plate mode, those with a higher aspect ratio got larger flapping amplitudes.Further, the phenomenon of drag reduction was found when the plates' aspect ratios were less than 1.0 or the higher aspect ratio plates flap in the symmetrical flapping or decoupled mode.For the energy conversion performance, the wider plate in the side-by-side configuration can capture more bending energy than in the isolated situation, while the slender plate captures less bending energy than in the isolated situation.For the purpose of improving the bending energy conversion efficiency, the plate in the side-by-side configuration is better than in the isolated situation.More specifically, the optimized selection for higher bending energy conversion efficiency is the plate flapping in single-plate mode, especially the plate with a higher aspect ratio.
This work has provided several physical insights into the understanding of aquatic animal swimming and is expected to inspire the development of bionic thrusters or renewable energy harvesters.However, the fluid-structure interaction problem in biological propulsion is far more complex and diverse than the model considered here.Future work will focus on the experiments.Furthermore, the hydrodynamic performances of multiple flexible plates with heaving motion in a 3D flow remain to be elucidated.

Figure 1 .
Figure 1.Schematic of the flow over two flexible plates in side-by-side arrangement: (a) axonometric view; (b) side view.

Figure 2 .
Figure 2. Comparison of the present result with previous data reported by Huang and Sung [31].

Figure 3 .
Figure 3. Envelop of an isolated flapping flexible plate: (a) 3D envelop of the flapping plate; (b) envelop of the plate's mid-line.

Figure 4 .
Figure 4. Vortical structure of the flow over an isolated flexible plate (Supplementary Video S1 is available online).

Figure 5
Figure 5 shows the force and displacement versus aspect ratio AR, including (a) mean drag coefficient C D and maximum lift coefficient C L,max and (b) flapping amplitude A f and Strouhal number St.In this study, the aspect ratio is changed from 0.5 to 2.0.As shown in Figure5a, both the mean drag coefficient and the maximum lift coefficient increase almost linearly with the increasing aspect ratio.As regards the increasing slope, the maximum lift coefficient is clearly higher than that of the mean drag coefficient.This indicates that the maximum lift coefficient is much more sensitive to the aspect ratio.Despite both of the force coefficients monotonically increasing with the aspect ratio, the slopes are less than 1.0.This means that the forces on the plate are not able to increase with the increasing forced area by an equivalent percentage.According to Equation (6), it seems that the force coefficients depend on the plate's forced area, but actually, the influence of the plate's geometry shape also has to be considered.As regards the value, the maximum lift coefficient is obviously higher than the mean drag coefficient, and the difference between them also increases with the increasing aspect ratio.This indicates that incoming flow has a stronger influence on the plate in the transverse direction than in the streamwise direction.

Figure 5 .
Figure 5. Force and displacement versus aspect ratio AR: (a) mean drag coefficient C D and maximum lift coefficient C L,max ; (b) flapping amplitude A f and Strouhal number St.
shows the snapshot of tracing particles from the top and side view, including (a) AR = 0.5 and (b) AR = 2.0.The particles are released from the sources at the upstream of the plate with a time interval of 10.It should be noted that the particles' motion is following the local fluid velocity, but does not affect the fluid motion.As shown in the figure, each pathline of the tracing particle shows a spiral trajectory in the streamwise direction, which is formed in the process of the plate's flapping.Due to the flexible plate deforming at most instants, the free flapping part of the structure always resists the incoming flow to move downstream in the streamwise direction.Hence, the particle driven by the flow is inclined to separate from the side edge when it encounters the resistance of the trailing edge.This separation behavior consequently induces the plate's spanwise deformation to some extent.Therefore, the symmetry motion of the tracing particles indicates the symmetry deformation of the plate spanwise.Compared to the results of AR = 0.5, those of AR = 2.0 are much more complicated, especially in the knot part of the two pathlines with the same coordinate in the Y-axis direction, as labeled by k 1 and k 2 in the figure.

Figure 6 .
Figure 6.Snapshots of the tracing particles from the top and side view: (a) AR = 0.5; (b) AR = 2.0 (Supplementary Videos S2 and S3 are available online).

Figure 9 .
Figure 9.Time evolutions of the transverse displacements of points B 1 and B 2 at G z = 0.1.

Figure 10 .
Figure 10.Envelops of the two flapping flexible plates at G z = 0.1: (a) envelops in 3D; (b) envelops of the plates' mid-lines.

Figure 11 .
Figure 11.Time evolutions of the vortical structures of two side-by-side flapping plates at G z = 0.1.

Figure 12 .
Figure 12.Time evolutions of the transverse displacements of points B 1 and B 2 at G z = 3.0.

Figure 14 .
Figure 14.Time evolutions of the vortical structures of two side-by-side flapping plates at G z = 3.0.

Figure 15 .
Figure 15.Time evolutions of the transverse displacements of points B 1 and B 2 at G z = 4.0.

Figure 16 .
Figure 16.Time evolutions of the vortical structures of two side-by-side flapping plates at G z = 4.0.

Figure 18 .
Figure 18.The relative force coefficient versus transverse distance: (a) maximum lift coefficient; (b) mean drag coefficient.

Figure 19 .
Figure 19.Performance of energy harvesting versus transverse distance: (a) bending energy; (b) efficiency of bending energy conversion.