Numerical investigation on the propulsive performance of flexible flapping fins using CFD/CSD method

2.1. Computational dynamics

Considering the unsteady flow structure and low-Reynolds-number characteristics of flexible flapping fins, arbitrary Lagrange Euler and preconditioned method was developed. The dual time-stepping three-dimensional preconditioned Navier-Stokes equations can be written in ALE form as:

where$\tau$and$t$are the pseudotime and physical time, respectively $\mathbf{W}=(\rho \rho u\rho v\rho w\rho e{)}^T$, $\mathbf{Q}=(puvwT{)}^T$. Preconditioning matrix is:

where $W={\left(\rho ,\rho u,\rho v,\rho w,\rho e\right)}^T$, $F\left(W\right)$ is the convective fluxes, $Fv$ is the viscous fluxes, $v_g=\dot{x}\cdot \overrightarrow{n}$ is the motion of the control volume $V$, $\dot{x}$ and $\overrightarrow{n}$ denote, respectively, the time-varying velocity and the normal vector of the control volume contour. In Eq. (1), $Q=(P,u,v,w,T{)}^T$, $\mathrm{\Gamma }$ is the preconditioning matrix which is employed to solve the system stiff problem, since the solver we primitively developed was to solve compressible flow. Here $\tau$ and $t$ designate the pseudo and physical time, respectively.

In the present study, $v$ is arbitrarily defined, corresponding to the so-called Arbitrary Lagrangian Eulerian (ALE) formulation. The KW-sst turbulence model is employed to close the N-S equations.

A detailed description of the in-house programmed URANS solver can be found in [13, 14].

2.2. Co-rotational finite element analysis

The flexible flapping fins exhibit highly nonlinear dynamic behavior. Total Lagrangian (TL), Update Lagrangian (UL) and Co-rotational (CR) formulation are three main methods for geometric nonlinear finite element analysis [15]. Currently, the co-rotational formulation has been extensively applied to geometric nonlinear analysis due to its efficiency and universality. The general form of in-house dynamic finite element equations is:

where $\mathbf{q}$ is the nodal DOF vector, $\mathbf{v}$ is the velocity vector and $\mathbf{a}$ the acceleration vector. $\mathbf{M}$, $\mathbf{C}$, $\mathbf{R}$ are mass stiffness matrix, damping matrix and tangent stiffness matrix, respectively. A fin-deform based three-node triangular shell element with 18 degree of freedom was developed by combining the optimal membrane element (OPT) and discrete Kirchhoff triangle (DKT) plate bending element, considering the geometric nonlinear analysis of thin flexible flapping fins with large rotation and small strain. Newton-Raphson iteration algorithms combined with an automatic load controlled technology were used for solving nonlinear equations. The detailed description of the co-rotational nonlinear finite element method can be found in [15-17].

3.1. Coupling strategy

A staggered second-order accurate coupling algorithm by Farhat [18] was used which allows exchanges of motion (the kinematics mapping from the structure solver to the fluid solver) and loads (from the fluid to the structure solver) at the fluid-structure interface. Fig. 1 illustrates the loose-coupling process. In view of the fact that the domains and meshes of the two solvers are non-conformal, the mapping of these incompatible domains and meshes was based on RBF (radial based function) method that interpolates data between the structural and the aerodynamic domains.

Fig. 1Coupling process and data exchange

3.2. Coupling strategy

In view of the inherent large flapping motions, an efficient and robust dynamic mesh strategy was required for the flapping motion and the mesh motion caused by the structure deformation. In this paper, a RBF-based deformable overset grid proposed by Deng et al [19] was employed and is roughly addressed here as:

1. Overset mesh strategy is used for large flapping motion

2. RBF method is used for local structure deformation

3.3. Outline of Fluid-structure interaction

To fully understand the process, a step-by-step solution procedure to solve the system of FSI equations is given as follows:

1. Meshes and calculation parameters are input and initialized at the beginning.

2. A new timestep is set up based on flow and structure characteristics. A detail selection of timestep is introduced in chapter 5.2.

3. Rigid motion is defined at current timestep and the position and orientation of the global/flapping frames are calculated. The previous timestep variables such as displacements and forces are then updated.

4. Aerodynamic and structure computations are carried out at current timestep using in-house fluid-structure interaction solver. Usually several iterations are needed in one timestep.

5. Displacements and forces are transferred at fluid/structure interfaces using RBF and coupling methods after each iteration.

6. After all the iterations in one timestep, check the convergence for the current timestep by evaluating the forces transferred form fluid solver and calculated in the structure solver using convergence tolerance. If the solution is not converged, return to the beginning of current step and choose a revised tiemstep; if it does, proceed to the next time step.

Key steps of the in-house CFD/CSD method is highlighted in Fig. 2.

Fig. 2Outline of fluid-structure interaction

4.1. 2-D chordwise flexible foil

The unsteady flow structure around a 2-D chordwise flexible foil was simulated in the in-house CFD/CSD solver. A benchmark validation was carried out based on the result from [8], with main parameters listed in Table 1. The computational overset mesh system (see Fig. 3) was created using an unstructured mesh with a refinement at the boundary layer region of $y_{plus}=$1. The mesh was moved entirely during the computation following the specific pitching motion.

Fig. 3a) Overset grid system and b) flexible foil at certain moment

a)

b)

A comparison of vortex distribution between Marais’s experiment and present FSI solver is presented in Fig. 4(a) and (b). As can be seen, the time varying vortex distributions are in reasonably good agreement with the experiment results. The slight differences of vortex distances suggest that the unsteady flow structure around flapping wings and wing stiffness distribution and deformation are difficult to precisely simulate. Overall, it can be seen that the computed result shows good agreements.

Fig. 4Comparison of vortex distribution between a) Marais’s experiment [8] and b) present solver

a)

b)

4.2. 3-D spanwise flapping wing

A spanwise 3-D flexible flapping-wing was simulated to further demonstrate the accuracy of the CFD/CSD solver. Experiments were conducted by Heathcote from [7] in a water tunnel on a straight, untapered wing, in pure plunging motion at Reynolds numbers ranging from 10,000 to 30,000. The main parameters are listed in Table 2.

Fig. 5 shows the overset grid system of the computational domain around the wing, with structured meshes about 6 million cells and a refinement at the boundary layer region. The overset mesh grid system moved according to the specified motion law and result comparisons are illustrated in Fig. 6(a) and 6(b).

Fig. 5a) Overset grid system for 3-D wing model and b) CSD meshes

a)

b)

Fig. 6a) Tip distance and b) thrust coefficient comparison

a)

b)

As shown in Fig. 6, both the normalized vertical displacements and the thrust coefficients are in good agreement with the measurements. Both phase and amplitude agreements demonstrated that the present FSI-based model is capable to reasonably predict the flexible wing.

It is seen that the above computed results showed good agreements with experiment data from both cases. However, slight difference in both lift and thrust coefficients especially the peak values suggested that wing stiffness distribution and deformation are difficult to extremely precisely simulate despite the complex flow structure. In general, agreements in both phase and amplitude demonstrate that the present FSI-based model is capable of reasonably predicting the flexion of flexible wings.

In the meanwhile, three main aspects can be considered for future investigations on the accuracy:

(1) Fluid turbulence models can be modified specifically for flapping motions to improve precision.

(2) Co-rotational finite elements can be further improved considering the large geometric nonlinear characteristics and dynamics.

(3) High accuracy algorithms can be developed at the interface of fluid-structure interaction instead of the current loose coupling.

5.1. Computational setup

So far we have addressed that the present FSI model is robust and efficient in modeling aeroelastic structural dynamics and aerodynamics of flexible fins. In this section, an extensive discussion was discussed on how the inherent flexibility of flapping fins affects its aerodynamic performance in terms of the flow structure, the force-production and efficiency. A high flexible fin with NACA0015 as its section was presented in the study to explore the effect of flexion on the propulsive performance and aerodynamic efficiency of fins, with schematic in Fig. 7 and main parameters in Table 3.

Fig. 7Schematic of high flexible fin (unit/mm, half model)

5.2. Motion law illustration

The kinematic model of flapping fins was constructed based on natural flyers and former experiments. Two main basic motions were considered to describe the flapping fin movement: (1) flapping about the $X$ axis, described by the fin plunging angle, and (2) rotating about the $Y$ axis, described by the fin pitching angle. The equations of motion were specified here in Eq. (3-4):

In detail, the preconditioned motion of flapping fins that contains eight main phases is illustrated in Fig. 8 with main features introduced as follows.

1. Flapping and twisting. Fins rotate about their spanwise and chordwise axis, generating flapping and twisting motions. Downstroke and upstroke are the basic motions the fin makes whilst flapping. Twisting is the primary method of controlling the fin’s angle of attack. It is also how the fin pronates and supinates as it goes through the full flap cycle.

2. Downstroke and upstroke. The downstroke consists of the fin moving from its highest positive fin flap angle, through to the most negative fin angle ($t=$1/4$T$-3/4$T$). The upstroke is the opposite of this motion ($t=$3/4$T$-1/4$T$). The combination of a full downstroke and upstroke makes a flapping cycle.

3. Pronation and supination. In natural flyers, the leading edges are always strong. Therefore, leading edges are always some phase angles ahead of rear parts. At the end of a downstroke, leading edges are rotated through a large angle such that the incidence angle remains constant through the subsequent upstroke. The top surface of the fin then becomes the bottom surface. This process is the supination of the fin. At the end of the upstroke, the fin rotates through a similar angle in the opposite direction, so that the leading edge is still leading and the surfaces have again switched; this is the pronation of the fin.

Fig. 8Illustration of motion of high flexible fin (half model)

To investigate and fully understand the effect of flexion on the propulsive performance and aerodynamic efficiency of flexible flapping fins, computations were conducted by measuring low/medium/high flexible models. An extensive investigation on the effect of flexion was then conducted. Different flapping frequencies and fin flexibilities were studied to find out the effect of these parameters on propulsive performance and aerodynamic efficiency.

Usually, a proper grid is chosen based on trade-offs between computational accuracy and computational costs. Considering the fact that meshes in CFD analysis are usually more refined than those in CSD meshes. In this study, a fine grid (30×24×10 elements), medium grid (60×48×10 elements) and coarse grid (15×12×5 elements) were initially chosen for CSD empirically. Also, a fine grid (60×48×10 elements), medium grid (120×96×20 elements) and coarse grid (30×24×5 elements) were initially chosen for CFD, respectively. The simulation results of these grids show that the result of coarse grid is less accurate than the others and the medium grid costs less than that of fine grid. Thus, medium grid, which is accurate enough, was chosen.

Fig. 9 shows the CFD (60×48×10 elements) and CSD (30×24×5 elements) meshes used in this paper.

Fig. 9a) CFD meshes and b) CSD meshes

a)

b)

In CFD analysis, only half model was used considering the symmetrical characteristics. Velocity inlet and pressure outlet boundary were adopted to simulate the farfield boundaries. In CSD analysis, all degrees of freedom were set stationary for inner part of the fin, while the outer part of the fin moved according to the motion law.

Different parameters and a time step of 0.001 were used for all the simulations.

First, flexion was studied in cases with different ${\theta }_{twist}^{}$ and ${\theta }_{flap}^{}$. In order to investigate the effect of flexion on thrust, computation is operated on ${\theta }_{twist}^{}=$ 10°, 15°, 20°, 25°, 30° and ${\theta }_{flap}^{}=$ 10°, 15°, 20°, 25°, 30°. Plots of thrust coefficient against ${\theta }_{twist}^{}$ and ${\theta }_{flap}^{}$ is show in Fig. 10 for different flexibilities. It is shown that $C_T$ increases with the increment of flapping angle and twisting angle. Among the two angles, the thrust increases faster with the same increment of twisting angle than that of flapping angle, which means that the twisting angle plays a more important role on thrust performance.

Fig. 10Thrust coefficient for different plunging and pitching angles

Fig. 11Efficiency for different flexibility and Strouhal number

Fig. 12Pressure coefficient distribution in one cycle (plunging angle = 25°, pitching angle = 25°, Strouhal number = 0.5)

a)$t=$0.125$T$

b)$t=$0.25$T$

c)$t=$0.375$T$

d)$t=$0.5$T$

e)$t=$0.625$T$

f)$t=$0.75$T$

g)$t=$0.875$T$

h)$t=$0/1$T$

A further study of the relationship between frequency, flexion and efficiency was carried out. Aerodynamic efficiency is plotted as a function of flapping frequencies. In this paper, $f=$0.05 Hz, 0.1 Hz, 0.15 Hz, 0.2 Hz, 0.3 Hz, 0.4 Hz and 0.8 Hz. It is seen that the aerodynamic efficiency firstly increases sharply with increasing flapping frequency and displays a decreasing trend afterwards. Peak efficiency is a matter of appropriate combination of frequency and flexion, which lies in $f=$0.2-0.4. Also, the efficiency increases and then decreases when Strouhal number goes up within the same flexible model and peaks within the range of 0.2 $<St<$0.4, which lies in the middle of the range observed in nature. The efficiency of the intermediate flexible fin is significantly higher than that with low flexibility and high flexibility. Therefore, it can be concluded that the flexibility can offer significant efficiency benefits, which agrees with other studies referred to in the introduction.

Distribution of pressure coefficient in one cycle with 25-degree flapping angle and twist angle is presented in Fig. 12. It is found from the figure that upper surfaces of the highly flexible wings present large lift loss, which partially explains the reason of thrust decreases.

Total deformation distribution at the end of one cycle is presented in Fig. 13. As can be seen, structures of fins deform under the pressure computed in CFD analysis. The maximum deformation is shown at the tip of fins, which is in accord with the structure stiffness distribution. Also, the maximum deformation suggests that medium flexion can maintain its structure under pressure and keep good aerodynamic and structural performances.

Fig. 13Total deformation distribution at the end of one cycle (plunging angle = 25°, pitching angle = 25°, Strouhal number = 0.5)