Study on Flow-induced Vibration Characteristics of 2-DOF Hydrofoil Based on Fluid-Structure Coupling Method

Yichen Jiang, Chuansheng Wang, Jingguang Li, Chunxu Wang and Qing Wang

Abstract The flutter of a hydrofoil can cause structural damage and failure, which is a dangerous situation that must be avoided.In this work, based on computational fluid dynamics and structural finite element methods, a co-simulation framework for the flow-induced vibration of hydrofoil was established to realize fluid-structure interaction.Numerical simulation research was conducted on the flow-induced vibration characteristics of rigid hydrofoil with 2-DOF under uniform flow, and the heave and pitch vibration responses of hydrofoil were simulated.The purpose is to capture the instability of hydrofoil vibration and evaluate the influence of natural frequency ratio and inertia radius on vibration state to avoid the generation of flutter.The results indicate that when the inflow velocity increases to a certain critical value, the hydrofoil will enter the flutter critical state without amplitude attenuation.The attack angle of a hydrofoil has a significant impact on the vibration amplitude of heave and pitch.Additionally, the natural frequency ratio and inertia radius of the hydrofoil significantly affect the critical velocity of the flutter.Adjusting the natural frequency ratio by reducing the vertical stiffness or increasing the pitch stiffness can move the vibration from a critical state to a convergent state.

Keywords Hydrofoil; Flutter; Flow-induced vibration; Fluid-structure interaction; Critical velocity

1 Introduction

Hydrofoil, as a typical lifting body, is widely used in the fields of ships and hydraulic machinery.Marine rudders,fin stabilizers, planing wings of hydrofoils, various control surfaces of underwater vehicles, and blades of propellers and hydraulic turbines are all hydrofoil structures.The use of hydrofoil technology also provides greater maneuvering stability, particularly with respect to changes in wave ele‐vation (Mahmud, 2015).Because a hydrofoil is an elastic body, it can vibrate and deform under the action of fluid loads.Due to the complex and variable flow field environ‐ment in which it works, its vibration response can also change significantly, with strong unsteady characteristics.If the design of hydrofoils is not reasonable, the fluid may stimulate hydrofoils to flutter as the speed increases, threat‐ening the safety of the structure.Flutter is a self-excited vibration phenomenon, in which the fluid excites the hy‐drofoil structure to vibrate, and the vibration in turn can disturb the flow field and change the distribution character‐istics of the load.If the changing fluid load happens to have a negative damping effect on the hydrofoil, the fluid will not contribute to structural stability.At this time, the system will continue to obtain energy from the fluid, and instability with gradually expanding amplitude will occur,ultimately causing structural damage and failure, which is a dangerous situation that needs to be avoided.Uncontrolled static and/or dynamic instabilities can lead to excessive deformation, material failure, accelerated fatigue, vibration,noise issues, and even catastrophic structural failures(Abramson, 1969; A.Kousen and O.Bendiksen, 1988; O.Bendiksen, 1992; Bendiksen, 2002).

The early research on flow-induced vibration of airfoil structures was mainly focused on the field of aviation en‐gineering, which was used to guide the aeroelastic design of aircraft and solve wing flutter problems.Researchers mainly used linear potential flow theory to analyze and study the aerodynamic characteristics of two-dimensional elastic wings, and had obtained some important theoreti‐cal results, especially the unsteady aerodynamic theory proposed by Theodorsen (Theodorsen, 1935; Theodorsen and Garrick, 1940), which laid the foundation for the aero‐elastic stability and flutter analysis of airfoil structures.Lottati et al.(Lottati, 1985; Smith and Chopra, 1990; Wu and Sun, 1991; White and Heppler, 1995; Banerjee, 2001)solved the flutter characteristics of box beams, Timoshen‐ko beams (vertically and horizontally curved beams), and thin-walled beams using spline methods based on The‐odorsen’s aeroelastic theory and wing theory, verifying the wide applicability of this method.Zhang (Zhanget al., 2019) proposed a low-cost and low-risk flutter bound‐ary prediction method, which requires only one upwind test at subcritical velocity.Even if the reference dynamic pressure is far from the flutter critical point, it can also have a certain prediction accuracy for flutter occurrence.Jonsson (Jonssonet al., 2019) reviewed the methods of flutter analysis in aircraft design optimization and summa‐rized the research status and main problems in this field.In the field of hydrodynamics, with the gradual develop‐ment of ships towards high speed and large scale, the hy‐drodynamic loads borne by structures continue to in‐crease, placing higher requirements for structural safety and reliability.The vibration and hydroelasticity of under‐water structures have also become noticeable, especially for various prominent hydrofoil appendages.Therefore,in recent years, more and more research has been conducted on the flow-induced vibration of hydrofoils.Currently, the main methods include experimental testing and numerical calculation.

In terms of experiments, the following scholars have conducted relevant research work.Taking the rudder shaft system as the research object, Jewell (Jewell and McCor‐mick, 1961) had demonstrated through a large number of experiments that the Theodorsen 2-DOF motion equation,which is widely used in aerodynamics, is also applicable to underwater lifting bodies.Ducoin (Ducoinet al., 2012)conducted a water tunnel test on the structural response of a cantilever flexible hydrofoil with or without cavitation,and measured the displacement deformation and surface vibration velocity of the free end of the hydrofoil using a high-speed camera and a laser doppler vibrometer.The re‐search showed that the response of the hydrofoil under non-cavitation flow is closely related to the hydrodynamic load affected by viscosity, and the unsteady hydrodynam‐ic force caused by cavitation can significantly increase the vibration level.At the same time, it will cause the at‐tached water mass to decrease and affect the response fre‐quency of some modes of the structure.Harwood (Har‐woodet al., 2019a; 2019b) conducted experiments on a rigid variant and two flexible variants of a vertical sur‐face-piercing hydrofoil in wetted, ventilating and cavitat‐ing conditions, describing the dynamic FSI response of a flexible surface-piercing hydrofoil in dry, wetted, ventilat‐ing and cavitating conditions.Smith (Smithet al., 2020a;2020b, Younget al., 2022) used high-speed camera and synchronous force measurements to study the effects of fluid-structure interaction on rigid and flexible hydrofoil cloud cavitation in a variable pressure water tunnel and explained the complex response.Herath (Herathet al.,2021) conducted cavitation tunnel tests on optimized shape adaptive composite hydrofoils and compared the re‐sults with other composite hydrofoils.The study found that hydrofoils can have significantly different responses for the same Reynolds number if the flow velocity is dif‐ferent.Liu (Liuet al., 2021; 2023) used a synchronous measurement system consisting of a high-speed camera, a hydrodynamic load cell, and a Laser Doppler Vibrometer to study the dynamics of sheet/cloud cavitation in fluidstructure interaction, and summarized the effects of bendtwist coupling on the cavitation behavior and hydroelastic response of composite hydrofoils.

Compared to experimental testing, using numerical methods to predict the hydroelastic response of hydrofoils can save manpower and material resources, and obtain more flow field and vibration information, which is helpful to reveal the flow-induced vibration laws of hydrofoils.With the emergence of workstations and large-scale computers,fluid-structure interaction numerical methods combining computational fluid dynamics (CFD) and computational structural dynamics (CSD) have become an effective means to solve vibration problems of underwater structures.Kin‐nas (Kinnas and Fine, 1993, Senocak and Wei, 2002) con‐ducted a large number of numerical studies on the cavita‐tion response of rigid hydrofoils, but their research on flex‐ible hydrofoils was limited.Ducoin (Ducoinet al., 2010)solved the fluid-structure interaction problem of a hydro‐foil with 2-DOF, bending and twisting, based on the com‐mercial computational fluid dynamics software STARCCM+.The results showed that the motion of the hydrofoil increases the angle of attack, leading to an increase in lift and drag coefficients and an increase in cavitation length.Chae (Chaeet al., 2013) coupled the unsteady RANS solv‐er with a two-degrees-of-freedom (2-DOF) structural model to explore the effects of relative mass ratio and viscous ef‐fects on the fluid-structure interaction response and stability of hydrofoils.The results showed that the instability modes of hydrofoils are divergent, dynamic divergent, and flutter in different relative mass ratio regions.Ducoin (Ducoin and Young, 2013) simplified the flexible hydrofoil into a 2-DOF system and calculated the static divergence velocity using CFX.The results showed that the viscous effect moves the pressure center toward the middle of the hydrofoil at large angles of attack, which helps to delay or suppress di‐vergence.Wang (Wang et al., 2009) examined the effect of different locations of the pitching axis on propulsive perfor‐mance of a two-dimensional flexible fin, the results showed that maximum thrust can be achieved with the pitching axis at the trailing edge.The static fluid-structure interaction simulation results of Huang (Huanget al., 2019) showed that compared to unidirectional fluid-structure interaction, the flexible hydrofoil has a lower lift and greater resistance when subjected to bidirectional fluid-structure interac‐tion.Liu (Liuet al., 2019) used dynamic mode decompo‐sition (DMD) and inherent orthogonal decomposition(POD) to analyze the coherent structure of cavitation flow over ALE-15 hydrofoils.The results showed that DMD is more effective in decomposing complex flow fields into uncoupled coherent structures with specific dy‐namic modes and corresponding frequencies than POD.Xu (Xu and Tang, 2020) simulated the flow field of a NA‐CA66-mod hydrofoil under different trailing edge vibra‐tion amplitudes and frequencies, and the results showed that the vibration amplitudes and frequencies have a sig‐nificant impact on trailing edge vortex shedding.Chae(Chae and Young, 2021) discussed the effects of key geo‐metric parameters, material parameters, and flow parame‐ters on the steady-state and dynamic responses of span‐wise flexible airfoils and hydrofoils.The results showed that the main instability of hydrofoils is divergence, while the main instability of airfoils is flutter.Hou (Hou et al.,2013) studied the effects of geometric parameters of the rudder on the hydrodynamic performance of the propellerrudder system, the calculation results showed that the rud‐der span has an optimal match range with the propeller di‐ameter.George (George and Ducoin, 2021) used the di‐rect numerical simulation (DNS) method to simulate three-dimensional incompressible flow, and the results showed that vibration has a significant impact on the pres‐sure gradient of hydrofoils, thereby changing the transi‐tion position and boundary layer characteristics that are proportional to the vibration amplitude and frequency ra‐tio.Negi (Negiet al., 2021) conducted a global linear sta‐bility analysis of the fluid-structure interaction problem of NACA0012 airfoil at transition Reynolds numbers and studied the spontaneous pitch oscillation phenomenon.The results showed that the aeroelastic instability is caused by the linear instability mode of the fluid-structure inter‐action problem.Wang (Wanget al., 2021) presented a multi-scale Eulerian-Lagrangian approach to simulate cavitating turbulent flow around a Clark-Y hydrofoil to study the bubble dynamics.The results showed that the pressure wave released as a bubble is compressed reaches 107Pa, which may cause cavitation erosion of the hydro‐foil surface.

Due to the tremendous damage that the flutter can cause to structures, many studies have been conducted on the flow-induced vibration of hydrofoils.However, most of the current research focuses on the flow-induced vibration char‐acteristics of flexible hydrofoils or the cavitation response of rigid hydrofoils, while there are few studies on the flowinduced vibration characteristics of rigid hydrofoils.At the same time, theoretical formulas are often used to derive the flutter critical velocity, and corresponding numerical stud‐ies are lacking.In response to these two shortcomings, the research goal of this article is to study the flow-induced vi‐bration characteristics of a rigid hydrofoil with 2-DOF un‐der uniform flow through numerical simulation, in order to avoid the occurrence of flutter.This work complements pre‐vious research and has a certain reference value.Firstly, the flow field model and structural model of this numerical sim‐ulation are introduced; Then, the dynamic responses under different velocities and angles of attack are analyzed, and the effects of natural frequency ratio and inertia radius on the vibration state of hydrofoils are compared; Finally, the work of this article is summarized, and the shortcomings and directions for further research are pointed out.

2 Numerical model

Based on the co-simulation of computational fluid dy‐namics (CFD) and structural finite element method(FEM), this work performs a bidirectional fluid-structure interaction calculation and conducts a numerical simula‐tion study on the flow-induced vibration of a three-dimen‐sional elastic hydrofoil to obtain the dynamic response characteristics of the hydrofoil at different velocities and angles of attack.At each time step, the three-dimensional hydrodynamic load is calculated by the CFD numerical model, and the pressure and shear force are imported into the FEM through the fluid-structure interface to solve the structural response.The three-dimensional structural dis‐placement is then imported into the flow field to achieve the update of the flow field grid, and the hydrodynamic load and structural response are recalculated.The calcula‐tion repeats the iteration until the solution variable con‐verges or the maximum internal iteration number is reached, and then the next time step is started.The fluidstructure coupling calculation process is shown in Figure 1.

Section 2.1 verifies the effectiveness of this calculation method using the NACA66 flexible hydrofoils.Subse‐quently, this method is further applied to calculate 2-DOF NACA66 rigid hydrofoils.Sections 2.2 and 2.3 respective‐ly describe the flow field model and structural model of the 2-DOF NACA66 rigid hydrofoil used in this work.

2.1 Numerical model validation

Figure 1 The fluid-structure coupling calculation process

Ducoin conducted flow-induced tests on NACA66 flexi‐ble hydrofoils at different cavitation numbers in a water tunnel (Ducoinet al., 2012).This section verifies the effec‐tiveness of the fluid-structure interaction method concern‐ing non-cavitation conditions.The chord length of hydro‐foil is 0.148 m, and the extension length is 0.191 m.One end is fixed to the wall of the water tunnel, and the other end is free.The cross-sectional dimension of the water tun‐nel is 0.192 m×0.192 m.Based on Star-CCM+ and Abaqus,a bidirectional fluid-structure interaction calculation was conducted for the NACA66 hydrofoil.The fluid calculation domain was 1 m long, with a velocity inlet at the left end, a pressure outlet at the right end, and a non-slip wall around.Large eddy simulation was used for calculation.The mesh is divided using a cutting volume method.The dimension‐less thickness of the first layer of mesh near the wall meetsy+<1, and the total number of meshes is 3.25 million.The flow field mesh scene is shown in Figure 2.

Figure 2 Flow field grid

The hydrofoil material is polyoxymethylene, with a den‐sity ofρ=1 480 kg/m3, Young’s modulus ofE=3.0 GPa, and a Poisson’s ratio ofϑ=0.35.The mesh is discretized using a second-order hexahedral element C3D20R.The finite ele‐ment model is shown in Figure 3.The natural frequency and modal calculation results of the hydrofoil are shown in Table 1 and Figure 4.It can be seen that the vibration modes of the hydrofoil in vacuum and in water are basically the same, with the first three vibration modes being first-order bending, first-order pitch, and second-order bending.How‐ever, due to the high density of water, the natural frequen‐cies in water are significantly lower than those in real air.Comparing the wet mode natural frequencies with the re‐sults given by Ducoin, it can be seen that the two modes are basically consistent.

Figure 3 NACA66 flexible hydrofoil finite element model

Table 1 First three natural frequencies (Hz)

After that, a fluid-structure interaction calculation was performed, with an incoming flow velocity of 5 m/s, and the structure was solved using implicit dynamic solutions.The vertical displacement of the monitoring point at the free end guide edge of the hydrofoil at 4°, 6°, and 8° an‐gles of attack is shown in Figure 5.It can be seen that the vibration gradually attenuates and eventually fluctuates slightly around a constant value.The larger the angle of attack, the higher the equilibrium position.The mean val‐ue of the stable section is in good agreement with the ex‐perimental values given in the literature (Ducoin and Young, 2013), which verifies the effectiveness of the fluidstructure coupling calculation method in this work.

Figure 4 First three modal shapes

Based on the aforementioned calculation method, this work also compared the lift coefficients with the experimen‐tal data.The experimental data were obtained from testing a hydrofoil with the NACA66 series, having a chord length of 0.150 m and a span of 0.191 m, at a nominal freestream velocity of 5.33 m/s and at 4°, 6°, and 8° angles of attack.The measured and predicted lift and drag coefficients of the hydrofoil at these angles are given in Table 2.The numeri‐cal results are consistent with the experimental data given in the literature (Lerouxet al., 2004), which verifies the ef‐fectiveness of the numerical model in this work.

Figure 5 Comparison between calculated displacement value and test value

Table 2 Comparison of lift and drag coefficients with test data

2.2 Flow field model

The hydrofoil has a section shape of NACA66 airfoil with a chord length ofc=0.148 m.The fluid density is 997.561 kg/m3and the kinematic viscosity coefficient is 1.012 5 mm2/s.The upstream and upper/lower boundaries of the flow field are set as velocity inlets, and the down‐stream boundary is set as a pressure outlet.The ratio of the spanwise width of the flow field to the chord length of the hydrofoil is 1∶3.Both sides are taken as symmetrical planes to reduce the number of grids in the flow field and structure and the computational cost of fluid-structure in‐teraction.The grid is divided using a cutting volume meth‐od and locally densified.The mesh near the wall meets the requirements ofy+<1.Change the basic size of the grid to divide three sets of flow field grids, with grid amounts of 2.08 million, 3.17 million, and 7 million, respectively.The 3.17 million grid case is shown in Figure 6.The grid inde‐pendence verification of the flow field is performed at a 6°angle of attack and a 5 m/s flow velocity.The time step is 4×10−4s.The SSTk−ωDES turbulence model is used for calculation, and the average value of the hydrofoil lift co‐efficient is monitored.The results are shown in Table 3.The difference percentage of the lift coefficient is defined as (Cl−Clm)/Clm, whereClmis the lift coefficient of the maximum grid quantity model.It can be seen that the lift coefficient basically converges when the grid quantity is 3.17 million and 7 million.

Figure 6 2-DOF hydrofoil flow field grid

Table 3 Results of lift coefficients for hydrofoils with different grid quantities

Based on a comparison of the calculation results of three time steps 2×10−4s, 4×10−4s, and 8×10−4s based on a 3.17 million grid, as shown in Table 4, it can be seen that the lift coefficient difference is small when the time steps are 2×10−4s and 4×10−4s.Considering the calcula‐tion accuracy and time length, a 3.17 million grid com‐bined with a time step of 4×10−4s is selected for subse‐quent calculation.

Table 4 Lift coefficient results of hydrofoils at different time steps

2.3 Structural model

As shown in Section 2.1, the numerical model is validated based on three-dimensional flow.In the following investi‐gations, the flow fields are also chosen to be three dimen‐sional.It is worth pointing out that the simulated hydrofoil has the same NACA66 cross section, however, it is rigid body and moves with heave and pitch motions only, as shown in Figure 7.It is simulated to represent a section of the 3-DOF hydrofoil in order to focus on the interaction between the hydrofoil vibration and the vortex shedding without the distraction of the three-dimensional deforma‐tion, as well as reduce the computational costs of multiple simulations.

Figure 7 2-DOF rigid hydrofoil

The elastic axis is located at one third of the chord length from the leading edge, and a center of gravity CG located at the geometric center between the elastic axis and the trailing edge.xis the horizontal dimensionless distance from the elastic axis,SandIare the static mo‐ment and moment of inertia of the hydrofoil about the elastic axis.The natural frequencies of pure heave and pure pitch motions of hydrofoils are the same, both are 9.5 Hz.The calculation formula for relevant parameters is as follows:

The relevant operating data and the geometrical of the hydrofoil are shown in Table 5.

Table 5 Operating data and geometrical data of hydrofoil

The vibration equation for a 3-dimensional rigid wing with 2-DOF can be expressed as:

whereMS,CS, andKSare, respectively, the structural iner‐tial, damping, and stiffness matrices, whileX,Ẋ, andẌ are, respectively, the structural displacement, velocity, and acceleration vectors withX=[h,θ]T.Furthermore,F=[L,M]Tis the fluid force vector.

wherem,Sθ, andIθare, respectively, the structural mass,the static imbalance, and the structural moment of inertia defined about the E.A., whileChandCθ,KhandKθare, re‐spectively, the structural damping values for the bending and twisting motions, the structural bending and torsional stiffness values.The overhead dots indicate Lagrangian time derivatives and henceh,ḣ, andḧ are, respectively, the bending displacement, velocity, and acceleration, whileθ,θ̇, andθ̈ are, respectively, the twisting displacement, veloc‐ity, and acceleration.LandMare the fluid lift and moment acting on the foil, defining the heave motion as positive upward and the pitch motion as positive counterclockwise.If the right side of the equation is 0, and the hydrofoil pa‐rameters are substituted into the equation, the natural fre‐quencies of the 2-DOF undamped system can be obtained to be 7.77 Hz and 12.95 Hz, respectively.

When hydrofoil flutter occurs under fluid excitation, its vibration state approximates simple harmonic motion.The‐odorsen has derived the unsteady forces and moments on two-dimensional thin wings that vibrate harmonically at the angular frequencyωin incompressible uniform flow.Based on this theory, the flutter critical velocities of 2-DOF hydrofoils are predicted.The derivation of the unsteady forces and moments is relatively complex.Here, the expres‐sions for the unsteady liftLand momentMof a hydrofoil with a span oflare directly given:

whereρis the density of the medium,Vis the inflow ve‐locity,ais the dimensionless distance of the hydrofoil midpoint from the elastic axis, backward is positive, the reduc‐tion frequencyk, andC(k) is the Theodorsen function.

whereJ0,J1,Y0,Y1is the first and second type of standard Bessel functions ofk, respectively.

Substituting Eq.(7) into (4) and (5) can obtain:

The expression ofLh,Lθ,Mh,Mθin the formula are:

To make the equations solvable, the determinant of the left coefficient matrix is 0.However, the coefficient matrix contains two unknown variables, the reduced frequencykand the vibration angular frequencyω.kexists in the coeffi‐cientLh,Lθ,Mh,Mθand cannot be directly solved.Here, the V-g method is used to solve the flutter determinant.Artifi‐cial damping is introduced into the vibration Eq.(3), so that the damping forces on the 2-DOF are −igkhhand −igkθθ, re‐spectively.At this time, the flutter determinant becomes:

Here, we can assume a series ofkvalues, for eachk, we can obtain a set ofLh,Lθ,Mh,Mθvalues, and solve the two complex rootsZ=ZR+iZIaboutZ.At the same time, we can determine the vibration angular frequencyω, artificial dampingg, and inflow velocityVbased on the following relationships:

Plot each set ofVandgvalues obtained, with the horizon‐tal axis beingVand the vertical axis beingg.Wheng<0, it is necessary to provide negative damping to maintain simple harmonic vibration of the hydrofoil, that is, to input energy to the system, indicating that flutter will not occur in the hydrofoil vibration at this flow velocity; Wheng>0, it is necessary to provide damping to consume system energy to generate simple harmonic vibration, indicating that flutter has occurred on the hydrofoil at this flow velocity; Wheng=0, the hydrofoil can maintain simple harmonic vibration only under fluid excitation, and the correspondingVis the critical velocity of flutter.

Next, the parameters of the hydrofoil structural model are substituted into the flutter Eq.(12) and calculated using the V-g method.Given a series ofkvalues, the corresponding V-g curve is obtained, as shown in Figure 8.It can be seen that as the inflow velocity increases, the value ofgchanges from a negative value to a positive value, and the trend of change decreases first and then increases.Wheng=0, the corresponding inflow velocityV=4.5 m/s, so the flutter criti‐cal velocity calculated by the V-g method is 4.5 m/s.

Figure 8 V-g figure

3 Results and discussion

This chapter focuses on analyzing the heave displace‐ment, pitch angle, hydrodynamic force, and real-time ener‐gy exchange between fluid and elastic systems, analyzing the dynamic response under different flow velocities and angles of attack, and comparing the effects of hydrofoil nat‐ural frequency ratios and inertial radius on vibration states.

3.1 Dynamic response of hydrofoil at different velocities

This section presents the dynamic response of a 2-DOF hydrofoil at velocitiesV=1 m/s, 3 m/s, 4.5 m/s, 5 m/s, and 6 m/s, with an initial angle of attack of 6°.WhenV=1 m/s,the curves of heave displacement, pitch angle, lift, and the torque of the fluid to the elastic axis are shown in Figure 9.It can be seen that the amplitude of heave and pitch motion attenuates with time.The spectrum results in Figure 10 shows that heave and pitch vibration contain two main fre‐quency components, 7.3 Hz and 12.2 Hz, respectively,which are close to the first-order and second-order natural frequencies of the undamped system.The vibration is greatly affected by the natural characteristics of the hydro‐foil.The non-dimensional Froude number, vibration state,maximum power and average work for each velocity are shown in Table 6.The maximum power and average work are taken from the last cycle of the critical state, and the calculation formula for Froude number is as follows:

whereVis the inflow velocity,gvalue is 9.81 m/s2, andcis the chord length of hydrofoil.

Figure 9 Vibration and hydrodynamic force time history results(1 m/s)

Figure 10 Heave and pitch motion spectrum results (1 m/s)

Table 6 Maximum power and average work for each velocity

To observe the flow field and the energy exchange of a 2-DOF elastic system, calculate the hydrodynamic powerP(t) and the total workW(t) it does to the system at differ‐ent times.The results are shown in Figure 11.It can be seen that the amplitude of hydrodynamic power and the to‐tal work done to the elastic system have a downward trend,indicating that the intensity of energy exchange between the fluid and the hydrofoil has decreased, and the energy of the elastic system has a downward trend, at which time the vibration of the hydrofoil is in a convergent state.

whereḣ andθ̇ are the heave velocities and pitch angular velocities, and the lift and fluid torque received by the elas‐tic axis areL(t) andM(t).

Figure 11 Hydrodynamic power and work history results (1 m/s)

The dynamic response of the hydrofoil atV=3 m/s is shown in Figures 12‒14.Compared toV=1 m/s, the ampli‐tudes of heave and pitch motions significantly increase,but the attenuation trend decreases.In addition, there is on‐ly one main frequency component remaining in the vibra‐tion spectrum, 12.2 Hz, and the amplitudes of other fre‐quency components are small.The vibration of the system starts to be dominated by a single frequency.At this point,the amplitude of the hydrodynamic power does not signifi‐cantly attenuate, while the maximum amplitude of the to‐tal work of the fluid still decreases from 0.4 seconds, indi‐cating that as the number of vibration cycles increases, the fluid will dissipate the energy of the system, and the vibra‐tion is still in a convergent state.

The time history results of hydrofoil vibration and hydro‐dynamic work atV=4.5 m/s are shown in Figures 15 and 16.The amplitudes of heave displacement, pitch angle, and lift and torque decay rapidly over time.According to the average time interval of vibration peaks, the vibration fre‐quencies of both degrees of freedom are 11.8 Hz.The nega‐tive work of hydrodynamic force in each vibration cycle is higher than the positive work, and the overall damping ef‐fect is negative.The system energy continuously decreases,and the vibration is in a convergent state.

Figure 12 Vibration and hydrodynamic force time history results(3 m/s)

Figure 13 Heave and pitch motion spectrum results (3 m/s)

When the inflow velocityVincreases to 5 m/s, the dy‐namic response of the hydrofoil is shown in Figures 17‒19.It can be seen that the amplitudes of heave and pitch mo‐tions continue to increase without a trend of attenuation,and the vibrations of the 2-DOF are still dominated by a single frequency, both 11.6 Hz.At the same time, the am‐plitude of hydrodynamic power no longer attenuates, and the total work of the fluid presents a slightly increasing trend with the increase of the vibration period.The energy input to the system begins to accumulate, which determines that the 5 m/s flow velocity of the hydrofoil is in an unsta‐ble state, that is, a critical state of flutter.It can be seen that the flutter critical velocity of the numerical simulation is higher than the theoretical calculation value of 4.5 m/s based on the Theodorsen unsteady force in Section 2.3.This is mainly because the derivation of the Theodorsen un‐steady force is based on the potential flow method, with‐out considering viscosity and specific hydrofoil shape, and the expression is a solution in the frequency domain.Al‐though the numerical simulation takes a long time to calcu‐late, it can fully consider the hydrofoil shape, In the time domain, the unsteady force and hydrofoil response are directly solved based on the CFD method.

Figure 14 Hydrodynamic power and work history results (3 m/s)

Figure 15 Vibration and hydrodynamic force time history results(4.5 m/s)

Figure 16 Hydrodynamic power and work history results (4.5 m/s)

Figure 17 Vibration and hydrodynamic force time history results(5 m/s)

Figure 18 Heave and pitch motion spectrum results (5 m/s)

Figure 19 Hydrodynamic power and work history results (5 m/s)

Continue to increase the incoming flow velocity, and the results atV=6 m/s are shown in Figures 20‒22.At this point, the vibration of the hydrofoil enters a flutter state,and the amplitudes of the heave and pitch motions gradually increase with time.The frequency spectrum shows that both the heave and pitch frequencies are 11 Hz, which is somewhat lower than the critical state.This may be due to the rapid increase in the amplitude of the hydrofoil in the flutter state and the greater impact of the attached fluid.In addition, the power amplitude and total work of this hydro‐dynamic force also multiply with the increase of the vibra‐tion period, and the system continuously obtains energy from the fluid until damage occurs.

Figure 20 Vibration and hydrodynamic force time history results(6 m/s)

The flow field vorticity results of the midspan cross section of the hydrofoil in the critical and flutter states are shown in Figure 23.Corresponding to the last complete cycleTin the hydrodynamic power curve, the hydrodynamic power is 0 at 0 andT, positive atT/5 and 2T/5, and nega‐tive at 3T/5 and 4T/5.The flow field vorticity results of the three-dimensional hydrofoil in the critical and flutter states at the final moment are shown in Figure 24.During the entire vibration cycle of the critical state, vortices are mainly generated in the middle and rear parts of the hydrofoil.With the constant change of the attack angle of the hydro‐foil, they are brought into the wake flow field by the main stream to form relatively fine shedding vortices.In the case of flutter, the vorticity of the suction surface is still concen‐trated in the rear region during the upstroke of the hydro‐foil.However, during the downstroke, a series of signifi‐cant vortices are generated in the middle and front of the suction surface, which are then dispersed by the incoming flow and merged into the vortex at the rear as the angle of attack decreases.This is mainly due to the large response of the hydrofoil to the heaving and pitching motions during the flutter state, and the rapid replenishment of surrounding fluid to the suction side surface during the downstroke.

Figure 21 Heave and pitch motion spectrum results (6 m/s)

Figure 22 Hydrodynamic power and work history results (6 m/s)

Figure 23 Flow field vorticity results of the midspan cross section of the hydrofoil (left: critical state, right: flutter state)

Figure 24 Flow field vorticity results of the three-dimensional hydrofoil

3.2 Dynamic response of hydrofoil at different angles of attack

Other conditions remain unchanged, changing the ini‐tial angle of attack to 4° and 8°, and calculating the flowinduced response at a flow velocity of 5 m/s.The vibra‐tion state, maximum power and average work for each an‐gle of attack are shown in Table 7.The vibration and hy‐drodynamic results at different angles of attack are shown in Figures 25 and 26.

Table 7 Maximum power and average work for each angle of attack

Figure 25 Vibration and hydrodynamic force time history results at different angles of attack (5 m/s)

After adjusting the angle of attack to 4° and 8°, the heave and pitch vibrations of the hydrofoil do not attenuate or expand, but maintain constant amplitude vibration.As the angle of attack increases, the average lift increases,and the oscillation amplitude and equilibrium position of the heave motion also significantly increase.Compared to 4° angle of attack, the pitch amplitude increases slightly at 6° angle of attack, while it increases significantly at 8°angle of attack, consistent with the variation of torque.In addition, the peak frequencies of heave and pitch vibra‐tions at 4°, 6°, and 8° attack angles are all around 11.6 Hz,indicating that changes in the initial attack angle have a small impact on the vibration frequency.Figure 27 shows the fluid work situation, and the results show that the oscil‐lation amplitude of the hydrodynamic power and work curve increases significantly with the increase of the angle of attack, and there is no attenuation trend and gradually tends to stabilize.Based on the vibration curve, it can be judged that the hydrofoil at 4° and 8° angles of attack is still in a critical state, that is, a small change in the initial angle of attack has no significant impact on the critical flutter speed of the hydrofoil.

Figure 26 Spectrum results of heave and pitch motions at different angles of attack (5 m/s)

The vorticity results of the flow field at angles of attack of 4° and 8° are shown in Figure 28, which also correspond to the last complete cycle of the hydrodynamic power curve.

It can be seen that at a 4° angle of attack, the vibration amplitude of the hydrofoil is relatively small, and there is a relatively obvious alternating trajectory of vortices in the wake field.Due to the influence of the hydrofoil motion and incoming flow, the tail vortices in the downstroke are longer and narrower than in the upstroke.At an angle of attack of 8°, the amplitude of the hydrofoil increases, and the vorticity intensity on the suction side increases and moves towards the leading edge.As the hydrofoil vibrates,the characteristics of alternating vortices emanating from the rear become weaker, and the vorticity distribution in the wake field is more finely distributed.

3.3 Influence of frequency ratio and inertia radius on vibration state

Figure 27 Hydrodynamic power and work history results at different angles of attack (5 m/s)

Figure 28 Vorticity results of the flow field at different angles of attack (left: 4°, right: 8°)

Based on the flutter critical velocity of the 6° angle of attack model, the effects of different natural characteristics of the elastic system on the vibration state are compared and calculated.The other parameters remain unchanged,and the heave and pitch stiffness of the elastic axis are adjusted respectively to change the natural frequency ratiofh/fθof the hydrofoil.Keeping the vertical stiffness con‐stant and increasing the pitch stiffness, a hydrofoil model offh/fθ= 0.845 is obtained; Keeping the pitch stiffness constant and reducing the vertical stiffness, a hydrofoil model offh/fθ= 0.707 is obtained.The vibration state, max‐imum power and average work for each natural frequency ratio are shown in Table 8.The flow-induced response re‐sults of hydrofoils with three frequency ratios at a flow velocity of 5 m/s are shown in Figure 29.

From the vibration curve, it can be seen that after increas‐ing the pitch stiffness to makefh/fθ=0.845, the vibration fre‐quencies of the hydrofoil heave and pitch motions at theinitial stage are higher than those of the original model, and then the amplitude attenuates and the vibration frequency decreases.The heave range of lift and torque also decreases,and the equilibrium position of the heave motions decreases.Whenfh/fθ=0.707, due to the decrease in vertical stiffness,the equilibrium position of the hydrofoil’s heave motion is lifted up, and the vibration periods at the 2-DOF are increased compared to the original model, the amplitude also shows an attenuation trend.Figure 30 shows that when the model frequency ratio is 0.707 and 0.845, the power amplitude of the hydrodynamic force and the total work it does gradually attenuate, and the vibrations of the two models are in a convergent state.It can be seen that increasing the pitch stiffness can effectively suppress flutter,and whenfh/fθapproaches 1, reducing the heave stiffness to stagger the heave natural frequencies from the pitch natural frequencies can also improve hydrofoil flutter.

Table 8 Maximum power and average work for each natural frequency ratio

Figure 29 Vibration and hydrodynamic force time history results at different frequencies (5 m/s)

Figure 30 Hydrodynamic power and work history results at different frequency ratios (5 m/s)

In addition, different structural designs can change the mass distribution of the hydrofoil, thereby changing the in‐ertia radius of the hydrofoil and affecting the inherent char‐acteristics of the system.The following maintains the mass, elastic axis position, stiffness, and centroid position of the hydrofoil unchanged at an angle of attack of 6°.By comparing and calculating the effects of different inertia radii on the flow-induced vibration of the hydrofoil, the dimensionless inertia radiusrof the three hydrofoil models are 0.45, 0.542, and 0.634, respectively.The correspond‐ing pure heave frequencyfhis 9.5 Hz, and the pure pitch frequencyfθis 11.4 Hz, 9.5 Hz, and 8.1 Hz, respectively.The vibration state, maximum power and average work for each dimensionless inertia radius are shown in Table 9.The response results of the three models at a flow velocity of 5 m/s are shown in Figures 31 and 32.It can be seen that whenr=0.45 and 0.634, the hydrofoil no longer maintains constant amplitude vibration and exhibits an attenuation trend.The power of fluid work and the energy input for the elastic system significantly decreases with the vibration of the hydrofoil, and whenr=0.45, the energy and motion of the system attenuate more rapidly, and the vibration enters a convergent state.

Table 9 Maximum power and average work for each dimensionless inertia radius

Figure 31 Vibration and hydrodynamic force time history results at different inertia radii (5 m/s)

The response results of three inertial radius hydrofoils when the flow velocity increases to 6 m/s are shown in Figures 33 and 34.

At a flow velocity of 6 m/s, the original model (r=0.542) was already in a flutter state, while the heave and pitch vibrations of the hydrofoil models withr=0.45 andr=0.634 continued to attenuate.The energy input by the fluid to the elastic system gradually decreased, and the vibration remained in a convergent state, but the amplitude signifi‐cantly increased compared to 5 m/s.

Continue to increase the flow velocity to 7 m/s, and the response results are shown in Figures 35 and 36.

Figure 32 Hydrodynamic power and work history results at different inertia radii (5 m/s)

Figure 33 Vibration and hydrodynamic force time history results at different inertia radii (6 m/s)

Figure 34 Hydrodynamic power and work history results at different inertia radii (6 m/s)

Figure 35 Vibration and hydrodynamic force time history results at different inertia radii (7 m/s)

Figure 36 Hydrodynamic power and work history results at different inertia radii (7 m/s)

At this time, the heave and pitch vibrations of the two models no longer attenuate, and the amplitude increases faster atr=0.634 than atr=0.45.However, due to the larger inertia radius, the vibration period is longer than atr=0.45,and the energy input by the fluid for the two models gradu‐ally increases.The hydrofoil enters a flutter state, and the critical velocity increases from 5 m/s to between 6 m/s and 7 m/s.It can be seen that the radius of inertia has a signifi‐cant impact on the amplitude and frequency of both heave and pitch motions, and can change the critical velocity.

4 Conclusion

Based on computational fluid dynamics and structural finite element methods, a co-simulation framework for the flow-induced vibration of hydrofoil in uniform flow is established and numerical simulation research is carried out.The conclusions drawn from this work can be summa‐rized as follows:

1) The flow velocity has a significant impact on the am‐plitude of the hydrofoil.The heave and pitch motions at low flow velocity are mainly affected by the second-order natural characteristics of the elastic system, and as the flow velocity increases, vibration begins to be dominated by a single frequency.When the input value of the fluid energy exceeds the dissipation value, the hydrofoil will enter the flutter critical state without amplitude attenuation.

2) The attack angle of the hydrofoil has a significant impact on the vibration amplitude, but has no apparent effect on vibration frequency and critical velocity.

3) Adjusting the natural frequency ratio by reducing the vertical stiffness or increasing the pitch stiffness can move the vibration from a critical state to a convergent state.To increase the critical velocity of the flutter, the natural fre‐quencies of heave and pitch can be staggered by adjusting the inertia radius.

Different geometric shapes will alter the excitation force characteristics of the flow field, thereby affecting the vibration of the hydrofoil.This work mainly calculates and analyzes the vibration characteristics of hydrofoils at different flow velocities and angles of attack, so it is possi‐ble to explore the impact of different trailing edge shapes on hydrofoil vibration in the future.

FundingThis work is supported by the National Natural Science Foundation of China (Grant No.52001043), the Chinese Academy of Sciences Youth Innovation Promotion Association (Grant No.2020205), the Fundamental Research Funds for the Central Universities (Grant No.DUT22GF202 and DUT20TD108) and Liaoning Revitalization Talents Program (Grant No.XLYC1908027).

Competing interestThe authors have no competing interests to de‐clare that are relevant to the content of this article.