The effect of dynamic twisting on the flow field and the unsteady forces of a heaving flat plate

Many marine animals can dynamically twist their pectoral fins while swimming. The effects of such dynamic twisting on the unsteady forces on the fin and its surrounding flow field are yet to be understood in detail. In this paper, a flat plate executing a heaving maneuver is subjected to a similar dynamic twisting. In particular, the effects of the direction of twist, non-dimensional heaving amplitude, and reduced frequency are studied using a force sensor and particle image velocimetry (PIV) measurements. Two reduced frequencies, k=0.105 , and 0.209 , and two twisting modes are investigated. In the first twisting mode, the plate is twisted in the direction of the heave (forward-twist), and in the second mode, the plate is twisted opposite to the direction of the heave (backward-twist). Force sensor measurements show that the forward-twist recovers some of the lift that is usually lost during the upstroke of flapping locomotion. Additionally, the forward-twist maintains a near-constant lift coefficient during the transition between downstroke and upstroke, suggesting a more stable form of locomotion. PIV results show that forward-twist limits circulation and leading-edge vortex growth during the downstroke, keeping Cd≈0 at the cost of the reduced lift. By contrast, backward-twist increases the circulation during the downstroke, resulting in large increases in both lift and drag coefficients. Force sensor data also showed that this effect on the lift is reversed during the upstroke, where the backward-twist causes a negative lift. The effects of each twisting mode are mainly caused by the changes in the shear layer velocity that occur as a result of twisting about the spanwise axis along the mid-chord. The twisting performed by forward-twist reduces the effective angle of attack through the upstroke and downstroke, resulting in a reduced shear layer velocity and lower circulation. The twisting performed by backward-twist does the exact opposite, increasing the effective angle of attack through the upstroke and downstroke and consequently increasing the shear layer velocity and circulation.


Introduction
Remotely operated aerial and marine vehicles mostly operate with rigid wings, which are often limited in terms of maneuverability, speed, and efficiency, unlike natural flyers and swimmers, who have compliant wings and fins. These natural flyers and swimmers rely on alternative methods of flight and propulsion, which can be found in the unsteady flapping motions of birds, bats, insects, and rays [1]. These creatures are proficient in controlling the unsteady flow field, solely through morphing and flapping, in a manner that allows them to maintain a high level of agility [2,3]. Recent studies into unsteady flapping locomotion have been primarily concerned with the heaving and pitching motion of rigid and passively flexible foils [4][5][6]. In comparison, there is a rarity of papers that focus on dynamic shape change of lifting and control surfaces.
Dynamic shape changes are observed in many instances of natural flight. In a study monitoring the bending of the propulsive surfaces of various natural animals, it was found that animals bend their propulsor's span at about two-thirds of the total span [7]. Butterflies produce more lift during forward flight with the help of time-varying wing twists [8]. Manta rays utilize the flexibility of their pectoral fin to execute turning maneuvers at a higher rate than other batoid rays [9]. The dexterity of natural flight often surpasses that of engineering flight vehicles in highly unsteady flows. In many such instances, engineering flights spin out of control, whereas birds and insects maintain stability [10,11]. For example, during rapid pitching of the wing, the flow becomes highly unsteady, affecting aerodynamic performances. During such rapid pitching, massive separation of the flow leads to the formation of a leading-edge vortex (LEV) [12,13]. In general, such vortex formation is detrimental to stable aerodynamic performance as it leads to severe lift oscillations. However, birds use these vortices during a controlled landing-the famous perching maneuver-to augment lift [14].
Flapping flight has various mechanisms of generating additional unsteady lift, such as wake capture, added mass, clap and fling motion, and LEV development [15]. The use of these mechanisms is not mutually exclusive, and natural flyers employ situation-specific combinations of them to enhance their flight. Of these methods, the ability to control the development of the LEV is one of the strongest. Research has shown that the formation, size, and strength of an LEV increase circulation and consequently increase the lift force [16]. Flapping also causes vortices to detach from the wing and become entrained in the wake, boosting circulation even further [16]. Several researchers used rigid wings to investigate the formation of LEVs at low Reynolds number flows typical of natural flights. LEVs form during a pure heaving motion or during a pitch-heave combined kinematics which is often used to emulate flapping wings [17][18][19][20][21][22][23]. The strength, growth rate of circulation, and position of an LEV are the main parameters influencing the aerodynamics around the wing when such vortices form [24]. In the case of the pitching-plunging wing, the LEV formation is controlled by several nondimensional parameters, such as the Strouhal number, the reduced frequency, the pitch-plunge phase angle, etc [25]. For a finite wing, the wing's aspect ratio adds to the parameter space [23]. Among these nondimensional parameters, reduced frequency is of particular interest. It characterizes the oscillatory flapping cycle as a measure of the flapping velocity relative to the flight speed. This relation makes it useful as a measure of the degree of unsteadiness in the flow, where higher values indicate higher unsteadiness. The reduced frequency is defined as: Here ω = 2πf is the circular frequency of flapping, b = c/2 is the semi-chord of the wing, and V = U ∞ is the freestream velocity. The Strouhal number, St, relates the heaving frequency and the forward flight speed and is defined as: Here f is the flapping frequency, and h is the heaving amplitude. Taylor et al studied this nondimensional parameter for various natural flyers and swimmers, discovering that the optimal St for propulsive efficiency lies within the range 0.2 < St < 0.4 due to increased LEV formation and growth [26].
The lift response of a maneuvering wing is further modulated in the case of a flexible wing. Several researchers have shown that the lift response improves due to flexibility during such LEV formation. This improvement was attributed to the increased stability of LEV [27,28]. Kim and Gharib showed that a flexible plate oriented at a 90 • angle of attack experienced a lower initial peak force, although later, the same plate generated higher forces compared to a rigid plate due to slower development of vortices [29]. Spanwise tip flexibility has been shown to lead to better flight performance for autorotating wings, which is attributed to flutter-induced drag reduction [30]. For flapping wings, twisting and bending in the flapping plane influence the lift only at specific points of the flapping cycle, whereas outof-plane bending has a more global effect [31]. For some birds with optimal aspect ratio wings, the flapping frequency is near the first spanwise natural frequency of the wing, which suggests that they utilize such resonance to generate thrust [32]. In a water tunnel study with flexible heaving wings, Heathcote and Gursul showed higher efficiencies for St between 0.2 and 0.4 [33].
More recent studies have shown that St is not a universal parameter for the evaluation of optimal flapping motion. Baik et al found that the effective angle of attack, α eff , and the reduced frequency, k, determines the flow evolution, while St determines the aerodynamic forces [17]. Gursul and Cleaver argued that St is only the dominant parameter for attached flow at small effective angles of attack; however, for detached flow at high effective angles of attack, the reduced frequency becomes dominant [4]. The effective angle of attack is defined by Gursul as shown in equation (3): Here α is the geometric pitch angle, U pl is the heaving velocity, and U ∞ is the freestream velocity. It is important to note that in Gursul and Cleaver's studies, the geometric angle of attack is kept constant. In another study, Chiereghin et al found that the effective angle of attack correlates with increased lift rather than St [4]. They also found that maximum lift in an oscillatory heaving cycle is determined by circulatory lift [4]. The circulatory lift is found to be highest near the maximum effective angle of attack, regardless of the presence of an LEV [4].
These previously mentioned studies highlighted the importance of effective angle of attack and circulation on flapping flight. As demonstrated by Baik and Bernal, pitching from the root is a useful mechanical method of modifying the effective angle of attack [34]; however, very few natural flyers have the joint mobility required to perform such maneuvers. Most flyers, such as birds and bats, employ an alternate method of controlling the effective angle of attack via active morphing [2,[35][36][37]. Since these flyers have anatomical limitations in the form of bones and ligaments, they perform twisting and bending motions along the wingspan instead. A 2014 study of the aerodynamic mechanisms of bat flight confirmed that the spanwise twisting of a bat wing performs the same function as the supination and pronation motions in insect flight [35].
This present paper addresses this issue by studying how a dynamically twisting foil affects the fluid-structure interactions of unsteady heaving locomotion. To accomplish this, a custom-built experimental setup was used to heave and morph a flexible hydrofoil simultaneously through a large water tank over a range of conditions. The hydrodynamic forces on the foil were measured with a 6-axis force sensor mounted near the root of the foil. particle image velocimetry (PIV) measurements were carried out at different stages of the heaving cycle to study the development of vorticity over the suction surface of the heaving hydrofoil. This research focuses particularly on the effects of spanwise twisting on the aerodynamic forces, LEV development, and circulation in the flow. Details on the experimental method and parameter space are presented in section 2, followed by the details of the analytical model in section 3 results and discussion in section 4. Finally, in section 5, we discuss the conclusion and the outlook of this research.

Plate and the twisting mechanism
A flat plate with a chord of 10 cm, a span of 30 cm (aspect ratio = 3), and a thickness of 6 mm, was used in this study ( figure 1(a)). It had a rounded leading edge and a tapered trailing edge. The flexible plate was manufactured from Agilus Black plastic (shore hardness 70 A) using a polyjet 3D printer. It had two internal cavities for housing a set of thin metal rods that were used for twisting the plate. These rods were curved at one end and routed internally through the plate such that the curved portion begins at a flexion ratio of 0.6. The flexion ratio is defined as the ratio between the span length of the beginning of curvature to the total span. Twisting is accomplished by rotating these metal rods via a pair of servo motors (SAVOX SA-1230SG) controlled by an Arduino microcontroller.
Two twisting modes were selected for testing and are denoted as follows: forward-twist and backwardtwist. Forward-twist is characterized as spanwise twisting in the direction of heaving motion, reducing the effective angle of attack throughout the cycle. Backward-twist is characterized as spanwise twisting in the direction opposite to the heaving motion, increasing the effective angle of attack throughout the cycle. In both cases, the tip of the plate was twisted to a maximum angle of 45 • . A diagram displaying the different twisting modes is shown in figure 2.

Towing tank and plate kinematics
Experiments were conducted inside a water tank with dimensions 6 × 1 × 1.5 m (length × width × height) and with 1.25 cm-thick walls ( figure 1(b)). The plate model was towed inside the water with the help of a linear traverse (HPLA80, Parker Motion) controlled by an AC Servo Controller (Applied Motion SV200) via its associated Applied Motion Q-Programmer software. In all cases, the freestream velocity is kept constant at U ∞ = 0.1 m s −1 . Assuming the kinematic viscosity of water is ν = 1 × 10 6 m 2 s −1 , and a chord of 10 cm, this resulted in a constant chord-based Reynold's number of Re = 10, 000. A skim plate was attached just above the root of the plate to suppress the formation of surface waves.
To implement the heaving motion, the plate was mounted on a shaft that was connected to a rotary stage motor (Zaber RSW60A-T3) through its midchord axis for pitch angle control. This motor was used to set the plate at a constant geometric angle of attack of α = 15 • . This rotary stage motor was connected in turn to a linear stage motor (Zaber LSQ150B-T3) via a connecting metal plate to provide the heaving motion. The heaving motion is given by equation (4), where A 0 is heaving amplitude in meters, f is the frequency of oscillation in cycles per second, and t is time in seconds. Heaving amplitude h is nondimensionalized by the chord length, c, and time t is nondimensionalized by the cycle frequency, as shown in equations (5) and (6). Here T is the total elapsed time of one cycle: Two separate heaving amplitudes, namely h * = 0.5 and 1.5 were implemented in this work. For each heaving amplitude, two separate reduced frequencies, namely k = 0.105 and 0.209, were tested. The reduced frequency k is defined as k = πfc U∞ .  The synchronization of the mechanical components to generate the prescribed kinematics is accomplished via an external trigger signal. This trigger is activated by a proximity sensor mounted on the linear traverse such that triggering occurs when a small metal piece on the carriage passes over the sensor. This trigger signal is then sent to Arduino and Zaber Console to begin the morphing and plunging cycles, respectively.

Twist profile measurements
Measurements of the twist profile were performed by placing markers near the leading and trailing edge of the plate surface at multiple spanwise locations, then tracking their position throughout the twisting motion. The twisting motion was recorded using two high speed cameras (PCO Edge 5.5, resolution 1920 × 1080 pixels) at a frame rate of 200 Hz. The recorded images were then post-processed using DLTdv, a MATLAB-based digitizing tool, to track the marker locations in 3D space. Further postprocessing was conducted in Python to compute the twist angle evolution at each spanwise location. The computed twist profiles are shown in figure 3.

Force sensor
A force sensor was mounted on the shaft that holds the plate such that its central axis is colinear with the axis of rotation. A 6-axis force/torque sensor (Nano 25, ATI) was used to sample the instantaneous forces. Force sensor data was acquired using a LabVIEW program.
All force sensor data was taken with a sampling frequency of 1000 Hz, with the initiation of acquisition occurring immediately upon triggering. The trigger signal used is the same one that is sent by the proximity sensor mentioned previously to keep all systems synchronized. The final force values were obtained after averaging over five runs, which were then low-pass filtered with a cutoff frequency of 7 Hz.
Due to the nature of the setup, the local coordinate system of the force sensor was not aligned with the local coordinate system of the plate. Since the angle of attack was kept at 15 • , the axial and normal force outputs of the force sensor (F x and F y ) did not directly represent the lift and drag forces. Additionally, the effective angle of attack, α eff , was constantly changing due to heaving motion. As such, it was necessary to first determine the distribution of the effective angle of attack throughout the plunging cycle as defined in equation (7), whereḣ is the heaving velocity: The heaving velocity profiles were obtained directly from the linear stage motor via the Zaber Console software. The instantaneous lift and drag forces could then be determined using the trigonometric relations shown in equations (8) and (9), where F x and F y are the axial and normal force outputs, respectively: The lift and drag measurements are nondimensionalized as C l and, as C d shown in equations (10) and (11). In these equations, F L is the lift force, F D is the drag force, ρ is the density, U ∞ is the freestream velocity, s is the span, and c is the chord:

Particle image velocimetry (PIV)
Planar PIV experiments were conducted at two spanwise locations (50% and 70% span) to measure the velocity field of the heaving plate (figure 4). The flow field was illuminated by a laser plane using a continuous wave and diode-pumped solidstate laser (CNI MGL-F-532, DragonLaser, China). The laser beam was passed through a set of cylindrical lenses (Thorlabs) that created a 2D laser sheet. Silver coated-glass spheres of diameter 100 µm (Conduct-O-Fil, Potter Industries) were used as tracer particles. Time-resolved images were recorded using a high-speed camera (PCO 12000HS, resolution 2048 × 2012 pixels) at a frame rate of 200 Hz and an exposure time of 5 ms. The field of view (FOV) of the camera at 50% span corresponded to a region of size 24.3 cm by 19.5 cm, which resulted in a spatial resolution of 0.19 mm/pixel. At 70% span, the FOV corresponded to a region of size 27.2 cm by 21.7 cm, which resulted in a spatial resolution of 0.21 mm/pixel. The camera was mounted to the linear stage motor, such that it moved with the traverse. The recording sequence was initiated via an Arduino controller and the previously mentioned trigger signal from the proximity sensor. The PIV images were post-processed using PIVlab, a MATLAB-based software. An FFT-based cross-correlation method was used to correlate two consecutive images. A multi-pass scheme comprising 64 × 64 pixels followed by 32 × 32 pixels was used for computing the flow vectors. This resulted in a flowfield resolution of 1 vector every 0.62 (x) by 0.65 (y) cm in the flow field at a 50% span. At 70% span the flow-field resolution is 1 vector every 0.69 (x) by 0.72 (y) cm. Further, post-processing was carried out in MATLAB to compute vorticity and circulation.
PIV measurements were taken on time interval 0.1 < t * < 0.3 during the downstroke of the heaving cycle. Figure 5 shows this interval in relation to the heaving position and the heaving velocity of the hydrofoil. For PIV experiments, vorticity ω is nondimensionalized by chord and freestream velocity, as shown in equation (12): Circulation is also computed for analysis. The vortex identification method used is referred to as ′ ′ Γ 1 ′ ′ and ′ ′ Γ 2 ′ ′ criteria, where Γ 1 identifies the vortex core location, and Γ 2 identifies the vortex boundary [17,38]. The two criteria are shown in equations (13) and (14), respectively: Here p is any point in the flow field, x is the position vector, u is the velocity vector, n is the unit vector in the z-direction, N is the total number of points in the control region, and u p is the average velocity of the same control volume. A vortex core is identified when |Γ 1 | > 0.9, and its boundary is characterized by values of |Γ 2 | > 2/π. Once the vortex boundary has been determined, the circulation is computed by executing a summation of the z-component vorticity enclosed by the vortex boundary. The circulation is then nondimensionalized using freestream velocity U ∞ and chord length c, as shown in equation (15):

Analytical model
Theodorsen's classical unsteady model is used to compare the unsteady forces of the heaving and twisting flat plate [39]. This model assumes a small angle of attack, attached flow, and a planar wake. Despite this, the model has been successfully used for airfoils undergoing sinusoidal plunging and pitching motion [40]. Unsteady lift per unit span is expressed as the summation of circulatory lift, C l,C , and noncirculatory lift, C l,NC : The circulatory and non-circulatory lift terms are expressed in equations (17) and (18), respectively: Here ρ is the density of water, b is the halfchord length (c/2), and U ∞ is the freestream velocity. Pitch angle, pitch velocity, and pitch acceleration are denoted as α,α, and .. α, respectively. Plunge velocity and plunge acceleration are denoted asḣ and .. h, respectively. C (k) is the Theodorsen function evaluated at the reduced frequency, k, and is expressed as: where H (2) 1 (k) and H (2) 0 (k) denote Hankel functions. Twisting is modeled as a spanwise change in the pitch angle. This is expressed in the lift terms as a spanwise integration of pitch angle, pitch velocity, and pitch acceleration, as shown in equations (20) and (21), respectively, where s * is the nondimensional span (s/s tot ): ds * (20) The Theodorsen-based model was computed by taking the first six planes (1.00 ⩽ s * ⩽ 0.85) of the twist angle measurements shown in figure 3 and performing a 1000-point linear interpolation in the time domain. Linear interpolation was chosen to match the linear slopes observed in figure 3. The first six planes where chosen because these are spaced evenly along the span (ds * = 0.03), whereas the last two are not. For each time step, a linear curve fit was then performed along the twisted spanwise region (1.00 ⩽ s * < 0.40). This fitted twisting region was then combined with the untwisted region (0.40 < s * ⩽ 0.00), where twist angle θ ≡ 0 • , to obtain a 100-plane spanwise twist profile of the full wing. The 2D Theodorsen solution for each of these planes was then integrated across the full span to obtain the full analytical solution.

Results and discussion
In this section, force sensor measurements are first used to investigate the effect of dynamic twisting on the unsteady forces of the plate. Next, PIV measurements are used to study the development of the flow field of the plate.

Unsteady lift forces
The variations in the coefficient of lift for k = 0.105, h * = 0.5 are shown in figure 6. The solid lines shown in figures 6-9 represent the experimentally measured force values, whereas the dotted lines represent estimated values obtained from the analytical model. The differences exhibited between the theoretical and experimental values can be attributed to various physical phenomena that are not captured by the analytical model. Firstly, the analytical model assumes attached flow along the surface of the plate, whereas flow separation is observed in experiments due to plunging and a high angle of attack. Next, the model is based on potential flow, which does not take into account the formation, growth, and shedding of the LEV present in experiments. Finally, since the analytical model is constructed as a spanwise integration of 2D planes, 3D effects-such as spanwise vorticity transport-are not captured.
The rigid case shown in figures 6-9, corresponds to the sinusoidal heaving motion of the plate without any twisting. The lift during both the forward-twist and backward-twist cases significantly deviates from the rigid case. In the forward twist case, the lift curve shows a downward peak during the downstroke around t * = 0.27 for each case. This effect is reversed during the upstroke, creating an upward lift peak around t * = 0.73. This upstroke effect is of particular interest as the lift is typically lost during this part of the flapping cycle, as evidenced by the rigid mode. For reduced frequency k = 0.105 , the peak-to-peak amplitude of the transient lift coefficient in the forward twist case is 75.3% and 4.6% greater than the rigid case for heaving amplitudes h * = 0.5 ( figure 6) and h * = 1.5 ( figure 7), respectively. For reduced frequency k = 0.209 , the peak-to-peak amplitude of the transient lift coefficient in the forward twist case is 32.6% greater and 55.5% lower than the rigid case for heaving amplitudes h * = 0.5 (figure 8), and h * = 1.5 ( figure 9), respectively. The forward twist case also maintains an entirely positive lift throughout the plunging cycle.
This suggests a more stable method of flapping locomotion in the forward-twist case, where the vehicle does not 'bounce' up and down as much as it moves forward. This stabilizing effect is beneficial as it makes the path of locomotion more linear and, thus, easier to control. The effect of backward-twist with respect to the rigid case is opposite to that of the forward-twist case. This twisting mode has an amplifying effect on the transient coefficient of lift, as it displays the same sinusoid pattern but with more extreme upward and downward peaks in the lift profile. In each case, the backward-twist increases the peak transient lift coefficient during the downstroke.  Figures 6-9 also show the effect of reduced frequency on the transient lift coefficient for each twisting mode. For both the rigid mode and backward-twist, we find similar sinusoidal trends, which proves that backward-twist amplifies the forces experienced during the rigid case. For the rigid case and the backwardtwist case, increasing the reduced frequency to k = 0.209 increased the peak-to-peak magnitude of the lift coefficient profile by 99.0% and 25.7% for h * = 0.5 (figure 8), compared to k = 0.105 ( figure 6). For  h * = 1.5, these increments were 109.8% and 50.2% respectively (figures 9 and 7). For a periodic plunging wing, such an increase in the lift coefficient due to increasing reduced frequency has also been reported by Gursul and Cleaver [41]. However, our results clearly show that twisting the wing backward can augment this lift much further from the rigid case during the downstroke part of the cycle.

Effect of reduced frequency on lift
The effect of increasing the reduced frequency is less prominent during the forward-twist case compared to the rigid and backward-twist case, especially in the lower plunging amplitude case, i.e. h * = 0.5 (figures 6 and 8). This clearly points to the force-canceling effect of forward-twist and shows that it is possible to reduce the lift excursions during an unsteady maneuver simply with dynamic  twisting. In the higher amplitude case of h * = 1.5, the effect of reduced frequency is unique during forward-twist. For this case, the lower reduced frequency k = 0.105 has a higher lift coefficient peak of C l = 1.5 at t * = 0.75 (figure 7) during the upstroke than the higher reduced frequency, k = 0.209 ( figure 9). This difference shows that the effect of lift cancellation by forward twist is more prominent in the higher plunging amplitude case accompanied by the higher reduced frequency. This observation warrants further exploration of forward-twist in unsteady flow, as it could be used to negate the effect of unsteady forces in such flows suggesting a more efficient method of flow control leading to increased flight stability.

Effect of heaving amplitude on lift
For most of the cases shown in figures 6 to 9, the nondimensional heaving amplitude (h * ) has similar effects on the transient coefficient of lift as the reduced frequency. In both the rigid mode and backward-twist, the higher nondimensional heaving amplitude, h * = 1.5, increases the peak-to-peak amplitude of the lift coefficient from the h * = 0.5 cases. When k = 0.105, increasing the heaving amplitude causes a 32.7% increase in the peak amplitude of the lift coefficient in the rigid case. These values were 12.4% and 17.0% in the backward-twist and forwardtwist cases, respectively (figures 6 and 7). Once again, a forward-twist shows a negating effect compared to the other two cases, with the effects of the nondimensional heaving amplitude being more prominent in the case of higher reduced frequency k = 0.209. At this reduced frequency, the lower nondimensional heaving amplitude case h * = 0.5 in the forward twist case displays a higher lift coefficient, C l = 1.65 (at t * = 0.8), during the upstroke compared to C l = 1.22 at the same t * in the case of h * = 1.5. Expanding upon the reduced frequency analysis of this case, lower heaving amplitude combined with low reduced frequency is observed to mitigate the negative impact on the lift coefficient during the upstroke, potentially due to vortex development. Due to the previously mentioned nature of heaving amplitude, this reduction in upstroke duration also applies to the downstroke; where its benefits of lift generation are also mitigated and completely nullified at the peak of the downstroke around t * = 0.23.

Discussion on analytical model
The analytical model for the rigid and backward-twist cases exhibits a sinusoidal lift profile, with maximum and minimum values observed at t * = 0.25 and t * = 0.75, respectively. By contrast, the forward-twist case is observed to have reversed peaks, with the minimum lift value observed at t * = 0.25, and the maximum lift value observed at t * = 0.75.
A similar lift profile is seen in the experimental data for these cases, but there are some key differences. For rigid and backward-twist, the maximum experimental lift is larger than the model prediction for all reduced frequencies and heaving amplitudes. The minimum experimental lift is also lower than the model prediction through these parameters. For k = 0.209 this discrepancy is more pronounced at both heaving amplitudes, as the high degree of flow separation in these cases is not captured by the model (figures 8 and 9).
By contrast, for both reduced frequencies and heave amplitudes of forward-twist, the model predicts a higher lift magnitude through the downstroke than the experimental data. The forwardtwist motion interrupts the formation of an LEV, cancelling its lift-enhancing effects. The leading-edge shear layer velocity that would otherwise feed the LEV gets reduced transported along the twisted part of the plate. This reduces the circulation in the region above the top surface of the plate, reducing the lift even further. This effect can also be seen in the upstroke of the forward-twist case, where reduced circulation on the lower surface results in lift enhancement.
For rigid and backward-twist at k = 0.105, the maximum experimental lift occurs earlier during the downstroke (t * < 0.2) than the model prediction for both heaving amplitudes (figures 6 and 7). This shift can be attributed to the influence of the LEV, which becomes more pronounced as the reduced frequency is lowered. The formation and growth of the LEV during the downstroke imparts additional lift, causing the peak to occur sooner while simultaneously increasing its magnitude. For forward-twist and backward-twist, this shift is reversed in the upstroke. Maximum lift for forward-twist occurs later during the upstroke than predicted (0.76 < t * < 0.82) for both reduced frequencies and heaving amplitudes. Similarly, the minimum lift for backward-twist also occurs later in the same time range.

Unsteady drag forces
The effect of twisting on the transient drag coefficient is explored in figure 10, which shows the evolution of the transient drag coefficient for each twisting mode as a function of non-dimensional time, t * . The rigid mode in figure 10 once again exhibits a sinusoidal trend corresponding to the heaving cycle, maintaining a transient drag coefficient C d < 1 in all cases except for the one in the case of k = 0.209 and h * = 1.5. By comparison, forward-twist displays a very low, almost constant trend during the downstroke and a drag peak during the upstroke. Its effect on the downstroke can be attributed to the change in the angle of attack along the span caused by the twisting. During the downstroke, the twisting action of forward twist lowers the spanwise distribution of the angle of attack from α = 15 • towards α = 0 • , such that the chord of the plate is more aligned with the direction of freestream velocity. This causes a reduction in the effective area facing the flow, thus reducing the drag coefficient such that C d ≈ 0. By contrast, backward-twist increases the drag coefficient during the upstroke. The same explanation given for the downstroke of the forward-twist also applies to the upstroke; however, during the upstroke, the spanwise distribution of the angle of attack is increased such that the chord of the hydrofoil is less parallel with the direction of freestream velocity. This change in the spanwise angle of attack distribution increases the effective area, which leads to increased drag and the shear layer velocity, thus increasing the circulation and transient drag coefficient. The backward-twist cases show drag coefficient peaks during both the downstroke and upstroke parts of the heaving cycle. This effect is much more pronounced during the downstroke, where the backward-twist exhibits a transient drag coefficient C d ⩾ 1 in all cases. Thus, while backward-twist shows great improvement in the lift coefficient during this part of the cycle, it comes at the cost of a much higher drag coefficient when compared to the other modes. The drag coefficient peak in the upstroke is smaller than that of the downstroke by a factor of two or more in all cases. Additionally, this upstroke peak is smaller than that of the forward-twist, despite an overall increase compared to the rigid mode. This upstroke effect of backward-twist is similar to the downstroke effect of forward-twist, albeit to a lesser degree, where the spanwise angle of attack distribution approaches zero during this part of the cycle.
Comparing the lift coefficient plots (figures 6 to 9) to the drag coefficient plot (figure 10), it can be seen that the increased lift coefficient effects of the twisting modes come at a trade-off with their drag coefficient effects. Further exploration is required to understand how these twisting modes affect the overall hydrodynamic efficiency of the heaving cycle.

Effect of reduced frequency on drag
The effect of reduced frequency on the transient drag coefficient is now explored. Figure 10 shows the transient drag coefficient for both reduced frequencies as a function of cycle-averaged time for each twisting mode. Figure 10 shows similar trends between the reduced frequencies for the rigid and backward-twist modes, where the lower reduced frequency case, k = 0.105, follows the same trend as the high reduced frequency case, k = 0.209, but at lower overall values of the transient drag coefficient. Once again, the effect of reduced frequency experienced by the forward-twist is unique, displaying a higher drag coefficient peak during the upstroke of the case with high amplitude, h * = 1.5, and high reduced frequency k = 0.209. This effect can also be seen in the high amplitude case of the rigid mode. Overall, the lower reduced frequency for this particular case of forward-twist exhibits a higher transient drag coefficient, again highlighting the trade-off between increased lift coefficient and increased drag coefficient.
Comparing the lift and drag results, the data suggests that it may be disadvantageous to maintain the same reduced frequency throughout the cycle. Instead, it may be more favorable to have a higher reduced frequency during the downstroke and a lower reduced frequency during the upstroke. This effect suggests that vortex development at lower reduced frequencies may play a key role in recovering some of the lost lift experienced during the upstroke. Further investigation into this phenomenon is needed to understand the mechanisms behind this effect.

Effect of heaving amplitude on drag
The effect of nondimensional heaving amplitude on the transient drag coefficient is now explored. Figure 10 shows the transient drag coefficient for both amplitudes as a function of cycle-averaged time for Figure 11. Evolution of the vorticity field during the downstroke, and at the 70% span location of the plate in the rigid case, with k = 0.105, and h * = 1.5. each twisting mode. Figure 10 shows that decreasing the nondimensional heaving amplitude also has a dampening effect on the drag coefficient distribution for the rigid mode and backward-twist, akin to its effect on the lift coefficient. Forward-twist again differs in its downstroke and upstroke effects. During the downstroke, the lower amplitude, h * = 0.5, dampens the drag coefficient similar to the other twisting modes. By contrast, at k = 0.105 the drag coefficient is very similar between heaving amplitudes, and at k = 0.209 the drag coefficient is increased. From this analysis of reduced frequency, and the twisting modes, it is evident that forward-twist has significant effects on the hydrodynamics of the heaving hydrofoil.

PIV analysis 4.3.1. Effect of twisting on the vorticity field
To evaluate the effect of morphing, we first analyze the evolution of the vorticity field at the 70% span location of the plate during the rigid, forward-twist, and backward-twist cases. The span length of 70% was chosen as the magnitude of twisting is higher at this location than at the 50% span, and as such, the effects of twisting on the vortex development within the flow field are more visible. Results are shown for 0.15 < t * < 0.30, with a reduced frequency, k = 0.105, and nondimensional heaving amplitude h * = 1.5. It is to be noted that the time range 0.15 < t * < 0.30 occurs during the downstroke part of the heaving cycle. The TEV was not captured during PIV measurements, as we focused only on the LEV development. Figure 11 shows the vorticity field at the 70% span of the rigid case with k = 0.105, andh * = 1.5. In figure 11, it is observed that the LEV is initially attached and increasing in size during0.15 < t * < 0.21. As the downstroke continues into the second interval, 0.24 < t * < 0.30, the heaving velocity keeps increasing until t * = 0.25. Hence, the effective angle of attack and consequently the shear layer velocity also increased until t * = 0.25, resulting in rapid growth and shedding of the LEV. Once shed, the LEV passes over the plate along the streamwise direction towards the wake.
The effects of forward-twist during the downstroke (0.15 < t * < 0.30) with k = 0.105, and h * = 1.5 are displayed in figure 12. It is evident that due to the forward-twist, the 70% span location almost reaches a geometric angle of attack, α = 0 • in the range 0.15 < t * < 0.24 and is reduced further into the negative range, α < 0 • , during 0.24 < t * < 0.30. Since this reduction in geometric angle of attack is occurring during the downstroke, the effective angle of attack α eff is also reduced, resulting slower shear layer velocity. This limits the growth of the LEV, allowing it to stay attached for a longer period of time. Figure 12 also shows how the reduced α eff prevents the vortices from being shed. Instead, the LEV stretches along the chord of the hydrofoil, sliding directly on the suction surface until it reaches the trailing edge.   This combination of reduced shear layer velocity and reduced LEV development severely limits circulation on the suction side of the plate. This limited circulation also results in less circulatory forces, which is reflected in the force data, where the forward-twist demonstrates low lift and drag coefficients during the downstroke.
The effect of backward-twist at the 70% span of the plate is displayed in figure 13. Since at this span location, the chord of the plate is being twisted upwards, the local effective angle of attack keeps increasing in this case during 0.15 < t * < 0.30. When combined with the heaving motion, this increased α causes an increase in shear layer velocity. This results in a larger LEV and vorticity region above the suction surface of the hydrofoil. Figure 13 also shows greater separation between the vortices in the second interval, suggesting that the increased spanwise twist and heaving velocity cause the LEV to be shed more quickly. However, some of this shed vorticity is trapped in the recirculation region above the suction surface of the plate. This recapture of shed leading-edge vorticity combined with a larger recirculation region results in highly increased values of lift and drag coefficients during the downstroke of the heaving cycle.
Next, the circulations inside the LEV were calculated during the first PIV interval for each twisting mode. The effects of nondimensional heaving amplitude on the circulation are studied for both 50% and 70% span lengths. Figures 14 and 15 display the transient nondimensional circulation of each twisting for nondimensional heaving amplitudes at both span lengths, for k = 0.105 and k = 0.209 respectively.
The circulation figures show us how the degree of spanwise twisting affects the circulation. The spanwise twisting movement is larger at 70% span, closer to the tip of the wing of the hydrofoil than at 50% span. As such, at 70% span the circulation profiles of the forward-twist and backward-twist modes deviate further from the rigid mode than at 50%. This effect is most noticeable for the h * = 0.5 case of forwardtwist. For this case, the 50% span measurement shows that the rigid mode and forward-twist share very similar circulation profiles for 0.14 < t * < 0.19. By contrast, at 70% span, forward-twist has a much weaker circulation than the rigid mode for the same time interval.
The circulation profiles shown in figures 14 and 15 are consistent with the observations made in the vorticity subsection of the PIV data analysis. Forward-twist exhibits decreased circulation across all demonstrated parameters, confirming the effects of reduced shear layer velocity that arise as a result of spanwise forward twisting. Backward-twist exhibits the opposite behavior, showing increased circulation among the demonstrated parameters. The effect of twisting is more pronounced in the higher amplitude h * = 1.5 cases, most notably for backwardtwist where circulation increases drastically from 0.19 < t * < 0.23. By contrast, the effect of twisting is limited at higher reduced frequency. This is seen most clearly in the forward twist-case, where the circulation deviates less from the rigid case at k = 0.209 than at k = 0.105. Similar results were also obtained by Rival et al (2009) for a plunging airfoil [42]. They showed that at higher reduced frequency cases, the circulation development of the LEV is limited by the period of oscillations.

Conclusion
Two spanwise twist morphing modes were investigated at two nondimensional heaving amplitudes and two reduced frequencies. The effect of twisting, heaving amplitude and frequency on the lift and drag coefficients were investigated. PIV measurements were taken during the downstroke of the heaving cycle to explore the fluid-structure interactions exhibited by the two twisting modes.
Heaving amplitude and reduced frequency had similar effects for both the rigid case and the backward-twist case. For both modes, increasing either parameter caused an amplifying effect in the force coefficients. This amplifying effect is characterized as an increase in the peak-to-peak amplitude of the force coefficients while maintaining the same overall trend. PIV measurements also showed that increasing the heaving amplitude resulted in increased circulation for all cases, lending support to this amplifying effect. Forward-twist stood as an outlier in this respect, with the lower frequency case k = 0.105 showing stronger effects than at the higher frequency case k = 0.209. PIV measurements also showed that increasing the heaving amplitude resulted in increased circulation for all cases.
Forward-twist displayed the most favorable traits at certain segments of the heaving cycle. During the upstroke, the spanwise twisting of forward-twist recovers some of the lift that is typically lost during this part of the cycle. Additionally, forward-twist exhibits a near-constant lift coefficient during the transition between upstroke and downstroke 0.4 < t * < 0.7. This suggests a more stable movement path during flapping locomotion, where the vehicle in question does not rise and fall as much throughout the heaving cycle. PIV analyses showed that morphing forward-twist limited circulation and LEV growth during the downstroke, resulting in a transient coefficient of drag C d ≈ 0 at the cost of the severely reduced lift. This limitation in circulation is attributed to the effect of spanwise twisting on the shear layer velocity of the hydrofoil. As the twisting in forward-twist shrinks the effective angle of attack, the shear layer velocity is reduced, limiting the amount of vorticity generated at the leading edge.
Backward-twist was demonstrated to affect the flow in a manner opposite to that of forward-twist morphing. It was shown to amplify the effects of rigid flapping, as its lift coefficient profile followed a sinusoidal trend similar to that of the rigid morphing mode; however, the lift coefficient was higher during the downstroke then dipped into the negative range C L < 0 during the upstroke. PIV analyses showed that backward-twist greatly increased the strength of circulation above the suction surface of the hydrofoil during the downstroke. The previously mentioned shear layer velocity mechanism applies here as well, but its implementation is opposite to that of forward-twist. The spanwise twisting in backwardtwist increased the effective angle of attack and, consequently the shear layer velocity. This caused stronger vorticity to form at the leading edge, resulting in a larger and stronger circulation region above the suction surface of the hydrofoil. Comparing the two morphing modes, it is hypothesized that executing a backward-twist during the downstroke and a forward-twist during the upstroke could provide a more favorable force coefficient distribution throughout the heaving cycle.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.