1. Introduction
For decades research has focused on studying the fluid dynamics of biological swimmers, both to better understand the underlying biology of aquatic animals, and to provide inspiration for developing innovative hydrodynamic propulsion technology (Smits Reference Smits2019). A salient feature of natural swimmers is the action of muscles, which are often distributed throughout an animal's propulsor (Flammang & Lauder Reference Flammang and Lauder2008; Adams & Fish Reference Adams and Fish2019). Through observations and measurements of animals, we know that swimmers can control their fin curvature, displacement and area as well as their stiffness (Fish & Lauder Reference Fish and Lauder2006). It stands to reason that swimmers may be able to utilize their muscles on the time scale of the oscillation of their propulsors to tune performance, e.g. by dynamically changing the stiffness of their propulsors; however, there is no definitive observation or consensus in the biological community that swimmers take advantage of their muscles in this way (Fish & Lauder Reference Fish and Lauder2021). In this work, we will show that, from a purely hydrodynamic perspective, time-varying stiffness leads to propulsive benefits over constant-stiffness propulsors.
The fluid dynamics of biological swimmers is rooted in the theory of rigid-wing flutter. Theodorsen (Reference Theodorsen1935) was the first to theoretically model the forces produced by oscillating foils in a fluid, which was later extended by Garrick (Reference Garrick1936) to analytically predict thrust and power for a rigid oscillating foil. Although these works focused on wing flutter, the connection to swimmers was clear. Later, Chopra & Kambe (Reference Chopra and Kambe1977) incorporated the effects of three-dimensionality via lifting line theory. Anderson et al. (Reference Anderson, Streitlien, Barrett and Triantafyllou1998) used particle image velocimetry to calculate thrust and power of harmonically oscillating foils and compared the results with inviscid theory predictions. More recently, Floryan et al. (Reference Floryan, Van Buren, Rowley and Smits2017) derived and experimentally validated propulsive scaling laws for rigid two-dimensional foils in pure heaving and pitching motions. This was extended to combined pitch-and-heave motions (Van Buren, Floryan & Smits Reference Van Buren, Floryan and Smits2019).
The passive flexibility (or elasticity) of an oscillating propulsor plays a key role in its propulsive performance. In a particularly influential work, Wu (Reference Wu1961) was one of the first to analytically consider passive flexibility, albeit through prescribed kinematics. Katz & Weihs (Reference Katz and Weihs1978, Reference Katz and Weihs1979) calculated the coupled fluid–structure interactions for a two-dimensional flexible foil. Since then, many analytical, experimental and computational studies have shown the influence of propulsor flexibility on characteristics such as the thrust and swimming efficiency, generally finding that flexibility can dramatically increase thrust and/or efficiency compared with a rigid propulsor (Alben Reference Alben2008; Michelin & Llewellyn Smith Reference Michelin and Llewellyn Smith2009; Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013; Moore Reference Moore2014, Reference Moore2015; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2014; Paraz, Schouveiler & Eloy Reference Paraz, Schouveiler and Eloy2016). In a particularly relevant work, Floryan & Rowley (Reference Floryan and Rowley2018) used small-amplitude theory to explore the role of resonance in constant-stiffness propulsors, finding that the benefits of resonance manifest in the thrust and power, but not necessarily in the propulsive efficiency. This analysis was extended in Goza, Floryan & Rowley (Reference Goza, Floryan and Rowley2020) to consider the role of nonlinearity in the fluid–structure system.
Few works have explored beyond simple uniform and passive flexibility. Floryan & Rowley (Reference Floryan and Rowley2020) considered the role of stiffness distribution, finding that concentrating stiffness toward the leading edge produced higher thrust, but lower efficiency compared with foils with stiffness concentrated away from the leading edge. Quinn & Lauder (Reference Quinn and Lauder2021) studied tuneable stiffness – that is, quasi-steady changes in stiffness – where it was shown that stiffness could be tuned to maximize desired performance parameters. To our knowledge, only two works consider time-varying stiffness at the same time scale as the kinematics. Hu et al. (Reference Hu, Wang, Wang and Dong2021) approximated a flexible plate by a rigid propulsor with a torsional spring at its leading edge heaving actively and pitching passively, respectively. By changing the stiffness of the torsional spring in time, both the thrust and propulsive efficiency could be enhanced in swimming. Shi, Xiao & Zhu (Reference Shi, Xiao and Zhu2020) used a finite-volume Navier–Stokes solver to study the effects of changing the flexibility of a nonlinear beam in the context of micro-air vehicles. This work provides valuable comparisons with the present study, as they also found (as we will show) that oscillating stiffness can strongly influence thrust and weakly influence the efficiency of the oscillating propulsor. However, they considered a smaller range of cases due to the time limitation of their solution method. Also, the stiffness oscillations considered were dramatic, changing stiffness by up to two orders of magnitude in a single cycle. Finally, the stiffness oscillations were non-sinusoidal, leading to to asymmetric motions and non-zero side forces. Non-zero side force is useful when considering flying as it leads to lift, however, here, our concern is rectilinear swimming, where asymmetric motions and non-zero side force would result in entering a manoeuvring condition.
In this work, we study the effects of time-periodic stiffness on the propulsive performance of a flexible plate. We solve a potential flow model that strongly couples the fluid and structural equations of motion. We use a pseudospectral method developed in Moore (Reference Moore2017) to solve the equations of motion, introducing time-varying flexibility through a Fourier series expansion in time. We then use Floquet theory to assess the stability of the system. We will show that the oscillating stiffness can strongly influence thrust and weakly influence the efficiency of the oscillating propulsor, indicating that animals or human-made swimmers would benefit from this capability. It should be noted that this paper uses external actuation at the leading edge of the plate. The model can be seen as an oscillating lifting surface attached to an external body or driving mechanism. This applies to thunniform swimmers, such as dolphins; flying animals and insects, such as bats; and also many uncrewed aerial/underwater vehicle applications. Thus, this model would not capture situations where the lifting surface is the entire body of the swimmer/flyer, such as an eel.
2. Problem statement and solution methods
Consider a two-dimensional, inextensible plate in a fluid flow, as sketched in figure 1. The plate has thickness 
$d$ and length 
$L$, and we assume that it is thin (
$d \ll L$) and that its maximum deflection is small, with its slope 
$|Y_x| \ll 1$.
A two-dimensional flat plate with time-varying Young's modulus moving through a fluid.
The deflection of the plate is then governed by the Euler–Bernoulli beam equation
\begin{equation} \rho_s d w Y_{tt} + E(t)I Y_{xxxx} = w {\rm \Delta} p, \end{equation}
where 
$Y$ is the transverse displacement of the plate, 
$\rho _s$ is its density, 
$E(t)$ is its time-varying Young's modulus, 
$I = w d^3/12$ is its second moment of area and 
$w$ is its width. The plate is immersed in a fluid which imparts a hydrodynamic load onto it, given by the pressure difference across the plate, 
${\rm \Delta} p$. Subscript 
$t$ denotes differentiation with respect to time, and subscript 
$x$ denotes differentiation with respect to the streamwise coordinate.
The fluid is inviscid and incompressible, with density 
$\rho _f$. Far from the plate, the fluid moves with a free-stream velocity 
$\boldsymbol U = U \boldsymbol i$, where 
$\boldsymbol i$ is the unit vector in the 
$x$ direction. Conservation of mass and momentum for the fluid leads to
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\rho_f ( \boldsymbol u_t + U \boldsymbol u _x) ={-} \boldsymbol{\nabla} p, \end{gather}$$
where 
$\boldsymbol u = u \boldsymbol i + v \boldsymbol j$ is the perturbation velocity induced by the motion of the plate, and 
$\boldsymbol j$ is the unit vector in the transverse direction. In obtaining (2.2b), we have assumed that 
$|\boldsymbol u| \ll U$, which is consistent with our assumption of small-amplitude deflections of the plate. The limitations of these small-amplitude assumptions are small angle deflections and amplitudes, and attached flow.
We non-dimensionalize the equations of motion using the half-length of the plate 
$L/2$ as the length scale, the free-stream velocity 
$U$ as the velocity scale and the convective time 
$L/(2U)$ as the time scale. The non-dimensional equation for the plate is
\begin{equation} 2RY_{tt}+ \tfrac{2}{3} S(t) Y_{xxxx} = {\rm \Delta} p, \end{equation}and the non-dimensional equations for the fluid are
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol u_t + \boldsymbol u _x = \boldsymbol{\nabla} \phi, \end{gather}$$where
\begin{equation} R = \frac{\rho_s d}{\rho_f L},\quad S(t) = \frac{E(t) d^3}{\rho_f U^2 L^3},\quad \phi = p_\infty - p. \end{equation}
The function 
$\phi$ is Prandtl's acceleration potential (Wu Reference Wu1961). Now 
$x$, 
$t$, 
$Y$, 
$\boldsymbol u$ and 
$p$ are non-dimensional, with 
$x = -1$ corresponding to the leading edge of the plate, and 
$x = 1$ corresponding to the trailing edge. The mass ratio 
$R$ is the ratio of a characteristic mass of the plate to a characteristic mass of fluid, and the stiffness ratio 
$S$ is the ratio of a characteristic bending force to a characteristic fluid force. The stiffness ratio, its inverse and variations of it are sometimes called the Cauchy number (Cermak & Isyumov Reference Cermak and Isyumov1999; de Langre Reference de Langre2008) or the elastohydrodynamic number (Schouveiler & Boudaoud Reference Schouveiler and Boudaoud2006).
To close the system of equations, we specify the boundary conditions. The fluid satisfies the no-penetration boundary condition and the Kutta condition
$$\begin{gather} v|_{x \in [{-}1,1],y=0} = Y_{t} + Y_x, \end{gather}$$
$$\begin{gather}|v||_{(x,y) = (1,0)} < \infty. \end{gather}$$
We specify heaving and pitching motions 
$h$ and 
$\theta$, respectively, at the leading edge of the plate, while the trailing edge is free (zero force and torque), resulting in the boundary conditions,
\begin{equation} Y({-}1,t) = h(t),\quad Y_{x}({-}1,t) = \theta(t),\quad Y_{xx}(1,t) = 0,\quad Y_{xxx}(1,t) = 0. \end{equation} Since we are interested in locomotion, we consider periodic actuation of the leading edge; that is, 
$h(t)$ and 
$\theta (t)$ are (zero-mean) periodic functions of 
$t$ with a non-dimensional angular frequency 
$\sigma = {\rm \pi}f L / U$, where 
$f$ is the dimensional ordinary frequency. Similarly, we consider a time-periodic stiffness. With our focus on forward propulsion, the upstroke and downstroke must be mirror images of each other in order for the mean side force to be zero. It can be shown that this requires that the stiffness vary at twice the frequency of the deflection of the plate. Intuitively, the stiffness must be the same during the upstroke and downstroke, which can only be true if it varies at twice the frequency of the plate's deflection.
To solve the system of equations for the kinematics of the plate, we use the method described by Moore (Reference Moore2017), adapting it to account for time-periodic stiffness; we describe the method in Appendix A. The method assumes that the kinematics are time periodic, with any transients decaying to zero. We will test this assumption by performing a Floquet analysis. More importantly, the Floquet analysis will provide physical insight into the problem. The Floquet analysis is adapted from the eigenvalue problem described by Floryan & Rowley (Reference Floryan and Rowley2018, Reference Floryan and Rowley2020) and uses a key idea from Kumar & Tuckerman (Reference Kumar and Tuckerman1994); the details are in Appendix B.
The motion of the plate induces a pressure difference across it. The projection of the pressure difference onto the horizontal direction contributes – together with a suction force at the leading edge – to a propulsive thrust force. Therefore, the energy put into the system by actuating the leading edge is converted into propulsive thrust, given by
\begin{equation} C_T = \int_{{-}1}^{1}{\rm \Delta} p Y_x\,{{\rm d}\kern0.06em x} + C_{TS}, \end{equation}
where 
$C_{TS}$ is the leading-edge suction, a formula for which is given by Moore (Reference Moore2017). The power input is
\begin{equation} C_P ={-}\int_{{-}1}^{1} {\rm \Delta} p Y_t\,{{\rm d}\kern0.06em x}. \end{equation}Finally, the Froude efficiency is defined as the ratio of the time-averaged thrust output to the time-averaged power input (i.e. how much of the power input is converted to propulsive thrust)
\begin{equation} \eta = \frac{\overline{C_T}}{\overline{C_P}}, \end{equation}where the overbar denotes averaging in time.
In this work, we restrict ourselves to leading-edge actuation and stiffnesses that are sinusoidal in time
$$\begin{gather} h(t) = \tfrac{1}{2}( h_0\,\exp(j \sigma t) + h_0^*\exp({-j \sigma t})), \end{gather}$$
$$\begin{gather}\theta(t) = \tfrac{1}{2}(\theta_0 \,\exp(j \sigma t) + \theta_0^* \exp({-j \sigma t})), \end{gather}$$
$$\begin{gather}S(t) = \bar S + \frac{\bar S}{2}(S_0 \exp({2j \sigma t}) + S_0^*\exp({-2j \sigma t})) , \end{gather}$$
where 
$h_0, \theta _0, S_0 \in \mathbb {C}$, 
$j = \sqrt {-1}$ and a superscript 
$*$ denotes complex conjugation. The formulation in the appendices, however, is valid for generic smooth periodic functions of time. We will make frequent use of the parameter 
$\phi _S = \text {arg}(S_0)$, the phase of the stiffness oscillation. Note that 
$|S_0|$ gives the amplitude of the stiffness oscillation as a fraction of the mean stiffness; for example, 
$S_0 = 0.5$ means that the stiffness oscillates with an amplitude that is 50 % of the mean stiffness. For a physically meaningful (i.e. positive) stiffness, we require 
$|S_0| < 1$.
Throughout, we fix the mass ratio to a low value of 
$R = 0.01$, appropriate for thin, neutrally buoyant biological swimmers. To build intuition for the effects of time-varying stiffness, we will extensively study the case with 
$\bar S = 20$; unless otherwise noted, this is the value we use for the mean stiffness. We will consider cases where the plate is in pure heave, in pure pitch and in combined heave and pitch. Throughout, we set 
$h_0 = 1$ and 
$\theta _0 = 0.5$, although we will add a phase offset between heave and pitch for combined motions.
A discrepancy the reader may have noted with our small-amplitude and deflection assumptions is the high values for 
$h_0$ and 
$\theta _0$ to 1 and 
$0.5$, respectively. These will produce very high deflections not consistent with the linearization of the system. The reason for these values is the ability to easily re-scale the outputs based on the input due to the linearity of the problem. If we want to find the kinematics for a heave input of 
$0.05$, then we take the deflection of this paper and multiply it by 
$0.05$ to find the new deflection at that input, and multiply thrust or power by 
$0.05^2$ to find the new thrust or power. For pitch, we would multiply our results by 
$0.05/0.5$ to find the new deflections, and by 
$(0.05/0.5)^2$ to find the new thrust or power. Efficiency stays the same for all scaled input cases. The results are only valid if the rescaling satisfies the small-amplitude assumption.
Kinematics of a heaving plate actuated at its first resonant frequency of 
$\sigma = 3.1$. Shown are a plate with constant stiffness (solid black), time-periodic stiffness in phase with the motion (
$\phi _S = 0$; dashed blue) and time-periodic stiffness out of phase with the motion (
$\phi _S = {\rm \pi}$; dashed red). The coloured arrows represent the difference between the instantaneous stiffness and the mean stiffness. The absence of arrows indicates the three plates have the same stiffness in the below snapshot. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$R = 0.01$ and 
$|S_0| = 0.5$.
3. Results and discussion
Heaving and pitching motions have fundamentally different thrust mechanisms (Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Van Buren, Floryan & Smits Reference Van Buren, Floryan and Smits2020), with thrust in heave coming from lift-based circulatory forces, and in pitch coming from added-mass acceleration forces. To start, we consider a periodically heaving plate moving through a fluid. After completing our heave-only analysis, we will analyse pitch-only and heave-and-pitch motions.
3.1. Heave-only motions
First, we familiarize ourselves with how time-varying flexibility impacts the plate's kinematics. Figure 2 shows the kinematics of three plates during one cycle of motion. The plates are all actuated at the same frequency – the first resonant frequency of the plate with constant stiffness, 
$\sigma = 3.1$. The reference case with constant stiffness is compared against a plate with time-varying flexibility that is in phase with the motion (
$\phi _S = 0$), meaning it is stiff at the turnaround and flexible at mid-stroke, and a plate with time-varying flexibility that is out of phase with the motion (
$\phi _S = {\rm \pi}$), meaning it is flexible at the turnaround and stiff at the mid-stroke. For a plate with constant stiffness, the peak deflection occurs around the mid-stroke, where the lateral velocity is highest. For the 
$\phi _S=0$ case, the pressure on the plate is reduced throughout the mid-stroke with more deflection, and the increase in stiffness towards the turnaround causes it to catch back up to the reference case. Conversely, the 
$\phi _S={\rm \pi}$ case has its stiffness reduced at the turnaround, ultimately achieving the largest trailing-edge deflection of all heave cases, and then recovering throughout the mid-stroke with increased stiffness. The stiffness oscillation essentially adds a phase lag relative to the motion of the constant-stiffness plate. When 
$\phi _S = {\rm \pi}$, the motion leads that of the constant-stiffness plate, whereas when 
$\phi _S = 0$, the motion lags that of the constant-stiffness plate.
In figure 3, we show how time-varying stiffness affects average thrust and efficiency. The frequency range is centred about the first resonant frequency of the plate with constant stiffness. Generally, adding periodic stiffness leads to a continuous and substantial increase in thrust as the stiffness oscillation amplitude increases, with up to a 35 % increase in thrust when 
$|S_0| = 0.5$. For 
$\phi _S = 0$, the resonance is shifted to higher frequencies, while for 
$\phi _S = {\rm \pi}$ the resonance is shifted lower compared with the constant stiffness case. Past the resonant frequency, the 
$\phi _S = {\rm \pi}$ case yields slightly lower thrusts than the baseline case.
Thrust (a) and efficiency (b) for a heaving plate with gradually increasing stiffness oscillation amplitude. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$ and 
$R = 0.01$.
For efficiency, we observe the opposite behaviour (figure 3b). Performance benefits for 
${\phi _S = 0}$ occur below the baseline resonant frequency, while for 
$\phi _S = {\rm \pi}$ they occur above the baseline resonant frequency. For both phases of the stiffness oscillation, there are frequencies yielding greater efficiency than the constant-stiffness plate, as well as frequencies yielding lower efficiency. The effect of time-varying stiffness on the efficiency is much more mild than it is for thrust, however, as time-varying stiffness only changes the efficiency by 
${\rm \Delta} \eta \leq 0.05$. This suggests that time-varying stiffness may be a strategy to substantially increase thrust without much, if any, penalty in efficiency, as was similarly shown in Shi et al. (Reference Shi, Xiao and Zhu2020).
To better understand the effects of time-varying stiffness, we turn to the time histories of the performance characteristics over the course of an actuation cycle. Figure 4 shows the instantaneous thrust, side force and power coefficients over one period of the motion. For reference, we also plot the leading-edge kinematics and plate stiffness. The plate with constant stiffness exhibits purely sinusoidal thrust, power and side force since frequencies are uncoupled in the small-amplitude limit. Time-periodic stiffness, however, causes cross-frequency coupling, leading to non-sinusoidal behaviour. Note that, because the stiffness oscillates at twice the frequency of motion, the side force remains symmetric, ensuring there is no mean side force (in a real system, this would lead to manoeuvring).
Instantaneous thrust, power and side force coefficients for a heaving plate. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$R = 0.01$, 
$\sigma = 3.1$ and 
$|S_0| = 0.5$. Shown are a plate with constant stiffness (solid black), time-periodic stiffness in phase with the motion (
$\phi _S = 0$; dashed blue) and time-periodic stiffness out of phase with the motion (
$\phi _S = {\rm \pi}$; dashed red). (a) Kinematic input at the leading edge (solid black, left ordinate), and the stiffness distribution (blue and red dashed, right ordinate) for one stroke. (b) Thrust coefficient as a function of time. (c) Power coefficient as a function of time. (d) Side force coefficient as a function of time.
The plate whose stiffness oscillates out of phase with the kinematics (
$\phi _S = {\rm \pi}$) produces the most thrust at the mid-stroke. During the turnaround the plate becomes the most flexible – this helps to relieve the power consumption which trends highest before the plate reverses direction. The opposite happens for the plate whose stiffness is in phase with the kinematics (
$\phi _S = 0$). It produces the most thrust closer to the turnaround, and the power relief comes when the plate is most flexible at the maximum plunge velocity. For both cases, higher thrust trends towards the times of higher stiffness, and lower power consumption trends towards the times of lower stiffness. The phase of the stiffness oscillation dictates when in the cycle the plate is most stiff (promoting thrust) or most flexible (relieving power). Thus, by timing these stiffness changes at opportune times during the cycle (e.g. becoming flexible at a time when power is highest), the performance can be significantly enhanced.
It is clear that the timing of the stiffness oscillation is important. However, until this point, only two stiffness oscillation phases have been considered: in phase and out of phase with the kinematics. We investigate whether there is a particular phase that maximizes thrust or efficiency. In figure 5, thrust and efficiency are plotted on a polar plot with frequency 
$\sigma$ on the radial axis and stiffness phase offset 
$\phi _s$ on the azimuthal axis. Both thrust and efficiency are shown relative to the constant-stiffness case. Time-varying stiffness increases thrust the most when 
$\phi _S$ is between 
${\rm \pi} /2$ and 
$2{\rm \pi} /3$, and it increases efficiency the most near 
$\phi _s = 3{\rm \pi} /2$. Again, the changes in efficiency are modest. It is important here to also highlight the importance of the kinematic frequency – for a given phase of the stiffness oscillation, there can be performance boosts and hindrances depending on the frequency of oscillation.
Impact of stiffness phase offset (
$\phi _S$; azimuthal axis) and kinematic frequency (
$\sigma$; radial axis) on thrust and efficiency for a heaving plate, relative to the constant-stiffness case. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$R = 0.01$ and 
$|S_0| = 0.5$.
We seek to understand why time-varying stiffness confers greater thrust than constant stiffness. For rigid and passively flexible plates, prior work has shown that the relevant velocity scale for thrust generation is the characteristic lateral velocity of the plate – not the free-stream velocity traditionally used in aerodynamics – and that thrust scales quadratically with the lateral velocity scale (Gazzola, Argentina & Mahadevan Reference Gazzola, Argentina and Mahadevan2014; Van Buren et al. Reference Van Buren, Floryan, Wei and Smits2018). The lateral velocity scale is usually taken to be the amplitude of the trailing edge's velocity, proportional to 
$\sigma A$, where 
$A$ is the amplitude of the trailing edge's displacement. In figure 6(a) we account for varying lateral velocity scales for the cases with time-varying stiffness by dividing the thrust by 
$A^2$ (the different cases share the same value of 
$\sigma$, so normalizing by 
$A^2$ captures the effects of the lateral velocity scale). By accounting for the lateral velocity scale, the curves become reasonably collapsed (cf. figure 3), indicating that time-varying stiffness increases thrust principally by increasing the characteristic lateral velocity. However, near the regions of plate resonance (
$3 < \sigma < 4$), the rescaled thrust curves still deviate from each other. Time-varying stiffness with phase 
$\phi _S = {\rm \pi}$ produces more thrust than one would expect from our scaling argument, whereas the opposite holds for 
$\phi _S = 0$. The lateral velocity scale, therefore, does not tell the whole story; other factors are afoot. (We tried other trailing-edge velocity scales as well, e.g. the maximal trailing-edge velocity, but none perfectly collapsed the thrust curves.)
(a) Total thrust, (b) thrust due to projected pressure difference and (c) thrust due to leading-edge suction, all scaled by the square of the maximum trailing-edge amplitude 
$A$.
To further explore the scaling breakdown, we decompose the thrust into its components (see (2.8)): (i) the projected pressure difference; and (ii) the leading-edge suction. We plot the two components, rescaled to account for the lateral velocity scale, in figure 6(b,c), respectively. The thrust due to the projected pressure difference stays completely collapsed throughout the resonance region for all cases, but the leading-edge suction deviates in the same region as in figure 6(a). Thus, for cases with time-varying stiffness, the deviation of thrust from the lateral velocity scaling is attributed to the leading-edge suction. The detailed solution for the leading-edge suction can be found in Wu (Reference Wu1961). Based on (21), (23b), (35) and (61) in Wu (Reference Wu1961), the leading-edge suction depends on the spatial distribution of the lateral velocity along the plate and on the circulation around the plate (which itself can be related to the spatial distribution of the lateral velocity along the plate). The lateral velocity scaling ignores the detailed spatial distribution of the lateral velocity. We conclude that time-varying stiffness changes thrust not only by changing the characteristic lateral velocity scale, but also by changing the detailed spatial distribution of the lateral velocity, which acts through the leading-edge suction to alter thrust.
Before moving to more in-depth analysis, we consider the implications of these results on real-world systems. We have shown to this point that time-varying stiffness leads to large increases in thrust with moderate changes in efficiency. Furthermore, the changes in thrust do not follow our traditional understanding of velocity scaling, i.e. one cannot assume that the dynamic stiffness is merely changing the trailing-edge amplitude and a passively flexible plate with equal trailing-edge amplitude would perform equally well. As a result, from a technological perspective, while it may be a greater challenge to design a system with actively changing flexibility, there may be intrinsic benefits to these types of propulsors. Furthermore, in systems where oscillation amplitude or frequency are limited, oscillating stiffness may be a realizable tool in biological or human-made underwater swimmers to improve swimming speed. When considering swimming efficiency, one must take into account that changing the stiffness of the system will require additional energy.
3.1.1. Higher resonance modes and instability
We now explore a much wider range of mean plate stiffnesses and frequencies encompassing higher-order resonant frequencies. The change in swimming performance due to time-periodic stiffness, with respect to the constant-stiffness case, is shown in figure 7 for 
$\phi _S = 0$ and 
${\rm \pi}$. The oscillation in stiffness modifies the resonant frequencies from those of the constant-stiffness plate, as we saw in figure 3, leading to sharp bands of increased and decreased thrust; these appear as adjacent narrow red and blue strips in figure 7. In a linear time-invariant system, resonant frequencies are related to the imaginary parts of the eigenvalues of the system, whereas in linear time-periodic systems, they are related to the imaginary parts of the Floquet exponents (Wereley Reference Wereley1990). In general, the two are different, leading to the modified resonant frequencies that we observe. Away from resonance, thrust is modified little by the oscillation in stiffness. The changes in efficiency are much broader over the stiffness–frequency plane and have features that align with resonant frequencies; they are, however, very mild. Generally, thrust increases where efficiency decreases, but increasing 
$|S_0|$ can greatly increase the thrust production with minimal impact on efficiency.
Impact of time-varying stiffness on the thrust (a,c) and efficiency (b,d), relative to the constant-stiffness case, for stiffness oscillations that are in phase with the motion (a,b), and 180
$^\circ$ out of phase with the motion (c,d). The parameters used are 
$h_0 = 1$, 
$R = 0.01$ and 
$|S_0| = 0.5$. (a) 
$\log {C_T/C_{T_{S_0=0}}}$, 
$\phi _{S} = 0$; (b) 
${\rm \Delta} \eta$, 
$\phi _{S} = 0$; (c) 
$\log {C_T/C_{T_{S_0=0}}}$, 
$\phi _{S} = {\rm \pi}$; (d) 
${\rm \Delta} \eta$, 
$\phi _{S} = {\rm \pi}$.
Perhaps the most interesting features in figure 7 are the hollowed peaks in the thrust, resembling the eye of a needle. They are more prominent at higher frequencies. To clarify this behaviour, in figure 8 we show the thrust and efficiency along a slice in the stiffness–frequency plane, taking 
$\bar {S} = 20$ and centring the frequency range about the second resonant frequency (this is the same parameter case as shown in figure 3). As the stiffness oscillation amplitude is gradually increased, there is a critical value of 
$|S_0|$ at which a single resonant peak in thrust bifurcates into two sharp peaks centred about the baseline resonant frequency, with the distance between the peaks increasing as 
$|S_0|$ increases. The bifurcated peaks are very sharp, indicative of natural frequencies with very little damping.
Values of 
$\log C_T$ (a) and efficiency (b) for a heaving plate with gradually increasing stiffness oscillation amplitude. As stiffness oscillation amplitude 
$|S_0|$ increases, a single resonant peak bifurcates into two. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$ and 
$R = 0.01$.
To investigate this possibility, we perform a Floquet analysis. Representative results from the Floquet analysis are shown in figure 9. There, we have plotted the neutral curves in the 
$\bar S$–
$|S_0|$ plane on top of contours of 
$\log C_T$. When the stiffness oscillation amplitude 
$|S_0|$ is below the neutral curves, the system is stable. Conversely, the system is unstable when 
$|S_0|$ is above the neutral curves, leading to unbounded growth in the plate's deflection. We see the emergence of tongues of instability, which are characteristic features of parametrically excited systems such as Mathieu's equation (Nayfeh & Mook Reference Nayfeh and Mook1995) and Faraday waves (Faraday Reference Faraday1831; Miles & Henderson Reference Miles and Henderson1990). The physical mechanism of the instability is the same as for parametrically excited oscillators: the time-varying stiffness creates an internal forcing. When this forcing is adding energy into the system at a higher rate than is being dissipated, this leads to a net increase of energy in the system, i.e. an instability. Because the system is damped (due to energy lost to the wake via a thin sheet of vorticity), the instability emerges at a non-zero value of 
$|S_0|$. For a fixed value of 
$|S_0|$, as we vary the mean stiffness 
$\bar S$, we cross into and then out of the unstable region. At the boundaries, the thrust exhibits sharp peaks, which are explained by the theory of forced linear time-periodic systems (Wereley Reference Wereley1990). This explains the double resonant peaks that we observed in figures 7 and 8. (Although we vary 
$\sigma$ at a fixed value of 
$\bar S$ in figure 8, we can see in figure 7 that we encounter double resonant peaks whether we vary 
$\bar S$ at a fixed value of 
$\sigma$, or we vary 
$\sigma$ at a fixed value of 
$\bar S$; physically, the origin of the double resonant peaks is the same.) No physical significance should be assigned to the thrust in the unstable region between the double resonant peaks since the thrust was calculated under the assumption of a stable system. We have verified for conditions other than those used to generate figure 9 that a splitting of one resonant peak into two coincides with the emergence of an instability; we conjecture that every region between double resonant peaks in figures 7 and 8 is actually unstable. While instability usually has a negative connotation in engineering applications, these unstable regions could potentially lead to greatly enhanced propulsive performance since an instability would produce large deflections from small actuation. Our linear method cannot capture the saturation of the instabilities, but we believe that exploring the unstable regions via experiments or nonlinear simulations is a promising direction.
Neutral stability curves overlaid on contours of 
$\log C_T$. The parameters used are 
$R = 0.01$ to calculate the neutral stability curves, and 
$h_0 = 1, \sigma = 20$ and 
$\phi _S = 0$ to calculate the thrust values.
3.2. Pitch-only and pitch-and-heave motions
We now consider how pitching the leading edge changes the results. As stated previously, pitching and heaving are fundamentally different in their thrust generation mechanism (Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Van Buren et al. Reference Van Buren, Floryan and Smits2020), with the former utilizing added-mass forces and the latter using lift-based or circulatory forces. Additionally, combined pitching and heaving motions are the most biologically relevant and also the best performing (Wu Reference Wu2011; Van Buren et al. Reference Van Buren, Floryan and Smits2019). In this section, we study two cases: (i) purely pitching, and (ii) combined pitching and heaving with pitch lagging heave by 
${\rm \pi} / 2$, the most ‘fish-like’ motion that is also generally the most efficient (Triantafyllou, Triantafyllou & Yue Reference Triantafyllou, Triantafyllou and Yue2000; Van Buren et al. Reference Van Buren, Floryan and Smits2019).
Figure 10 shows the kinematics for both cases during one cycle of motion. The kinematics are qualitatively the same as for the purely heaving case in figure 2. Adding an oscillatory component to the stiffness causes a lead/lag effect on the kinematics which dictates at what point in the stroke thrust is generated. When the stiffness is in phase with the motion, rapid trailing-edge motion occurs near the turnaround where the plate transitions from flexible to stiff, whereas when the stiffness is out of phase with the motion, the rapid trailing-edge motion occurs through the midpoint of the cycle. Even for the pitching and heaving case – which specifically reduces the side force on the plate (Van Buren et al. Reference Van Buren, Floryan and Smits2019) – we see the stiffness oscillations have similar impact in both magnitude and timing.
Analogous to figure 2, but for a pitching plate (a) and a pitching and heaving plate (b). Shown are a plate with constant stiffness (solid black), time-periodic stiffness in phase with the motion (
$\phi _S = 0$; dashed blue) and time-periodic stiffness out of phase with the motion (
$\phi _S = {\rm \pi}$; dashed red). The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$\theta _0 = -0.5j$, 
$R = 0.01$ and 
$|S_0| = 0.5$.
Thrust and efficiency of a purely pitching plate (a,b) and a pitching and heaving plate (c,d) as a function of driving frequency 
$\sigma$. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$\theta _0 = -0.5j$ and 
$R = 0.01$: (a) 
$C_T$, 
$h = 0, \theta = 0.5$; (b) 
${\rm \Delta} \eta$, 
$h = 0, \theta = 0.5$; (c) 
$C_T$, 
$h = 1, \theta = 0.5$; (d) 
${\rm \Delta} \eta$, 
$h = 1, \theta = 0.5$.
Impact of stiffness phase offset (
$\phi _S$; azimuthal axis) and kinematic frequency (
$\sigma$; radial axis) on thrust and efficiency for a pitching (a,b) and a pitching and heaving (c,d) plate, relative to the constant-stiffness case. Note the white unit circles on the efficiency contours. They are whited out due to large negative efficiency values, which detract from the regions of interest. The parameters used are 
$\bar S = 20$, 
$h_0 = 1$, 
$\theta _0 = -0.5j$, 
${R = 0.01}$ and 
$|S_0| = 0.5$. (a) 
$\log {C_T/C_{T_{S_0=0}}}$, 
$h = 0, \theta = 0.5$; (b) 
${\rm \Delta} \eta$, 
$h = 0, \theta = 0.5$; (c) 
$\log {C_T/C_{T_{S_0=0}}}$, 
$h = 1, \theta = 0.5$; (d) 
${\rm \Delta} \eta$, 
$h = 1, \theta = 0.5$.
Figure 11 shows the thrust and efficiency for the two cases in a frequency range centred about the first resonant frequency of the plate with constant stiffness. For both the pitch and pitch-and-heave cases, we see that the oscillating stiffness has a strong impact on the peak thrust near resonance. The efficiency is less impacted by the stiffness oscillation and switches from being greater than the efficiency of a constant-stiffness plate to less than it across the resonant frequency – depending on the phase of the stiffness oscillation – which is similar to the behaviour of the purely heaving plate in figure 3. For the purely pitching plate, thrust is much more impacted by the oscillating stiffness when it is out of phase with the kinematics, 
$\phi _S = 0$. For pitching, the blue plate lags behind the red and black plates at the turnaround, and becomes the least stiff at this moment. The blue plate accelerates the fastest between frames 2 and 4 in figure 10, which corresponds to the fastest angular velocity at the leading edge. The blue plate had the highest trailing-edge amplitude, because the phase difference between the leading-edge input and the trailing-edge deflection and its lower stiffness allows the blue plate to reach a higher trailing-edge amplitude before feeling the effects of the acceleration at the leading edge. The blue plate then becomes stiffer as the effect of the acceleration at the leading edge reaches the trailing edge. This combination of high acceleration, high trailing-edge amplitude and high stiffness creates high thrust. The red plate's stiffness distribution is misaligned with the phase difference caused by pure pitching, and therefore reaches a lower maximum trailing-edge amplitude, and benefits less from the high acceleration of the leading edge between frames 2 and 4 in figure 10. For a better visualization of the plate dynamics and how it ties into the performance enhancements, refer to the online supplemental movies available at https://doi.org/10.1017/jfm.2023.166 where we show the heave, pitch and heave plus pitch kinematics.
Finally, we consider the role of the phase of stiffness oscillation in more detail. Figure 12 shows the change in thrust and efficiency relative to the constant-stiffness case. For the heaving plate, the ideal phase for thrust was approximately between 
${\rm \pi} /2$ and 
$2{\rm \pi} /3$; for the pitching plate, however, the ideal phase is approximately 
${\rm \pi} /3$. The ideal phase for the combined pitching and heaving plate is approximately 
${\rm \pi} /2$, directly between the isolated pitching and heaving cases. Generally, the efficiency has opposing behaviour to the thrust, with increased thrust coinciding with decreased propulsive efficiency.
4. Conclusions
In this work, we explore the impact of time-periodic stiffness on an oscillating plate in a free stream – a model for swimmers. The stiffness oscillations were varied in amplitude and phase and consistently compared with the baseline constant-stiffness case. The stiffness oscillation frequency is fixed at twice the kinematic frequency to enforce zero-mean side force, which is required for rectilinear swimming. The fluid–structure interaction model is solved using a spectral method first introduced in Moore (Reference Moore2017), modified to account for time-periodic stiffness.
It is shown that the oscillation in stiffness has a significant impact on the thrust and kinematics. The impact on thrust is most prominent, with time-periodic stiffness increasing thrust by up to 35 % at the first resonant peak. The impact on efficiency is, on the other hand, mild, with the efficiency changing by 
${\rm \Delta} \eta \le 0.05$. The changes in thrust and efficiency are often negatively correlated. These performance changes are consistent across heaving, pitching and combined pitch-and-heave motions.
The performance alteration due to time-varying stiffness is strongly linked to the phase of the stiffness oscillation. This is because thrust and side forces are produced at different stages throughout the kinematic cycle, and whether the plate is more or less stiff at those stages either yields enhanced thrust or power reduction. This is especially evident when comparing pure heaving motions with pure pitching motions, which have different phase differences between the leading-edge input and the trailing-edge deflection. Further analysis of the thrust behaviour with respect to the trailing-edge velocity scaling used in rigid and passively flexible systems indicates that the oscillating stiffness leads to additional physics that are not completely captured by the trailing-edge velocity scaling. This is due to the oscillating stiffness changing the spatial distribution of the plate's lateral velocity, which directly leads to a difference in the circulation around the plate and consequently the leading-edge suction.
When the kinematic frequency and amplitude of stiffness oscillation are large enough, instabilities emerge in regions of parameter space centred about resonant frequencies of the constant-stiffness plate. Our linear method cannot capture the saturation of the instabilities, but we anticipate that the unstable regions of parameter space may yield enhanced propulsive performance. Exploring the unstable regions via experiments or nonlinear simulations is a promising future research direction.
While there may not yet be concrete evidence of biological swimmers actively changing their stiffness on the time scale of their kinematic frequency, we have presented strong evidence that there would be a hydrodynamic benefit to doing so. This may be a promising avenue to pursue when designing robotic swimmers.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.166.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Method of solution
We assume the kinematics 
$Y$, the hydrodynamic load 
$q = {\rm \Delta} p$ and the stiffness ratio 
$S$ are periodic in time and can be expressed as Fourier series
$$\begin{gather} Y(x,t)= \sum_{m={-}\infty}^{\infty} \hat{y}_{m}(x) \exp({\, j m \sigma t}), \end{gather}$$
$$\begin{gather}q(x,t) = \sum_{m={-}\infty}^{\infty}\hat{q}_m(x)\exp({\, j m \sigma t}), \end{gather}$$
$$\begin{gather}S(t) = \sum_{m={-}\infty}^{\infty}\hat{S}_m \exp({\, j m \sigma t}). \end{gather}$$Substituting these expressions into (2.3) gives
\begin{align} & \sum_{m={-}\infty}^{\infty} \hat{q}_m(x) \exp({\, j m \sigma t}) + 2R \sum_{m={-}\infty}^{\infty}m^2 \sigma ^2 \hat{y}_m (x)\exp({\, j m \sigma t}) \nonumber\\ &\quad -\tfrac{2}{3}\sum_{m={-}\infty}^{\infty}\sum_{k={-}\infty}^{\infty} \hat{S}_{m - k} \hat{y}_k''''(x) \exp({\, j m \sigma t}) = 0, \end{align}
where a prime denotes differentiation with respect to 
$x$. Separating the Fourier components yields the following system of ordinary differential equations:
\begin{equation} \hat{q}_m (x) + 2m^2 \sigma^2 R \hat{y}_m (x) - \tfrac{2}{3}\sum_{k={-}\infty}^{\infty} \hat{S}_{m-k} \hat{y}''''_{k} (x) = 0,\quad \forall\ m \in \mathbb{Z}. \end{equation}Expanding the heaving and pitching motions into Fourier series as
$$\begin{gather} h(t) = \sum_{m={-}\infty}^{\infty} \hat{h}_m \exp({\, j m \sigma t}), \end{gather}$$
$$\begin{gather}\theta(t) = \sum_{m={-}\infty}^{\infty} \hat{\theta}_m \exp({\, j m \sigma t}), \end{gather}$$
and substituting them and (A1a) into (2.7a–d) gives boundary conditions for 
$\hat y_m$,
\begin{equation} \hat{y}_m({-}1) = \hat h_m,\quad \hat{y}'_m({-}1) = \hat{\theta}_m,\quad \hat{y}''_m(1) = \hat{y}'''_m(1) = 0, \quad \forall\ m \in \mathbb{Z}. \end{equation}We also require that all the functions are real, leading to the reality conditions
\begin{equation} \left.\begin{gathered} \hat y_m(x) = \hat y_{{-}m}^*(x),\quad \hat q_m(x) = \hat q_{{-}m}^*(x),\quad \hat S_m = \hat S_{{-}m}^*,\\ \hat h_m = \hat h_{{-}m}^*,\quad \hat \theta_{m} = \hat \theta_{{-}m}^{*},\quad \forall\ m \in \mathbb{Z}, \end{gathered}\right\} \end{equation}
where the superscript 
$*$ again denotes complex conjugation.
Given the hydrodynamic load 
$q$, the stiffness ratio 
$S$ and the heaving and pitching actuation at the leading edge, we can solve for the kinematics 
$Y$ by solving the system of coupled ordinary differential equations given by (A3) along with their boundary conditions given by (A5a–c). What remains is to couple the solid and fluid mechanics by expressing the hydrodynamic load in terms of the kinematics, which we do next.
The pressure difference across the plate creates the hydrodynamic load, and is related to Prandtl's acceleration potential by
\begin{equation} {\rm \Delta} p = \phi_{bottom} - \phi_{top}. \end{equation}
Taking the divergence of (2.4b) and using incompressibility shows that 
$\phi$ is harmonic, implying the existence of a harmonic conjugate 
$\psi$. We define the complex acceleration potential
\begin{equation} F(x,y,t) = \phi(x,y,t) + i \psi(x,y,t), \end{equation}
where 
$\textrm {i} = \sqrt {-1}$. It will be convenient to work with the complex variable 
$z = x + \textrm {i}y$. Note that 
$F$ is analytic in 
$z$. We conformally map the physical 
$z$ plane to the exterior of the unit disk in the 
$\zeta$ plane using
\begin{equation} z = \frac{1}{2}\left( \zeta + \frac{1}{\zeta}\right). \end{equation}
Because the acceleration potential 
$F$ is conformally invariant under (A9), 
$F$ is analytic in 
$\zeta$ and can be represented by a multipole expansion
\begin{equation} F = {\rm i}\left(\frac{a_0(t)}{\zeta + 1} + \sum_{k = 1}^{\infty}\frac{a_k(t)}{\zeta ^k} \right), \end{equation}
where the coefficients 
$a_k$ are real (Wu Reference Wu1961). The first term represents the singularity at the leading edge, and the infinite series represents an analytic function that is regular on and outside the unit circle, decaying in the far field (Wu Reference Wu1961).
Continuing with our assumption of a time-periodic flow, we expand the coefficients in the multipole expansion into Fourier series
\begin{equation} a_k(t) = \sum_{m={-}\infty}^{\infty} \hat{a}_{k,m}\exp({\, j m \sigma t}), \quad \forall\ k \in \mathbb{W}. \end{equation}Substituting (A11) into (A10) yields
\begin{equation} F = {\rm i}\left(\frac{1}{\zeta + 1}\sum_{m={-}\infty}^{\infty} \hat{a}_{0,m}\exp({\, j m \sigma t}) + \sum_{k = 1}^{\infty} \frac{1}{\zeta ^k} \sum_{m={-}\infty}^{\infty} \hat{a}_{k,m}\exp({\, j m \sigma t})\right), \end{equation}
which we can rewrite as a Fourier expansion of 
$F$
\begin{equation} F = \sum_{m ={-}\infty}^{\infty} \hat F_m \exp({\, jm\sigma t}) = \sum_{m={-}\infty}^{\infty} {\rm i}\left(\frac{\hat{a}_{0,m}}{\zeta + 1} + \sum_{k = 1}^{\infty}\frac{\hat{a}_{k,m}}{\zeta ^k} \right) \exp({\, jm\sigma t}). \end{equation}
Evaluating (A13) on the unit circle (
$\zeta = \textrm {e}^{\textrm {i} \theta }$), which corresponds to the surface of the plate in the 
$z$ plane, and separating the real and imaginary parts yields
$$\begin{gather} \phi|_W = \sum_{m={-}\infty}^{\infty} \left(\frac{1}{2}\hat{a}_{0,m} \tan \frac{\theta}{2} + \sum_{k = 1}^\infty \hat{a}_{k,m} \sin{k\theta} \right) \exp({\, j m \sigma t}), \end{gather}$$
$$\begin{gather}\psi|_W = \sum_{m={-}\infty}^{\infty} \left( \tfrac{1}{2} \hat a_{0,m} + \sum_{k = 1}^{\infty} \hat{a}_{k,m}\cos{k\theta}\right) \exp({\, j m \sigma t}) = \sum_{m={-}\infty}^{\infty} \varPsi_m(x) \exp({\, j m \sigma t}). \end{gather}$$
We recognize the expression for 
$\varPsi _m(x)$ as a Chebyshev series in 
$x$
\begin{equation} \varPsi_{m}(x) = \tfrac{1}{2}\hat{a}_{0,m} + \sum_{k = 1}^{\infty} \hat{a}_{k,m} T_{k} (x), \quad \forall\ m \in \mathbb{Z}, \end{equation}
where 
$T_k(x) = \cos (k \arccos x)$ is the 
$k$th Chebyshev polynomial, and 
$x = \cos \theta$.
The complex acceleration potential can be related to the complex velocity 
$w = u - \textrm {i}v$ through the momentum equation (2.4b), by
\begin{equation} \frac{\partial F}{\partial z} = \frac{\partial w}{\partial t} + \frac{\partial w}{\partial z}. \end{equation}
Evaluating the imaginary part on the plate's surface (
$z = x$) and substituting the no-penetration boundary condition (2.6a), yields
\begin{equation} \left.\frac{\partial \psi}{\partial x}\right|_{y=0} ={-}\left(\frac{\partial}{\partial t} + \frac{\partial}{\partial x} \right)^{2} Y. \end{equation}Substituting (A1a) and (A14b) into (A17) gives
\begin{equation} \sum_{m={-}\infty}^{\infty} D \varPsi_m (x) \exp({\, j m \sigma t}) ={-}\sum_{m={-}\infty}^{\infty} (\, j m \sigma + D)^{2}\hat{y}_m (x)\exp({\, j m \sigma t}), \end{equation}
where 
$D = {\textrm {d}}/{\textrm {d}\kern 0.06em x}$. Separating Fourier components gives
\begin{equation} D\varPsi_m (x) ={-}(\, j m \sigma + D)^{2} \hat{y}_m (x), \quad \forall\ m \in \mathbb{Z}. \end{equation}
Given 
$\hat y_m$, we expand it into a Chebyshev series and use (A19) to express the Chebyshev coefficients of 
$\varPsi _m$ – that is, 
$\hat a_{k,m}$ for 
$k \ge 1$ – in terms of the Chebyshev coefficients of 
$\hat y_m$. To determine 
$\hat {a}_{0,m}$, we expand the vertical velocity on the plate in a Fourier series
\begin{equation} v|_{W} = \sum_{m ={-}\infty}^{\infty} \hat{v}_m(x) \exp({\, j m \sigma t}), \end{equation}and the spatial coefficients in Chebyshev series
\begin{equation} \hat{v}_m(x) = \tfrac{1}{2} \hat{V}_{0,m} + \sum_{k = 1}^{\infty} \hat{V}_{k,m}T_k(x),\quad \forall\ m \in \mathbb{Z}. \end{equation}We can express the no-penetration boundary condition (2.6a), as
\begin{equation} \hat{v}_m(x) = (\, j m \sigma + {D})\hat{y}_m(x),\quad \forall\ m \in \mathbb{Z}. \end{equation}
Given the Chebyshev coefficients of 
$\hat y_m$, (A22) gives the Chebyshev coefficients of 
$\hat v_m$. The coefficient 
$\hat a_{0,m}$ is then given by
\begin{equation} \hat{a}_{0,m} ={-}C(\, j m\sigma)(\hat{V}_{0,m} + \hat{V}_{1,m}) + \hat{V}_{1,m}, \quad \forall\ m \in \mathbb{Z}, \end{equation}where
\begin{equation} C(\, j m \sigma) = \frac{K_1(\, jm\sigma)}{K_0(\, jm\sigma) + K_1(\, jm\sigma)} \end{equation}
is the Theodorsen function, and 
$K_\nu$ is the modified Bessel function of the second kind of order 
$\nu$. The expression for 
$\hat a_{0,m}$ is derived in Wu (Reference Wu1961).
With all 
$\hat a_{k,m}$ in hand, the hydrodynamic load is given by
\begin{equation} \hat q_m(x) = \hat a_{0,m} \sqrt{\frac{1-x}{1+x}} + 2 \sum_{k=1}^{\infty} \hat a_{k,m} \sin k\theta, \quad \forall\ m \in \mathbb{Z}. \end{equation}
The hydrodynamic load 
$\hat q_m$ depends linearly on 
$\hat y_m$.
To summarize, given the kinematics 
$\hat y_m$, we can calculate the coefficients 
$\hat a_{k,m}$, with which we can calculate the hydrodynamic load, which alters the kinematics via (A3). The coupled fluid–structure problem must be solved numerically. The pseudospectral numerical method for constant stiffness is given by Moore (Reference Moore2017), which is relatively straightforward to adapt to account for time-periodic stiffness; in it, the kinematics are expanded into Chebyshev series. All infinite series are truncated to finite series; for the results presented in this work, we have used 16 Fourier modes and 64 Chebyshev modes. The method is fast and accurate, pre-conditioning the system with continuous operators to avoid errors typically encountered when using spectral methods to solve high-order differential equations. Formulas to calculate the thrust and power coefficients are given in Moore (Reference Moore2017).
Appendix B. Floquet analysis
To test the stability of the solutions computed using the method in Appendix A, we perform a Floquet analysis of the problem with homogeneous boundary conditions. This analysis is adapted from the eigenvalue problem described by Floryan & Rowley (Reference Floryan and Rowley2018, Reference Floryan and Rowley2020). Following the preceding analysis, but not assuming a form for the time dependence, we arrive at
$$\begin{gather} 2R Y_{tt} + \tfrac{2}{3} S(t) Y_{xxxx} = {\rm \Delta} p, \end{gather}$$
$$\begin{gather}Y(x,t) = \tfrac{1}{2} y_0(t) + \sum_{k=1}^\infty y_k(t) T_k(x), \end{gather}$$
$$\begin{gather}{\rm \Delta} p(x,t) = a_0(t)\sqrt{\frac{1-x}{1+x}}+2\sum_{k=1}^{\infty}a_k(t)\sin{k \theta}, \end{gather}$$
$$\begin{gather}\sum_{k=1}^{\infty}a_k(t) T_k'(x) ={-}\tfrac{1}{2}\ddot{y}_0(t) - \sum_{k=1}^{\infty} [\ddot{y}_k(t) T_k (x) + 2 \dot{y}_k(x) T_k'(x) + y_k (t)T_k''(x)], \end{gather}$$
$$\begin{gather}Y({-}1,t) = 0, \quad Y_{x}({-}1,t) = 0, \quad Y_{xx}(1,t) = Y_{xxx}(1,t) = 0. \end{gather}$$
Above, a dot denotes differentiation with respect to 
$t$. In (B1b), we have written 
$Y$ as a Chebyshev series in 
$x$. The expression in (B1d) derives from the no-penetration condition (A17).
In what follows, we set 
$S(t) = \bar S + ({\bar S S_0}/{2})(\exp(j \omega t) + \exp ({- j \omega t}))$, as in (2.11c). Here, 
$\omega$ is the frequency of the stiffness oscillation. Assuming the Floquet solution form
$$\begin{gather} y_k(t) = {\rm e}^{\lambda t} \sum_{m ={-}\infty}^{\infty} \hat y_{k,m}\exp({\, j m \omega t}), \end{gather}$$
$$\begin{gather}a_k(t) = {\rm e}^{\lambda t} \sum_{m ={-}\infty}^{\infty} \hat a_{k,m} \exp({\, j m \omega t}),\quad \forall\ k \in \mathbb{W}, \end{gather}$$
where 
$\lambda$ is the Floquet exponent, gives, after separating the Fourier components, the following set of equations:
$$\begin{gather} 2R(\lambda + j m \omega)^2 \hat{\boldsymbol y}_m + \tfrac{2}{3} \bar{S} {{\boldsymbol{\mathsf{D}}}}^4 \hat{\boldsymbol y}_m + \tfrac{1}{3} \bar S S_0 {{\boldsymbol{\mathsf{D}}}}^4 \hat{\boldsymbol y}_{m-1} + \tfrac{1}{3} \bar S S_0 {{\boldsymbol{\mathsf{D}}}}^4 \hat{\boldsymbol y}_{m+1} = \hat{\boldsymbol p}_m, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol p}_m = {{\boldsymbol{\mathsf{A}}}} \hat{\boldsymbol a}_m, \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{D}}}} \hat{\boldsymbol a}_m ={-}(\lambda + j m \omega)^2 \hat{\boldsymbol y}_m - 2(\lambda + j m \omega) \boldsymbol{D} \hat{\boldsymbol y}_m - {{\boldsymbol{\mathsf{D}}}}^2 \hat{\boldsymbol y}_m, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol V}_m = (\lambda + j m \omega) \hat{\boldsymbol y}_m + {{\boldsymbol{\mathsf{D}}}}\hat{\boldsymbol y}_m, \end{gather}$$
$$\begin{gather}\hat a_{0,m} ={-}C(\lambda + j m \omega)(\hat V_{0,m} + \hat V_{1,m}) + \hat V_{1,m}, \quad \forall\ m \in \mathbb{Z}. \end{gather}$$
Above, 
$\hat {\boldsymbol y}_m$ is a vector of the Chebyshev coefficients corresponding to index 
$m$ in (B2a): 
$\hat {\boldsymbol y}_m = [\hat y_{0,m} \enspace \hat y_{1,m} \enspace \hat y_{2,m} \enspace \cdots ]^\textrm {T}$; analogous expressions hold for 
$\hat {\boldsymbol p}_m$ (pressure), 
$\hat {\boldsymbol a}_m$ (potential) and 
$\hat {\boldsymbol V}_m$ (vertical velocity). The equality in (B3b) states that the Chebyshev coefficients of the pressure are linear combinations of the Chebyshev coefficients of the potential, and 
$\boldsymbol {D}$ is the differentiation operator in Chebyshev space. Putting these equations together gives
\begin{align} & {{\boldsymbol{\mathsf{A}}}}[\boldsymbol{e_1}(\boldsymbol{e_2}-C(\lambda_m) \boldsymbol{e_1}-C(\lambda_m)\boldsymbol{e_2})^{\rm T}(\lambda_m \boldsymbol{\hat y}_m + {{\boldsymbol{\mathsf{D}}}} \boldsymbol{\hat y}_m) - \lambda_m ^2 \boldsymbol{D}^- \boldsymbol{\hat y}_m - 2 \lambda_m {{\boldsymbol{\mathsf{D}}}}^- {{\boldsymbol{\mathsf{D}}}}\boldsymbol{\hat y}_m \nonumber\\ &\quad - {{\boldsymbol{\mathsf{D}}}}^-{{\boldsymbol{\mathsf{D}}}}^2 \boldsymbol{\hat y}_m] - 2R\lambda_m^2 \boldsymbol{\hat y}_m - \frac{2}{3}{{\boldsymbol{\mathsf{D}}}}^4 \bar{S} \boldsymbol{\hat y}_m = \frac{\bar S S_0}{3}{{\boldsymbol{\mathsf{D}}}}^4 [ \boldsymbol{\hat y}_{m-1} + \boldsymbol{\hat y}_{m+1}], \quad \forall\ m \in \mathbb{Z}, \end{align}
where 
$\lambda _m = \lambda + j m \omega$, 
$\boldsymbol{\mathsf{D}}^-$ is the Chebyshev-space representation of the integration operator that makes the first Chebyshev coefficient zero, and 
$\boldsymbol {e}_k$ is the 
$k$th Euclidean basis vector.
Because of the presence of the Theodorsen function, solving for the Floquet exponents requires solving a nonlinear generalized eigenvalue problem. Since we are mainly interested in delineating regions of parameter space where solutions are unstable, we instead follow the idea of Kumar & Tuckerman (Reference Kumar and Tuckerman1994) to find the marginal stability curves. Rather than calculating the Floquet exponents for given values of the parameters, we set the value of the Floquet exponent such that 
$Re(\lambda ) = 0$ and solve for 
$S_0$. Physically, we are trying to find the strength of the stiffness oscillation that borders regions of stability and instability. This leads to a linear generalized eigenvalue problem where 
$S_0$ is the eigenvalue. Note that we must also set a value for 
$\textrm {Im}(\lambda ) \in [0, \omega ]$; 
$\textrm {Im}(\lambda ) = 0$ is called the harmonic case, and 
$Im(\lambda ) = \omega /2$ is called the subharmonic case (Kumar & Tuckerman Reference Kumar and Tuckerman1994).
As can be seen in (B4), the time dependence of the stiffness causes cross-frequency coupling in the kinematics. We write (B4) compactly as
\begin{equation} {{\boldsymbol{\mathsf{B}}}}_m \boldsymbol{\hat y}_{m} = \frac{\bar S S_0}{3} {{\boldsymbol{\mathsf{C}}}}_m\boldsymbol{\hat y}_{m}, \quad \forall\ m \in \mathbb{Z}, \end{equation}where
\begin{align} {{\boldsymbol{\mathsf{B}}}}_m &= \boldsymbol{\mathsf{A}}[\boldsymbol{e_1} (\boldsymbol{e_2}-C(\lambda_m)\boldsymbol{e_1}-C(\lambda_m)\boldsymbol{e_2})^{{\rm T}} (\lambda_m {{\boldsymbol{\mathsf{I}}}} + {{\boldsymbol{\mathsf{D}}}}) - \lambda_m ^2 {{\boldsymbol{\mathsf{D}}}}^-- 2 \lambda_m {{\boldsymbol{\mathsf{D}}}}^- {{\boldsymbol{\mathsf{D}}}} - {{\boldsymbol{\mathsf{D}}}}^-{{\boldsymbol{\mathsf{D}}}}^2 \nonumber\\ &\quad -2R \lambda_m^2 {{\boldsymbol{\mathsf{I}}}} - \tfrac{2}{3}{{\boldsymbol{\mathsf{D}}}}^4 \bar{S}, \quad \forall\ m \in \mathbb{Z}. \end{align}
The operator 
$\boldsymbol{\mathsf{C}}_m$ maps 
$\boldsymbol {\hat y}_m$ to 
$\boldsymbol{\mathsf{D}}^4 (\boldsymbol {\hat y}_{m+1} +\boldsymbol {\hat y}_{m-1})$. This is a linear generalized eigenvalue problem where the eigenvalues are the amplitudes of the stiffness oscillation, 
$S_0$.
To proceed, we truncate all Chebyshev expansions to the upper limit 
$N$, and truncate all Fourier expansions so that they contain frequencies up to 
$M\omega$. Doing so makes 
$\boldsymbol {\hat y}_m$ a vector of length 
$N+1$. We must simultaneously solve the system (B5) for all 
$\boldsymbol {\hat y}_m$, leading to a system of size 
$2M(N+1) \times 2M(N+1)$. We choose 
$N = 16$ and 
$M = 11$. We reduce the order of the Chebyshev and Fourier modes so the Floquet problem can be solved in reasonable time.
In the harmonic (
$\textrm {Im}(\lambda ) = 0$) and subharmonic (
$\textrm {Im}(\lambda ) = \omega$) cases, 
$\boldsymbol {\hat y}_m$ must obey the reality conditions: 
$\boldsymbol {\hat y}_{-m} = \boldsymbol {\hat y}_m^*$ in the harmonic case, and 
$\boldsymbol {\hat y}_{-m} = \boldsymbol {\hat y}_{m-1}^*$ in the subharmonic case. The reality conditions allow us to rewrite the Fourier expansions in terms of only non-negative indices. For 
$0 < \textrm {Im}(\lambda ) < \omega$, positive and negative frequencies are independent of each other, and (B2a) must be added to its complex conjugate to form a real field.
Explicitly writing the real and imaginary components in (B5) yields
\begin{equation} ({{\boldsymbol{\mathsf{B}}}}_m^r + {\rm i}{{\boldsymbol{\mathsf{B}}}}_m^i) (\boldsymbol{\hat y}_m^r + {\rm i}\boldsymbol{\hat y}_m^i) = \frac{\bar S S_0}{3} {{\boldsymbol{\mathsf{D}}}}^4 ( \boldsymbol{\hat y}_{m-1}^r + {\rm i}\boldsymbol{\hat y}_{m-1}^i + \boldsymbol{\hat y}_{m+1}^r + {\rm i}\boldsymbol{\hat y}_{m+1}^i), \quad \forall\ m \in \mathbb{Z}, \end{equation}from which we can separate the real and imaginary components to get
$$\begin{gather} {{\boldsymbol{\mathsf{B}}}}_m^r \boldsymbol{\hat y}_m^r - {{\boldsymbol{\mathsf{B}}}}_m^i \boldsymbol{\hat y}_m^i = \frac{\bar S S_0}{3} {{\boldsymbol{\mathsf{D}}}}^4 (\boldsymbol{\hat y}_{m-1}^r + \boldsymbol{\hat y}_{m+1}^r), \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{B}}}}_m^r \boldsymbol{\hat y}_m^i + {{\boldsymbol{\mathsf{B}}}}_m^i \boldsymbol{\hat y}_m^r = \frac{\bar S S_0}{3} {{\boldsymbol{\mathsf{D}}}}^4 (\boldsymbol{\hat y}_{m-1}^i + \boldsymbol{\hat y}_{m+1}^i), \quad \forall\ m \in \mathbb{Z}. \end{gather}$$This can be rewritten in matrix form as
\begin{equation}
\begin{bmatrix} {{\boldsymbol{\mathsf{B}}}}_m^r & -
{{\boldsymbol{\mathsf{B}}}}_m^i \\
{{\boldsymbol{\mathsf{B}}}}_m^i &
{{\boldsymbol{\mathsf{B}}}}_m^r \end{bmatrix}
\begin{bmatrix} \boldsymbol{\hat y}_m^r \\ \boldsymbol{\hat
y}_m^i \end{bmatrix} = \frac{\bar S S_0}{3} \begin{bmatrix}
{{\boldsymbol{\mathsf{C}}}}_m & \boldsymbol{0} \\
\boldsymbol{0} & {{\boldsymbol{\mathsf{C}}}}_m
\end{bmatrix} \begin{bmatrix} \boldsymbol{\hat y}_m^r \\
\boldsymbol{\hat y}_m^i \end{bmatrix}, \quad \forall\ m \in
\mathbb{Z},
\end{equation}
where 
$\boldsymbol{\mathsf{C}}_m$ is as before.
Finally, we incorporate the boundary conditions into (B9). The formula to evaluate a Chebyshev series at the endpoints is
\begin{equation} \hat y_m({\pm} 1) = \tfrac{1}{2} {\hat y}_{0,m} + \sum_{k = 1} ^N ({\pm} 1)^k {\hat y}_{k,m}, \end{equation}which can be split into its real and imaginary parts
$$\begin{gather} \hat y_m^r({\pm} 1) = \tfrac{1}{2} {\hat y}_{0,m}^r +\sum_{k = 1} ^N ({\pm} 1)^k {\hat y}_{k,m}^r, \end{gather}$$
$$\begin{gather}\hat y_m^i({\pm} 1) = \tfrac{1}{2} {\hat y}_{0,m}^i +\sum_{k = 1} ^N ({\pm} 1)^k {\hat y}_{k,m}^i. \end{gather}$$
This allows us to express the boundary conditions in (B1e) in terms of the Chebyshev coefficients. To enforce the boundary conditions, we replace the last four rows of the first block in (B9) by the four boundary conditions on the real part, and the last four rows in (B9) by the four boundary conditions on the imaginary part. Combining the equations for all values of 
$m$ into one large system leads to a generalized eigenvalue problem of the form 
$(3/\bar S)\boldsymbol{\mathsf{B}} \boldsymbol {\hat y} = S_0 \boldsymbol{\mathsf{C}}_H \boldsymbol {\hat y}$, which is linear in 
$S_0$. In the harmonic case, the matrices 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{C}}_H$ take the form
$$\begin{gather}
{{\boldsymbol{\mathsf{B}}}} = \begin{bmatrix}
{{\boldsymbol{\mathsf{B}}}}_0^r & -
{{\boldsymbol{\mathsf{B}}}}_0^i & \boldsymbol 0 &
\boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 & \cdots \\
{{\boldsymbol{\mathsf{B}}}}_0^i &
{{\boldsymbol{\mathsf{B}}}}_0^r & \boldsymbol 0 &
\boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 & \cdots \\
\boldsymbol 0 & \boldsymbol 0 &
{{\boldsymbol{\mathsf{B}}}}_1^r & -
{{\boldsymbol{\mathsf{B}}}}_1^i & \boldsymbol 0 &
\boldsymbol 0 & \cdots \\
\boldsymbol 0 & \boldsymbol 0 &
{{\boldsymbol{\mathsf{B}}}}_1^i &
{{\boldsymbol{\mathsf{B}}}}_1^r & \boldsymbol 0 &
\boldsymbol 0 & \cdots \\
\boldsymbol 0 & \boldsymbol 0 &
\boldsymbol 0 & \boldsymbol 0 &
{{\boldsymbol{\mathsf{B}}}}_2^r & -
{{\boldsymbol{\mathsf{B}}}}_2^i & \cdots \\
\boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 &
{{\boldsymbol{\mathsf{B}}}}_2^i &
{{\boldsymbol{\mathsf{B}}}}_2^r & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} &
\mathbf{0} & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots
\end{bmatrix},
\end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{C}}}}_H
= \begin{bmatrix} \mathbf{0} & \mathbf{0} & 2
{{\boldsymbol{\mathsf{D}}}}^4 & \mathbf{0} & \mathbf{0} &
\mathbf{0} & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0}
& \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\
{{\boldsymbol{\mathsf{D}}}}^4 & \mathbf{0} & \mathbf{0} &
\mathbf{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \mathbf{0} & \cdots \\
\mathbf{0} & {{\boldsymbol{\mathsf{D}}}}^4 &
\mathbf{0} & \mathbf{0} & \mathbf{0} &
{{\boldsymbol{\mathsf{D}}}}^4 & \cdots \\
\mathbf{0} & \mathbf{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \mathbf{0} &
\mathbf{0} & \mathbf{0} & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \mathbf{0}
& \mathbf{0} & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & {{\boldsymbol{\mathsf{D}}}}^4 &
\mathbf{0} & \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix},
\end{gather}$$with the boundary conditions incorporated as previously described, and the Chebyshev coefficients for all frequencies are stored in the vector
\begin{equation} \boldsymbol{\hat y} = [(\boldsymbol{\hat y}_0^r)^{\rm T} \quad (\boldsymbol{\hat y}_0^i)^{\rm T} \quad (\boldsymbol{\hat y}_1^r)^{\rm T} \quad (\boldsymbol{\hat y}_1^i)^{\rm T} \quad \cdots \quad (\boldsymbol{\hat y}_M^r)^{\rm T} \quad (\boldsymbol{\hat y}_M^i)^{\rm T}]^{\rm T}. \end{equation} For the subharmonic case, 
$\boldsymbol{\mathsf{B}}$ is identical to the harmonic case, but by using the reality condition 
$\boldsymbol {\hat y}_{-m} = \boldsymbol {\hat y}_{m-1}^*$, 
$\boldsymbol{\mathsf{C}}_{SH}$ takes the form
\begin{equation} {{\boldsymbol{\mathsf{C}}}}_{SH} = \begin{bmatrix} {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \cdots \\ \boldsymbol{0} & -{{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \cdots \\ {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \cdots \\ \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \cdots \\ \boldsymbol{0} & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \cdots \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \boldsymbol{0} & \cdots \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & {{\boldsymbol{\mathsf{D}}}}^4 & \boldsymbol{0} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}. \end{equation}Appendix C. Modelling verification and validation
There are two aspects to consider: validation and verification. Validation answers the question: ‘Are we solving the correct equations?’ That is, it answers whether the equations we solve reflect physical reality. Verification asks: ‘Are we solving those equations correctly?’
First, we want to validate the physical model. To model stiffness variation, we take the Young's modulus of an Euler–Bernoulli beam to be a function of time while leaving other properties constant. The form of the partial differential equation (PDE) governing the deflection of the beam is unchanged from the Euler–Bernoulli equation, but the term corresponding to resistance to bending now has a time-varying coefficient due to the time-varying Young's modulus. To see that the form does not change, we start by forming the Lagrangian of the system
\begin{equation} \mathcal{L} = \tfrac{1}{2} \rho_s d w Y_t^2 - \tfrac{1}{2}EIY_{xx}^2 + w{\rm \Delta} p Y. \end{equation}The first term in the Lagrangian corresponds to kinetic energy, the second term corresponds to strain energy (analogous to the potential energy stored by a stretched spring) and the third term corresponds to work done by the external loading. Above, the Young's modulus is allowed to be a function of time. The corresponding Euler–Lagrange equation is
\begin{equation} \frac{\partial \mathcal{L}}{\partial Y} - \frac{\partial}{\partial t} \left(\frac{\partial \mathcal{L}}{\partial Y_t}\right) + \frac{\partial^2}{\partial x^2}\left(\frac{\partial \mathcal{L}}{\partial Y_{xx}}\right) = 0. \end{equation}The Euler–Lagrange equation yields the same PDE as for a beam with constant properties. This is because the term in the Lagrangian involving the time-varying Young's modulus does not enter the term in the Euler–Lagrange equation that involves differentiation with respect to time. This PDE was derived from the same Lagrangian used to derive the classical Euler–Bernoulli beam equation, whose validity is well established. Thus, our model for time-varying stiffness is equally valid and subject to the same limitations as the classical Euler–Bernoulli beam equation.
With the validity of our equations established, we next address verification (i.e. whether our numerical method is correctly solving the equations). This can be done in several ways. Firstly, we check whether our solution methodology reproduces previously published results We modify the fast Chebyshev method introduced by Moore (Reference Moore2017) to account for time-periodic stiffness. Due to the Fourier series expansion method and the time-periodic stiffness used here we have cross-coupling of Fourier modes. When the plate has constant stiffness in time there is no coupling of modes; in fact, only the frequency that the plate is actuated at is active. We verify that the results produced from our code with constant stiffness are identical to those produced by Moore (Reference Moore2017). These plots are seen in figure 13.
Reproduction of figure 7 in Moore (Reference Moore2017) using our solution method: (a) 
$C_T$; (b) 
$C_P$; (c) 
$\eta = C_P/C_T$.
Secondly, we check whether our numerical solutions indeed satisfy our equations; that is, starting from a computed solution, we calculate each term that appears in the governing equation and check whether all the terms balance. We have done so for numerous cases and have found that our solutions indeed satisfy the governing equations.
Lastly, we check our numerical method against a different numerical method. Specifically, we check whether solutions obtained via our harmonic balance approach match those obtained via time stepping. This check cannot be performed on the fluid–structure interaction problem we study here (because the modelling of the wake is performed in the frequency domain, and there is no analogous time-domain wake model that has been published); however, the numerical method we use (and its implementation in code) is generic, so we can verify its correctness by checking it against any problem. We have considered an analogous problem of a mass–spring–damper system with a time-varying spring constant, solving for the dynamics using our harmonic balance code and a time stepper. After transients decay, the numerical solutions obtained using the two different numerical methods agree.
Open access
