2.1 Fluid analysis

The fluid analysis is based on a 3D preconditioned Navier–Stokes equation in order to simulate the low-Mach-number flows around the plunging wings. By introducing a preconditioned pseudo-time derivative term, an integral form of the non-dimensional governing equations under a free-stream condition is expressed as follows:

(1) $${{\rm d} \over {{\rm d}t}}\mathop{\int}_V \underline{W} {\rm d}V{\plus}\underline{{\underline{\Gamma } }} _{a} {{\rm d} \over {{\rm d\tau }}}\mathop{\int}_V \underline{Q} _{p} {\rm d}V{\plus}\mathop{\int}_S \overleftarrow{{\underline{F} }} \cdot \underline{{\hat{n}}} {\rm d}S{\equals}\mathop{\int}_S \overleftarrow{{\underline{F} _{v} }} \cdot \underline{{\hat{n}}} {\rm d}S,$$

where $$\underline{W} $$ is the conservative solution vector, and $$\underline{Q} _{p} $$ , $$\overleftarrow{{\underline{F} }} $$ and $$\overleftarrow{{\underline{F} _{v} }} $$ , are the primitive solution vector, inviscid flux vector and viscous flux vector, respectively. Each vectors are defined by the following relation:

(2) $$\underline{W} {\equals}\left\{ {\matrix{ {\rm \rho } \cr {{\rm \rho }v_{x} } \cr {{\rm \rho }v_{y} } \cr {{\rm \rho }v_{z} } \cr {{\rm \rho }E} \cr } } \right\},\quad \overleftarrow{{\underline{F} }} {\equals}\left\{ {\matrix{ {{\rm \rho }v} \cr {{\rm \rho }vv_{x} {\plus}p\hat{i}} \cr {{\rm \rho }vv_{y} {\plus}p\hat{j}} \cr {{\rm \rho }vv_{z} {\plus}p\hat{k}} \cr {{\rm \rho }E} \cr } } \right\},\quad \overleftarrow{{\underline{F} _{v} }} {\equals}\left\{ {\matrix{ 0 \cr {{\rm \tau }_{{xi}} } \cr {{\rm \tau }_{{yi}} } \cr {{\rm \tau }_{{zi}} } \cr {{\rm \tau }_{{ij}} v_{j} {\plus}q} \cr } } \right\}$$

$$\overleftarrow{{\underline{F} _{v} }} $$ is discretised using Roe’s approximated Riemann solver ⁽ Reference Roe ^{31 )} with MUSCL extrapolation ⁽ Reference van Lear ^{32 )}. High spatial accuracy is realised by using Van Albada’s limiter. By applying the gradient theorem over an auxiliary cell, the spatial derivatives of the solution vectors are obtained at the cell interfaces. And the derivatives are employed in order to compute the viscous flux vector, which is equivalent to the second-order central difference method on a regular grid.

In Equation (1), $$\underline{{\underline{\Gamma } }} _{a} $$ is the preconditioning matrix proposed by Weiss and Smith⁽ Reference Weiss and Smith ^{33 )}. By using the preconditioning matrix, it is possible to reduce the condition number of the system for low Mach number flows and to enhance the convergence of the iterative solution. And t and τ denote physical time and pseudo-time used in the time-marching procedure, respectively. For the unsteady flow analysis, the semi-discretised governing equations in fictitious time are analysed using a dual-time stepping method in conjunction with the approximate factorisation-alternate direction implicit method. The real-time derivative term is discretised by the second-order backward difference formula. After the convergence of the dual-time sub-iterations, the fictitious time term becomes zero, which ensures the second-order accurate solutions in time.

Details of the flow analysis are included in the literature⁽ Reference Yoo and Lee ^{34 )}. To consider structural deformation, a radial basis function is employed, and a geometric conservation law is utilised when computing the volume of the deforming grid.

2.2 Coupling strategy

A description of the current FSI framework is described in Refs 27 and 28. An implicit coupling approach is used to couple the structural analysis with the fluid dynamics analysis. Thus, updated coupled solutions are obtained for the same time-step at the end of the sub-iteration routine.

The current improvement of the FSI framework, compared with the authors’ previous works, is as follows. In the present fluid analysis, distributed loads are obtained by exerting pressures over the surface grids. A 3D shell analysis is used for the structural dynamics to develop improved data exchange between the fluid and structural boundaries. In the present numerical prediction, a wing stiffened by a metallic plate is employed. Thus, the following procedure for the interface is considered. First, the pressures acting on the upper and lower wing surfaces are collected. Then, those are transformed to the force vectors at grid points by integrating them with each area of the grid. And those are redefined to the force vector with respect to the intermediate plate at the centre of the wing. For the shell analysis, the relevant distributed loads are interpolated and transferred in the form of the nodal forces. A radial basis function is used for interpolation between the unmatched meshes in the fluid and structure. The present FSI framework is described in Algorithm 1.

3.1 CR framework

The CR formulation can be established by tracking a co-ordinate-based elemental motion.

The CR co-ordinate is an additional intermediate configuration between the undeformed and deformed frames. The local system can be defined between the CR and local co-ordinates. Figure 1 shows the co-ordinates and related transformations obeying the elemental kinematics for a three-node shell element. Beginning with the elemental fixed frame $$\underline{{\underline{E} }} _{I} $$ , the rotational matrix $$\underline{{\underline{R} }} _{o} $$ can be defined by tracking the elemental initial state. The rotational matrix $$\underline{{\underline{R} }} _{G} $$ can be defined by elemental rotational displacement referring to an undeformed configuration $$\underline{{\underline{E} }} _{G} $$ . The complete rotation can be decomposed into rigid body rotation referring to the CR frame $$\underline{{\underline{E} }} _{C} $$ and elastic deformational rotation referring to a deformed configuration $$\underline{{\underline{E} }} _{L} $$ . The origin of each co-ordinate is taken at the centroid. Each transformation can be accomplished by rotational operators $$\underline{{\underline{R} }} _{G} $$ , $$\underline{{\underline{R} }} _{r} $$ and $$\underline{{\underline{R} }} _{L} $$ , which are composed of the global, rigid-body and local rotations θ _G , θ _r and θ _L , respectively. Here, the local rotational variables are obtained by multiplying the operators: $$\underline{{\underline{R} }} _{L} {\equals}\underline{{\underline{R} }} _{r}^{T} \underline{{\underline{R} }} _{G} \underline{{\underline{R} }} _{o} $$ .

Figure 1 Co-ordinates and elemental kinematics in the CR formulation for 3-node shell.

It is noted that the origin of the co-ordinate is different from that used in the previous FSI study⁽ Reference Chimakurthi ^{14 )} and the present co-ordinates follows Battini’s approach⁽ Reference Battini ^{40 )}. Thus, it becomes possible to obtain stable solution without the effect of sequence in elemental nodes. Details of the advantages regarding the present approach are described in Battini⁽ Reference Battini ^{40 )}.

In the present formulation, the rotational operator can be defined by Rodrigues’ formula:

(3) $$\underline{{\underline{R} }} {\equals}\underline{{\underline{I} }} {\plus}{{\sin {\rm \rtheta }} \over {\rm \rtheta }}\underline{{{\tilde {\rtheta }}}} {\plus}{{1{\minus}\cos {\rm \rtheta }} \over {{\rm \rtheta }^{2} }}\underline{{{\tilde{\rtheta }}}} ^{2} $$

where a tilde denotes the skew-symmetric matrix. The spatial variation can be described using the following relation:

(4) $${\rm \delta }\underline{{\tilde{\phi }}} {\equals}\underline{{\underline{T} }} _{s} (\underline{{\rm \rtheta }} ){\rm \delta }\underline{{{\tilde{\psi }}}} $$

where $${\rm \delta }\underline{{{\tilde{\psi }}}} $$ and $${\rm \delta }\underline{{\tilde{\phi }}} $$ denote material and spatial angular variations (i.e. infinitesimal rotations superimposed onto the rotational matrix), respectively. The relevant variations obey the relation $${\rm \delta }\underline{{\underline{R} }} {\equals}({\rm \delta }\underline{{\tilde{\phi }}} )\underline{{\underline{R} }} {\equals}\underline{{\underline{R} }} ({\rm \delta }\underline{{{\rm \tilde{\psi }}}} )$$ . Further, $$\underline{{\underline{T} }} _{s} $$ is expressed as

(5) $$\underline{{\underline{T} }} _{s} (\underline{{\rm \rtheta }} ){\equals}\underline{{\underline{I} }} {\plus}{{1{\minus}\cos {\rm \rtheta }} \over {{\rm \rtheta }^{2} }}\underline{{{\tilde{\rtheta }}}} {\plus}{{{\rm \rtheta }{\minus}\sin {\rm \rtheta }} \over {{\rm \rtheta }^{3} }}\underline{{{\tilde{\rtheta }}}} ^{2} .$$

The local displacement $$\underline{u} _{L}^{i} $$ of the element is given by the following equation:

(6) $$\underline{u} _{L}^{i} {\equals}\underline{{\underline{R} }} _{r}^{T} \left\{ {\underline{r} ^{i} {\plus}\underline{u} _{G}^{i} {\minus}\underline{r} ^{c} {\minus}\underline{u} _{G}^{c} } \right\}{\minus}\underline{{\underline{R} }} _{o} \left\{ {\underline{r} ^{i} {\minus}\underline{r} ^{c} } \right\}$$

where the superscript i denotes an arbitrary node of the element. The geometric and rigid body rotations consisting of $$\underline{{\underline{R} }} _{o} $$ and $$\underline{{\underline{R} }} _{r} $$ are defined with respect to the elemental reference frame.

Let $$\underline{e} _{C} $$ denote the orthonormal basis vectors on the element. Then, the orthogonal matrix defining the rigid rotations is simply given by

(7) $$\matrix{ {\underline{{\underline{R} }} _{r} {\equals}[\underline{e} _{{C,1}} } & {\underline{e} _{{C,2}} } & {\underline{e} _{{C,3}} ].} \cr } $$

The first co-ordinate axis referring to the CR frame $$\underline{e} _{{C,1}} $$ is selected as the basis of the vector parallel to the side connected between nodes 1 and 2:

(8) $$\underline{e} _{{C,1}} {\equals}{{\underline{r} _{G}^{2} {\plus}\underline{u} _{G}^{2} {\minus}\underline{r} _{G}^{1} {\minus}\underline{u} _{G}^{1} } \over {\parallel\underline{r} _{G}^{2} {\plus}\underline{u} _{G}^{2} {\minus}\underline{r} _{G}^{1} {\minus}\underline{u} _{G}^{1} \parallel}}$$

Using an auxiliary vector, $$\underline{r} _{{3c}} {\equals}\underline{r} _{G}^{3} {\plus}\underline{u} _{G}^{3} {\minus}\underline{r} _{G}^{c} {\minus}\underline{u} _{G}^{c} $$ , the third co-ordinate axis, $$\underline{e} _{{C,3}} $$ , is assumed to be orthogonal to the elemental plane:

(9) $$\underline{e} _{{C,3}} {\equals}{{\underline{e} _{{C,1}} {\times}\underline{r} _{{3c}} } \over {\parallel\underline{e} _{{C,1}} {\times}\underline{r} _{{3c}} \parallel}}$$

The second co-ordinate axis, $$\underline{e} _{{C,2}} $$ , is then obtained by the following expression: $$\underline{e} _{{C,3}} {\times}\underline{e} _{{C,1}} $$ .

Thus, it is possible to establish the rotational operators, $$\underline{{\underline{R} }} _{o} $$ and $$\underline{{\underline{R} }} _{r} $$ . By obeying the elemental kinematics, the stiffness matrix and internal load vector for the triangular shells can be derived. The virtual energy can be defined by equating the local/global internal load vectors and the variation of the corresponding nodal displacement vectors:

(10) $$V{\equals}{\rm \delta }\underline{q} _{L}^{T} \underline{f} _{L} {\equals}{\rm \delta }\underline{q} _{G}^{T} \underline{f} _{G} $$

By equating the relationship between the local and global internal load vectors, Equation (10), the local internal load vector is defined as follows:

(11) $$\underline{f} _{L} {\equals}\underline{{\underline{B} }} \underline{f} _{G} $$

The global tangent stiffness matrix is defined by taking the variations of Equations (11):

(12) $$\underline{{\underline{K} }} _{G} {\equals}\underline{{\underline{B} }} ^{T} \underline{{\underline{K} }} _{L} \underline{{\underline{B} }} {\plus}{{\partial (\underline{{\underline{B} }} ^{T} \underline{f} _{L} )} \over {\partial \underline{q} _{G} }}{\rfloor}_{{\underline{f} _{L} }} $$

For the local system, $$\underline{{\underline{K} }} _{L} $$ , the OPT-DKT facet shell element⁽ Reference Khosravi, Ganesan and Sedaghati ^{42 )} is used in the current shell analysis. The transformation matrix $$\underline{{\underline{B} }} $$ can be defined by introducing the rigid body and deformational components:

(13) $$\underline{{\underline{B} }} {\equals}\underline{{\underline{E} }} \underline{{\underline{P} }} ^{T} $$

where $$\underline{{\underline{E} }} $$ is determined using the rotational operator $$\underline{{\underline{R} }} _{r} $$ , and $$\underline{{\underline{P} }} $$ is the projector matrix, considering the deformational components. Details of the derivations related to the stiffness matrix and internal load vector are presented in Refs 39 and 43.

3.2 Inertial terms based on CR framework

The derivation of the inertial force vector and elemental inertial matrices, i.e. the mass and gyroscopic matrices, is presented. The relevant formulation is based on Ref. 41. The kinetic energy of the element is expressed by considering the velocity vector:

(14) $${\cal K}{\equals}{1 \over 2}\mathop{\int}_\Omega {\rm \rho }\underline{{\dot{q}}} ^{T} \underline{{\dot{q}}} {\rm d\Omega }$$

where $$\underline{{\dot{q}}} $$ denotes the velocity vector. The kinetic energy can be re-expressed using the translational and rotational components:

(15) $${\cal K}{\equals}{1 \over 2}\mathop{\int}_\Omega \underline{{\dot{u}}} ^{T} m_{\rho } \underline{{\dot{u}}} {\plus}\underline{{{\dot{\rtheta }}}} ^{T} \underline{{\underline{I} }} _{{\rm \rho }} \underline{{{\dot{\rtheta }}}} {\rm d\Omega }$$

where m _ρ =ρ t, and $$\underline{{\underline{I} }} _{{\rm \rho }} $$ is the spatial dyadic tensor given by $${\rm diag}\left( {{{{\rm \rho }t^{3} } \over {12}},{{{\rm \rho }t^{3} } \over {12}},0} \right)$$ .

The kinetic energy variation can be expressed in terms of the material velocity and acceleration:

(16) $${\rm \delta }{\cal K}{\equals}\mathop{\int}_\Omega {\rm \delta }\underline{u} ^{T} m_{\rho } \underline{{\ddot{u}}} {\plus}{\rm \delta }\underline{{\rm \rtheta }} [\underline{{\underline{I} }} _{{\rm \rho }} \underline{{{\ddot{\rtheta }}}} {\plus}\underline{{{\tilde{\dot{\rtheta }}}}} \underline{{\underline{I} }} _{{\rm \rho }} \underline{{{\dot{\rtheta }}}} ]{\rm d\Omega }.$$

The material velocity and acceleration with respect to the local co-ordinate can be defined using elemental kinematics:

(17a) $${\rm \delta }\underline{{\dot{q}}} _{L} {\equals}\underline{{\underline{E} }} ^{T} {\rm \delta }\underline{{\dot{q}}} _{G} $$

(17b) $${\rm \delta }\underline{{\ddot{q}}} _{L} {\equals}\underline{{\underline{E} }} ^{T} {\rm \delta }\underline{{\ddot{q}}} _{G} .$$

Consequently, by employing matrix $$\underline{{\underline{N} }} _{s} $$ , which is composed of element nodal shape functions, the nodal quantities can be re-expressed in terms of those functions at the element level:

(18) $$\underline{q} _{L} {\equals}\underline{{\underline{N} }} _{s} \underline{q} _{L} {\equals}\left\{ {\underline{u} ^{T} \quad \underline{{\rm \rtheta }} ^{T} } \right\}^{T} .$$

The inertial load vector is derived from the relation $${\rm \delta }{\cal K}{\equals}\underline{f} _{{k,G}}^{T} \underline{q} _{G} $$ . Recalling the relation in Equation (11), the expression of the inertial load vector is obtained as follows:

(19a) $$\underline{f} _{{k,G}} {\equals}\left\{ {\underline{f} _{{k,G}}^{1} ^{T} \quad \underline{f} _{{k,G}}^{2} ^{T} \quad \underline{f} _{{k,G}}^{3} ^{T} } \right\}^{T} $$

(19b) $$\underline{f} _{{k,G}}^{i} {\equals}\underline{{\underline{E} }} \underline{f} _{{k,L}}^{i} {\equals}\underline{{\underline{E} }} \left\{ {\matrix{ {\mathop{\int}_{\Omega _{e} } N^{i} m_{{\rm \rho }} \underline{{\ddot{u}}} {\rm d\Omega }_{e} } & {} \cr {\mathop{\int}_{\Omega _{e} } N^{i} (\underline{{\underline{I} }} _{{\rm \rho }} \underline{{{\ddot{\rtheta }}}} {\plus}\underline{{{\tilde{\dot{\rtheta }}}}} \underline{{\underline{I} }} _{{\rm \rho }} \underline{{\rm \rtheta }} ){\rm d\Omega }_{e} } & {} \cr } } \right\},\quad i{\equals}1,2,3$$

Linearisation of the variation form of the kinetic energy yields the mass and gyroscopic matrices:

(20) $$\underline{{\underline{M} }} _{G}^{{i,j}} {\equals}\underline{{\underline{E} }} ^{T} \left[ {\matrix{ {\left( {\mathop{\int}_{\Omega _{e} } N^{i} N^{j} m_{{\rm \rho }} {\rm d\Omega }_{e} } \right)\underline{{\underline{I} }} _{3} } & {\underline{0} _{3} } \cr {\underline{0} _{3} } & {\mathop{\int}_{\Omega _{e} } N^{i} N^{j} \underline{{\underline{I} }} _{\rho } \underline{{\underline{T} }} _{s}^{{{\minus}1}} (\underline{{\rm \rtheta }} ){\rm d\Omega }_{e} } \cr } } \right]\underline{{\underline{E} }} $$

(21) $$\underline{{\underline{C} }} _{{k,G}}^{{i,j}} {\equals}\underline{{\underline{E} }} ^{T} \left[ {\matrix{ {\underline{{\underline{0} }} _{3} } & {\underline{0} _{3} } \cr {\underline{0} _{3} } & {\mathop{\int}_{\Omega _{e} } N^{i} N^{j} (\tilde \dot \underline {\rtheta } \underline{{\underline{I} }} _{{\rm \rho }} {\minus}\widetilde{{\underline{{\underline{I} }} _{{\rm \rho }} \dot \underline{{\rm \rtheta }} }})\underline{{\underline{T} }} _{s}^{{{\minus}1}} (\underline{{\rm \rtheta }} ){\rm d\Omega }_{e} } \cr } } \right]\underline{{\underline{E} }} $$

3.3 Governing equation for time-transient analysis

To solve the non-linear dynamic equation, the Hilbert Hughes Taylor (HHT)-α method is used. The resulting algorithm is of the predictor-corrector type, similar to that used in an earlier study⁽ Reference Crisfield ^{44 ,} Reference Le, Battini and Hjiaj ^{45 )}. In non-linear time-transient structural analyses, the displacements, velocities and accelerations are obtained iteratively at each time-step, i.e. by a Newton–Raphson scheme.

Linearisation of the inertial load vector is required to obtain the tangent inertial matrix. The final form of the tangent dynamic stiffness matrix is expressed as

(22) $$\underline{{\underline{K} }} _{{Dyn,G}}^{{\rm \alpha }} {\equals}{1 \over {{\rm \beta }h^{2} }}\underline{{\underline{M} }} _{G} {\plus}{{\rm \gamma } \over {{\rm \beta }h}}\underline{{\underline{C} }} _{{k,G}} {\plus}\underline{{\underline{K} }} _{G} $$

where $$\underline{{\underline{M} }} _{G} $$ , $$\underline{{\underline{C} }} _{{k,G}} $$ and $$\underline{{\underline{K} }} _{G} $$ respectively denote the mass, and gyroscopic and stiffness matrices, with respect to the global co-ordinate.

The governing equation of a moving structure can be expressed as

(23) $$(1{\plus}{\rm \alpha })\underline{f} _{e}^{{n{\plus}1}} {\minus}(1{\plus}{\rm \alpha })\underline{f} _{G}^{{n{\plus}1}} {\minus}\underline{f} _{{Km,G}}^{{n{\plus}1}} {\plus}{\rm \alpha }\left( {\underline{f} _{G}^{{,n}} {\minus}\underline{f} _{e}^{n} } \right){\equals}\underline{0} $$

where $$\underline{f} _{e} $$ , $$\underline{f} _{G} $$ and $$\underline{f} _{{Km,G}} $$ denote the external, internal and inertial force vectors, respectively. The prescribed motion is now considered to be the relevant acceleration of the motion:

(24) $$\underline{f} _{{Km,G}} {\equals}\underline{f} _{{k,G}} {\minus}\underline{{\underline{M} }} _{G} \underline{a} _{m} $$

where $$\underline{a} _{m} $$ denotes the acceleration of the motion. At the element level, $$\underline{a} _{m} $$ can be expressed as follows:

(25) $$\underline{a} _{m} {\equals}\left\{ {\matrix{ {\underline{a} _{{m,u}}^{T} } & {\underline{a} _{{m,{\rm \rtheta }}}^{T} } \cr } } \right\}^{T} $$

where $$\underline{a} _{{m,u}} $$ and $$\underline{a} _{{m,{\rm \rtheta }}} $$ are defined by the relevant acceleration of the translational ( $$\underline{{\ddot{u}}} _{p} $$ ) and rotational ( $$\underline{{{\ddot{\rtheta }}}} _{p} $$ ) motions, respectively:

(26a) $$\underline{a} _{{m,u}} {\equals}\underline{{{\tilde{\ddot{\rtheta }}}}} _{p} \underline{X} _{p} {\plus}\underline{{\ddot{u}}} _{p} $$

(26b) $$\underline{a} _{{m,\rtheta }} {\equals}\underline{{{\ddot{\rtheta }}}} _{p} $$

In the following sections, the FSI analysis is validated using the results of Heathcote’s experiment.

4.1 Description of experiment

The present FSI framework is validated by comparison with the results observed in the experiment conducted by Heathcote et al⁽ Reference Heathcote, Wang and Gursul ^{6 )}. Heathcote’s experiment considered rectangular cantilevered wings with different degrees of flexibility. The three types of wings (rigid, flexible and highly flexible) possessed different degrees of spanwise flexibility under a prescribed plunging motion.

Table 1 lists the operating conditions used in the present analysis. The reduced frequency k _G was varied from 0 to 1.82. k _G is defined by $$k_{G} {\equals}{{\pi f_{m} c} \over {U_{\infty} }}$$ , accounting for plunging frequency f _m , chord length c and free-stream velocity U _∞. The rigid, flexible and highly flexible wings were designed to be stiff, intermediately flexible and extremely flexible, respectively. In this study, the flexible and highly flexible wings are analysed in the present FSI analysis. The flexible and highly flexible wings were fabricated using polydimethylsiloxane rubber (PDMS) cast in an NACA0012 mold. The flexible wing was additionally reinforced using s 1-mm stainless steel plate, whereas 1-mm aluminium plate was used for the highly flexible wing. Cross-sections of the relevant configurations are shown in Fig. 2.

Table 1 Operating conditions

Figure 2 Cross-section of the flexible and highly flexible wings⁽ Reference Heathcote, Wang and Gursul ^{6 )}.

After the grid sensitivity analysis is completed, 389,400 grid cells and 29,000 surface boundary nodes are chosen in the current fluid analysis, and 936 elements are employed in the present shell analysis. A complete 3D wing configuration is considered in the fluid analysis, although the 3D configuration is reduced to a 2D plate for the shell analysis. The detailed discretisation in fluid and structural analyses is shown in Fig. 3, and the material properties of the flexible and highly flexible wings are provided in Table 2.

Figure 3 Discretisation information in the present FSI analysis.

Table 2 Structural properties used in the previous studies⁽ Reference Chimakurthi ^{14 ,} Reference Gordnier, Chimakurthi, Cesnik and Attar ^{19 )}

4.2 Validation with previous experimental results (k _G =1.82)

This subsection compares the present results with those obtained from previous studies and experiments. Mainly, the flexible and highly flexible wings are considered. First, the previous numerical and experimental results are compared with the present result when k _G is equal to 1.82. Figure 4 shows a comparison of the thrust coefficient histories. The average values for the thrust coefficient are then compared, and a summary of the relevant values is presented in Table 3

Figure 4 Comparison of the thrust coefficient history.

Table 3 Summary of averaged values of the thrust coefficient history

. The differences are obtained by subtracting the relevant physical quantities (Table 3).

For the flexible wing, the present thrust coefficient history shows good agreement with that observed in the experiment (Fig. 4). Further, the present prediction regarding the highly flexible wing shows a trend similar to that predicted by Gordnier et al using a higher-order fluid dynamic analysis (Fig. 4). However, the high-frequency response observed in the experiment was not captured in either prediction. In neither the present nor previous studies⁽ Reference Chimakurthi ^{14 ,} Reference Gordnier, Chimakurthi, Cesnik and Attar ^{19 )}, physical reason regarding the high-frequency response is obtained. Extensive investigation of such FSI phenomena in the highly flexible wing will still be required.

The wing tip displacement history is illustrated in Fig. 5. The tip displacement z _tip is normalised by the amplitude of the plunge motion z _o . For the flexible wing, the present prediction exhibits a phase difference compared with the previous numerical and experimental results. However, the present result still correlates reasonably with the experimental results. Regarding the structural flexibility, a phase variation in the tip history is induced as a result of inertial and aerodynamic forces, which becomes more significant for the highly flexible wing. The relative differences in the peak-to-peak values and degrees of phase variation in the tip displacement history are summarised in Table 4.

Figure 5 Comparison of the wing tip displacement history.

Table 4 Summary of peak-to-peak value and phase difference in the wing tip displacement history

Then, the pressure coefficient C _P on the wing surface is presented in order to evaluate the influence of the wing flexibility in detail. The relevant comparison is conducted in situations with 0.0 t/T, 0.25 t/T, 0.5 t/T and 0.75 t/T. The relevant comparisons between the present results for the flexible and highly flexible wings are presented in Figs 6 and 7, respectively. Figure 6 shows the pressure distribution of the flexible wing, and it is found that the pressure distribution is uniform in the spanwise direction in all situations. For the highly flexible wing, the variation of the pressure distribution in the spanwise direction is significant (Fig. 7). This phenomenon is due to the phase-lagging motion of the wing tip. Thus, lower and higher pressure regions exist simultaneously around the leading edge. Specifically, the difference in the pressure distribution along the spanwise direction is significant for 0.0 t/T and 0.5 t/T. As a result, degradation of the thrust coefficient is shown in the highly flexible wing. This phenomenon was also predicted by the previous studies, and it was found that an intermediate flexibility may be beneficial for increasing the thrust.

Figure 6 Pressure coefficient, C _P, on the surface of the flexible wing (k _G=1.82).

Figure 7 Pressure coefficient, C _P, on the surface of the highly flexible wing (k _G=1.82).