Fluid–structure interaction analysis on motion control of a self-propelled flexible plate near a rigid body utilizing PD control

Inspired by a previous experimental study of fish swimming near a cylinder, we numerically investigate the swimming and station-holding behavior of a flexible plate ahead of a circular cylinder whose motion is controlled by a proportional–derivative (PD) controller. Specifically, the deformation of this two-dimensional plate is actuated by a periodically varying external force applied on the body surface, which mimics the fish muscle force to produce propulsive thrust. The actuation force amplitude is dynamically adjusted by a feedback controller to instruct the plate to swim the desired distance from an initial position to a target location and then hold the station there. Instead of directly using the instantaneous position signal, an average speed measured over one force actuation period is proposed with the inclusion of instantaneous position information to form the tracking error for the PD control. Our results show that the motion control of swimming and station holding has been achieved by this simple but effective feedback control without large overshoot when approaching the target at different flow conditions and actuation force formulas. Although the swimming distance remains the same, a plate whose initial position is closer to the cylinder requires less energy expenditure to swim to the target location and hold the station there. This is because the low-pressure zone near the trailing edge of the plate is reduced in size, which provides drag reduction, contributing to reduced swimming energy. A higher Reynolds number also leads to energy savings. Under the same control strategy, the swimming performance is more affected by the force-frequency while the phase shift of the actuation force has a less significant impact.


Introduction
Aquatic animals provide ample prototypes for the design of efficient artificial underwater robots and vehicles. This motivates biologists, physicists, and engineers to study the mechanism of fish swimming through experimental and theoretical modeling. Most of the existing research focuses on the scenarios when fish swims in the still fluid or uniform flow (Liu et al 2017, Maertens et al 2017, Mittal et al 2006. For example, straightforward swimming and turning maneuvers of a 3D eel-like swimmer in turbulent flow were simulated by Leroyer and Visonneau (2005) which provided quantities that can hardly be measured on live species, like drag and efficiency. Recently published papers reported pioneering applications of large eddy simulations to study more fish-like selfpropelled virtual swimmers in turbulent flow conditions (Bottom Ii et al 2016, Borazjani 2015, Ogunka et al 2020). These studies improve our understanding of fundamental hydrodynamics during the locomotion of swimming fish. However, the majority of fishes commonly experience altered flows, a scenario different from the still or steady flow conditions (Liao 2007). Fishes present different swimming behaviors in these altered flow conditions in contrast to still or uniform flow. For instance, biological observations by Shuler et al (1994) revealed that brown trouts in rivers prefer a feeding site that is near wing dams, midchannel boulder clusters, and natural bank cover, but avoid areas without structures. Another example is that dolphins are seen to ride the bow wake wave behind ships (Scholander 1959).
In addition to the above biological observations in the field, laboratory studies, in which the fluid conditions are subject to elaborated control, have also been conducted to explore the swimming behaviors and underlying hydrodynamic mechanism. One example is the experimental study of rainbow trout swimming near a bow cylinder in the uniform flow by Liao et al (2003). They compared the axial swimming kinematics in response to three different fluid conditions, i.e. in the Karman vortex street behind the cylinder, in the cylinder bow wake, and the free stream. The trout fish behind the cylinder was found to adopt a distinctive gait pattern to hold the station, which was termed the Karman gait. Their measurement revealed that the body oscillation amplitude and curvatures are much larger than that when the fish swims in the absence of the upstream cylinder. Besides, the fish also voluntarily matches its tail-beat frequency with that of the vortex shedding. This could be the evidence that fish can actively capture energy from the vortices produced by the environment during the station-holding. Another phenomenon discovered by this experiment is that the fish tends to hold the station in the bow wake (in front) of the cylinder rather than behind the cylinder. Following this study, Beal et al (2006) found that even a dead fish can be propelled upstream when placed in the wake of a bluff cylinder, for which the energy harvest is purely passive due to the resonance oscillation of its flexible body under the action of incoming vortices. Inspired by the above experimental studies of live fish, some numerical simulations were conducted using computational fluid dynamics (CFD) techniques in an attempt to elucidate the underlying hydrodynamic mechanism when fish swims near a rigid body, e.g. a cylinder. It started with the simplification that the model kinematics was prescribed mathematically in many studies (Bao and Tao 2013, Wu et al 2014, Xiao et al 2012b, in which the dynamic interaction between the flexible deformation and the surrounding flow is excluded. This is distinct from real biological scenarios where the fish body and fins are characterized as flexible structures. A few studies considered the fluid-structure interaction (FSI) when the swimmer is swimming in the altered flow. For example, by assuming the fish body as a flexible filament, Tian et al (2011) studied the hydrodynamic force and the deformation pattern of the filament in the Karman gait region and the entrainment region termed by Liao et al (2003). A subsequent study reported the swimming behavior of a self-propelled flexible fin in the wake of a circular cylinder . They observed three patterns of the plate in the cylinder wake: propulsion upstream, drift downstream (DD), and holding station at an equilibrium position. The station-holding mode is dependent on the initial longitudinal position of the fin relative to the cylinder and its heave frequency. Similar phenomenon was reported in a following study (Wang et al 2019), in which two tandem cylinders were considered to yield more complex hydrodynamic environment.
Nevertheless, fishes were found to favor holding the station in the bow wake (Liao et al 2003). Most of the existing research concentrated on the swimming behind the cylinder, while few on the behavior in front of the cylinder. It is thus necessary to study the latter scenario which would help to elucidate the mechanism of station-holding or the socalled 'flow refuging' behavior. One such example is an FSI simulation of the interaction between a flexible filament and a downstream rigid body by Tian et al (2010). This model was fixed at a location without self-propulsion, thus it is unclear how the interaction would influence free-swimming behavior. Carling et al (1998) suggested that the force calculated in this tethered mode could differ remarkably from the force developed by a freely swimming model. Although their work demonstrated that drag reduction is possible as the fluid immediately upstream of a bluff body is being 'pushed', the deformation of the filament was purely passive, distinct to the fish swimming where both active actuation and passive deformation are involved to achieve the objectives of swimming forward or holding station.
Indeed, very limited research considered dynamically changing the swimming gait or actuation scheme, as real fishes do. A real fish can constantly and actively tune its body and/or fin oscillation for acceleration/deceleration to swim towards a specific target, as well as holding a position in an unsteady fluid environment. A recent study reported the effective implementation of a feedback proportional-integral-derivative (PID) controller to manipulate a soft plate to follow specific target trajectories (Hess et al 2020). In their work, the swimming speed and direction were dynamically controlled by actively changing the muscle-like contractive actuation strain. It demonstrated the feasibility to employ a simple control to fish swimming study.
In this work, we numerically study the swimming performance and station-holding capability of a 2D flexible plate in front of a fixed cylinder in a uniform flow. The flexible plate with uniform stiffness along the centerline is actuated via an external force, which is periodically and uniformly applied along the lateral direction. In the following, this force is referred to as actuation force. The amplitude of the actuation force (control variable) is dynamically changed by a feedback proportional-derivative (PD) controller that neglects the integral (I) term of a typical PID control in which all the proportional, integral, and derivative terms of the tracking error are considered to calculate and reduce the difference between the desired setpoint and the measured process variable. Similar to real fish swimming scenarios (Liao et al 2003), the control objective is to instruct the plate to swim a certain distance to a target position and hold the station there, subject to different initial locations relative to the cylinder. Through the coupling of PD control with the FSI numerical solver, we aim to investigate the effectiveness of real-time feedback control in achieving swimming and station-holding involving complex FSI. We also aim at exploring energy consumption under controlled actuation, which cannot be measured in a biological experimental environment.
The remainder of this paper is organized as follows. The geometry, actuation mechanism, fluid environment, and control scheme of the self-propelled plate are introduced in section 2. The metrics to evaluate the swimming and station-holding performance are also defined in this section. In the next section, the numerical models and approaches to discretize and solve the governing equations are described, followed by the implementation of the PD control and selfconsistency study. Section 4 presents the numerical results including swimming performance, flow field, and the structural deformation under the closed-loop PD control. Finally, the conclusions are drawn in section 5.

Problem statement
The model considered in this work is a twodimensional flexible plate, placed in front of a fixed cylinder, as shown in figure 1. This 2D modeling is chosen to reduce computational cost because of numerous feedback control loops that serve to reduce tracking errors. The initial distance from the plate to the cylinder is denoted by d 0 . Different from most of the previous studies (Dai et al 2016, the deformation of the plate is considered to be manipulated by an external force applied in the transverse direction, which is uniformly distributed on its surface that mimics the muscle forces of the fish fin. This actuation mechanism by adding an external force on the body surface has been successfully implemented for a squid-like jet model (Luo et al 2020b) and a ribbon tail ray-like bio-inspired underwater robot (Shi and Xiao 2021).
In structural dynamics, the leading edge of the plate is fixed (pinned) while the trailing edge is set as free. Following previous similar studies of selfpropelled swimmers near a cylinder , Shao et al 2010, Wang et al 2019, the plate is restricted to move in the horizontal direction because the effect of the presence of the cylinder on the plate is dominant by the horizontal distance between them. Besides, our numerical tests showed that the lateral movement is not substantial under the present actuation method and control scheme. Thus, the plate is only allowed for free-swimming in the horizontal direction. The induced horizontal motion along the x-direction is given by Newton's second law where m b is the mass of the plate, u b is the swimming velocity in x-direction, and F x represents the component of the overall fluidic force in the x-direction. A feedback controller is implemented to dynamically adjust the actuation force and therefore, control the motion of the plate. The control scheme, detailed in section 3.2, is derived to lead the plate to swim from an initial position x i = p 0 (measured at the leading edge of the plate) to a target position x t = p 0 − L (see in figure 1), and stay there, balancing thrust and drag forces over one plate deflection cycle period. As a result, the actuation force amplitude is given by where α denotes the control variable ranging between [0, 1] which serves as the adjustment factor of the actuation force amplitude, f ef is the actuation frequency of the force, and ϕ represents the phase shift of the actuation force. F ef0 is the maximum amplitude of the actuation force, given by F ef0 = 0.5ρ f U 2 L · C ef0 , where ρ f is the fluid density, U is the free stream velocity, L is the chord length of the plate which is equal to the diameter of the cylinder D (L = D), and C ef0 denotes the actuation force coefficient. The adjustment factor α is determined by the feedback controller and its calculation will be described in the next section. The Reynolds number is defined by where ν is the kinematic viscosity of the fluid. The instantaneous thrust coefficient C T is given by To evaluate the swimming performance at different initial positions d 0 , the time T s , which the plate takes to travel the distance of L to reach the target location, is recorded. The overall energy expenditure coefficient C Ps required by the plate during T s is defined as where P in (t) stands for the instantaneous power input as P in (t) = Γ f F ef (s, t) · u (s, t) ds with u (s, t) being the velocity of the grid node s where the actuation force is applied along the body surface Γ f . Besides, the mean energy expenditure coefficient required to dynamically maintain near the target point is calculated by where T h is the selected reference period when the plate stays around the target after reaching the position x t = x i − L, which is a constant for all the cases considered. Through the evaluation of equations (5) and (6), we can assess energy consumption during swimming and station holding. Key dimensionless parameters are as follows: the mass ratio m * = ρ s h/ρ f L with ρ s denoting the density of the plate and h being the thickness of the body; the non-dimensional actuation force frequency f * = f ef L/U; the dimensionless stiffness K = EI/ ρ f U 2 L 3 , where E is Young's modulus and I = h 3 /12 denotes the area moment of inertial of the cross-section; and the Poisson ratio ν s .

Numerical model and approach
The FSI numerical solver involves the coupled resolution of fluid and structural dynamics. The coupling between the fluid and structural solver is implemented in a partitioned manner. In this section, the governing equations and numerical methods of the fluid and structural dynamics are described, respectively, followed by a brief introduction of the coupling between them.
The fluid around the plate is governed by unsteady, compressible, and viscous Navier-Stokes equations. Under the assumption of no internal source term and by ignoring the body force, the governing equations can be written in an integral form as where G = (ρ f , ρ f v, ρ f E) T is the conservative variable vector, v is the velocity vector in the Cartesian coordinate system, and E denotes the total energy. Ω f represents the fluid control domain with the boundary as Γ f . F is the flex vector consisting of the convective and pressure fluxed fluxes arising from the viscous shear stress and thermal diffusion. The above governing equations are discretized and numerically resolved based on the cell-centered finite volume method. Specifically, with a structured grid method, the fluid domain Ω f is divided into an array of hexahedral grid cells, and each of these cells is indexed by i, j, k . The conservation laws are then applied to these cells, and their semi-discretized form for a cell can be obtained as where ΔΩ f is the control volume of the cell i, j, k , and H i,j,k is the total convective and diffusive fluxes entering the hexahedral cell through its surfaces. An artificial viscosity term D i,j,k is introduced to stabilize the scheme and eliminate the spurious numerical oscillations (Jameson et al 1981).
For the present unsteady simulations, we employ a dual-time stepping algorithm for temporal integration, in which equation (8) can be reformulated as a steady-state problem by introducing a pseudo-timet where where the superscripts n − 1 and n denote the two previous time steps, yielding a second-order accuracy. A hybrid multistage Runge-Kutta scheme is implemented to integrate equation (9). Parallel computation is achieved by domain decomposition based on message passing interface to enable large-scale computation in this numerical flow solver. The local time-stepping and multigrid methods are implemented to increase convergence speed. More details of this fluid solver are described in Sadeghi (2004). Besides, a multi-block overset grid system is utilized in this study to handle the relative motion between the plate and cylinder. This overset grid method is based on the implicit hole cutting technique proposed by Lee and Baeder (2003) and improved by Liao et al (2007). Specifically, there are two clusters of structured grids in our simulations, as shown in figure 4, one is the background grid containing the cylinder (cluster 1), and the other one is the body-fitted grid around the body (cluster 2). In this way, the latter cluster can move relative to the former one by avoiding large mesh distortion or even the occurrence of negative cell volume. The details of this overset grid assembler are described in Shi et al (2019). The grid deformation of the inner points within cluster 2 is obtained using a spring-analogy method (Batina 1990) and arc-lengthbased trans-finite interpolation method (Sadeghi et al 2003). To retain the original grid quality after deformation, Hermite polynomials as blending functions are used to maintain the grid angles of the original grid near the wall by specifying the derivatives of the grid coordinates at the boundaries. The details of the mesh smoothing algorithm can be found in Tsai et al (2001).
The self-propulsion of the plate in the horizontal direction within the computational domain is solved based on Newton's second law as mentioned in section 2. To determine the instantaneous lateral location and velocity, an explicit time marching scheme is implemented to discretize equation (1) as where n and n − 1 denotes the current and the previous time level, and Δt is the time step. The location of the plate is then obtained by integrating equation (11). It is necessary to ensure that the compressibility is negligibly small to simulate an incompressible flow problem using the present compressible flow solver. As we did in the previous applications to incompressible simulations, the freestream Mach number, defined as Ma ∞ = U/a ∞ with a ∞ denoting sound speed, is chosen as 0.06. This value is chosen to be far below the critical value of 0.3, where the compressibility effect becomes pronounced but still sufficiently large to ensure numerical stability. Given that the boundary is moving during the deformation of the plate, the actual Mach number may become larger than Ma ∞ near the body. Therefore, for computation accuracy, the local Mach numbers in the whole computational domain are monitored during computation to ensure that they are always below the critical value. This compressible flow solver has been successfully applied to various incompressible flow simulations involving flexible deformation in our previous work (Liu et al 2016, Luo et al 2020a, Luo et al 2020b, Xiao and Liao 2010, Xiao et al 2012a. This numerical solver is validated in Luo et al (2020a) and Luo et al (2020b). The quantitative validations include simulations of flow over a flexible plate attached behind a stationary square cylinder, the bending of a flexible plate in the uniform flow, the response of a flexible plate in forced heave motion (Luo et al 2020a), and flow over a NACA 0012 airfoil (Luo et al 2020b). Good agreements between the simulation results and the counterparts reported in the literature are obtained. The capability of this solver to resolve the self-propulsion of a flexible plate is further validated in appendix.
Regarding the structural dynamics, the deformation of the plate is governed by the weak form of momentum balance, which can be written in the differential form as where the left-hand side of the equation denotes the acceleration of a material point, calculated by the second derivative of the displacement vector S of the structure. P is the second Piola-Kirchoff stress tensor which models the surface forces, and the body forces per unit mass are represented by f. A constitutive equation that describes the correlation between the stress and the strain is established to close up equation (12), as follows. Specifically, for a Saint Venant-Kirchhoff material, the second Piola-Kirchoff stress tensor P is given by where C is the elastic tensor, E denotes the Green-Lagrange strain tensor, J represents the deformation gradient and δ is the unit tensor. The general governing equations of structural dynamics are discretized based on the finite element method (Dhondt 2004). With the application of the virtual work method, a linear algebraic equation system as a result of the discretization in the whole solid domain is given as  (2004) is implemented for temporal discretization. Denoting the velocity vector as {U} := Ṡ and the acceleration vector as {A} := S , the solution at the time step n + 1 can be obtained by where the predictor of velocity and displacement vector at time step n + 1 is denoted by Û n+1 and Ŝ n+1 , which are only dependent on the values at the time step n. γ and β are constants determined by a chosen constant α s . The structural solver is based on CalculiX and more details of the numerical methods can be found in Dhondt (2004).
Another main ingredient of the present FSI solver is the coupling between the fluid solver and the structural solver. In this work, the two solvers are coupled based on preCICE (Bungartz et al 2016) in a partitioned approach. An implicit coupling is implemented, in which sub-iterations and convergence criteria are introduced during each time step to ensure numerical stability and accuracy. An improved IQN-ILS method is also applied to stabilize and accelerate coupling iterations . The data mapping at the interface between the fluid and structural solver is established via the radial basis functions based interpolation. More details of this FSI solver are provided in our previous work (Luo et al 2020a).

Feedback controller
Control objectives are set for the swimming motion, i.e. swimming towards a target point from an initial point in front of the cylinder, and remaining at that target position. A feedback control system is developed for this purpose. Previous numerical studies have shown that PID control is robust and effective for fish swimming motion control Triantafyllou 2018, Hess et al 2020). A general form of a PD controller is given as where u (t) is the control action corresponding to the adjustment factor α of the force amplitude here, e (t) is the tracking error, and k p and k d are the tuning gains. The plate's body deformation and global motion are actively influenced by the actuation force. To achieve the control objectives, a control variable α is added to dynamically adjust the actuation force amplitude, as shown in equation (2). The feedback control diagram with PD control and the FSI solver is shown in figure 2.
In a previous study on swimming motion control, the tracking error is calculated using the sign from the dot-product of the swimming velocity vector and the distance vector (Hess et al 2020). This can avoid the overshoot that the plate moves past the target position without stopping at the right time. Different from the above method, we propose to use the swimming speed for control, rather than the position tracking error used in many previous studies (Kopman et al 2014, Yen et al 2018. This is because our numerical tests show that a tracking error that only considers the distance between the target position and the current position leads to a large overshoot in position control. Specifically, a PD controller is used to dynamically adjust the control variable α in order to drive the plate's velocity towards the desired velocity u set (t). u set is dependent on the gap between the current plate position and the target position, given by G (t) = x s (t) − x t with x s (t) denoting the instantaneous position of the plate's leading-edge and x t being the target position, as shown in figure 3. Given the periodic deformation of the oscillating plate, we use the mean velocity of the plate averaged over one actuation period in the tracking error evaluation instead of instantaneous value. This enhances numerical stability by avoiding high fluctuation of the error signal, which is essential to maintain stable control of the plate through balancing the thrust and drag forces over one cycle period, especially when approaching the target position.
Following the above logic, the tracking error e (t) is defined as where u b (t) is the time-average swimming velocity of the plate during one force actuation period T given by and the set velocity u set associated with the gap G is defined as where T set denotes a constant period used to tune the magnitude of this set velocity. Since the desired swimming speed towards the target, relative to the initial location, is negative (along the negative x-direction) in the coordinate system in figure 1, a negative sign is added to u set . The value of T set is obtained from numerical tests to yield stable control of swimming speed. Based on the proposed scheme in equation (18), the plate would swim at a high velocity when the distance G(t) is large and slow down when approaching the target position. In case of an overshoot in position, reversing the sign of u set would lead to reducing the tracking error, and subsequent deceleration due to the reduced actuation force. And the plate would even be pushed back by the incoming flow if the produced thrust is overwhelmed by the drag. Nevertheless, the gap value will become positive after being pushed back, driving the plate to move forward again. When the plate is approaching the target and the value of u set is very small, the controller instructs the plate to move at a nearly zero mean speed during each actuation period. We employ an incremental calculation of the adjustment factor α (t), and the current value of α (t n ) is obtained by where α (t n−1 ) denotes the value of the last time step and Δα (t n ) is the error increment which is calculated by the PD control action where t s is the sampling period. The value of α calculated by equation (21) is limited to [0, 1] so that the amplitude of the actuation force stays between 0 and the specified maximum bound, while the force direction will not be changed. The values of k p and k d (listed in table 1) are obtained from numerical studies based on the Ziegler-Nichols method (Ziegler and Nichols 1942) to yield an agile response and a small overshoot.

Computational domain and self-consistency study
The computational domain is presented in figure 4(a). On the flexible plate surface, the non-slip/no-flux condition is imposed, while for the other boundaries, the non-reflective far-field boundary condition is applied. A self-consistency study is conducted to assess the sensitivity to the mesh density and time step size by solving the laminar flow around a tethered flexible plate at Re = 500, m * = 0.5, K = 0.5, f * = 2.5, and d 0 = 1.0L. Three meshes are generated: a coarse mesh with 82 574 cells, a mediumsize mesh with 103 888 cells, and a fine mesh with 130 533 cells. The first grid height of the three meshes is 0.005L. The structural mesh contains 132 fifteennode wedge elements. First, the mesh convergence test is performed in which the three generated meshes are used along with the non-dimensional time step size, defined as Δt = ΔtU/L, Δt = 0.003 33. The results of C T yielded by the three meshes are presented in figure 5(a). It can be seen that the medium-size mesh yields result very close to that of the fine mesh. Thereafter, the effect of the time step size is also studied when the medium-size mesh is used, as shown in figure 5(b). These results suggest that the simulation results are not sensitive to the mesh size and time step size if they are sufficiently small. Therefore, the  medium mesh is used for the following simulations to ensure computational accuracy, and meanwhile, reduce the computational cost.

Results
Considering that a fish tends to hold station in front of the cylinder in a uniform flow, as reported in an experimental study by Liao et al (2003), we focus on the swimming of a flexible plate whose actuation force amplitude is dynamically adjusted in front of a cylinder via a feedback controller, as shown in figure 1. The swimming and station-holding in the absence of a cylinder in a uniform flow condition are also examined for comparison. Before we start the FSI simulations involving motion control of the plate, the fluid field around a rigid and static plate near the cylinder is first studied. Then the effect of Reynolds number (Re = 500, 1000, and 2000) and the initial distance between the plate and cylinder, d 0 , is considered with the fluid-only simulation results as the initial fluid field. Next, we explore how the actuation frequency f * and the phase shift ϕ of the actuation force would influence the swimming performance and the control action. Some parameters remain unchanged for the following simulations, which are listed in table 1. It is noted that the maximum magnitude of the actuation force C ef0 is obtained through numerical tests to bend the plate to produce a large enough thrust to swim forward, but meanwhile avoid a too large magnitude that yields exaggerated deformation which may cause numerical issues.

Flow over the rigid plate at rest in front of a cylinder
Before we initiate the actuation force to deform the flexible plate, the flow over the cylinder and plate at rest is firstly simulated. Taking Re = 1000 as an example, the mean drag coefficiencies C d of the cylinder and the static plate at different d 0 is shown in figure 6. It can be seen that the presence of a cylinder in the wake results in a reduced drag force of the static rigid plate compared with the cases when it is placed solely in a uniform flow. The closer the plate is to the cylinder, the more remarkable this drag decrease is. The cylinder also gains a similar drag reduction when it is placed behind the plate. Nevertheless, this benefit does not present a monotonic variation with respect to the distance between the two bodies, and a maximum drag decrease is seen at d 0 = 2L. After this minimum value, C d of the cylinder grows again as the plate is placed farther away from it. This may be related to the unsteady vortex shedding of the cylinder which is influenced by the upstream plate in a complex way at different gaps. A similar variation pattern of drag force of the cylinder ahead of a plate was reported in Wu et al (2014).
To get an insight into the surrounding fluid field, we plot the pressure distribution around the plate and cylinder in figure 7. It can be seen that the trailing part of the plate is exposed to the high pressure attributed to the cylinder when it is placed near the cylinder, as shown in figure 7(a), which helps to reduce the drag force. In comparison, when the plate is located further upstream of the cylinder, it loses the relatively high pressure which pushes it along the negative x-direction, as shown in figure 7(b). Meanwhile, the flow-incident surface of the cylinder also withstands larger pressure, which is similar to the scenarios when it is solely in a uniform flow as shown in figure 7(c). This high pressure on the left-side surface of the cylinder leads to a larger overall resistance compared to that when it is near the wake of the plate. It is interesting to further investigate how this hydrodynamic interaction between the plate and cylinder would influence the flexible deformation and free-swimming of the plate when the actuation force is applied under the feedback control.

The effect of Re and initial distance d 0
After the above simulations of the static plate without flexible deformation, we now explore the selfpropulsion performance with the control objectives, i.e. swimming to a target position and stationholding, being implemented with k p = 7 and k d = 0.05. The results in section 4.1 when the plate is rigid and static are taken as the initial fluid field for the following simulations involving flexible flapping and self swimming. The following results are shown in figure 8: the overall energy expenditure coefficient C Ps and time T s required to travel the distance of one body length L from different initial positions, the mean energy expenditure coefficient C Ph , and the converged force factor α h when the plate holds the station near the. The results for an isolated plate in the absence of the cylinder in the flow are also included for comparison. It can be found that under the present control scheme, generally, the less overall energy is consumed by the flexible plate to travel to the target position when it has an initial location closer to the cylinder behind. When the initial position is closer to the cylinder, it is also more energy-saving for the plate to stay there, reflected by a smaller C Ph at a smaller d 0 . The curves of swimming time T s to reach the target position present more complex variation patterns as d 0 and Re vary.
A closer inspection on figure 8(a) reveals that the variations of C Ps are more sensitive to the initial distance d 0 before d 0 exceeds 1.0L, after which the curves of C Ps are flat as d 0 increases. This trend is more noticeable for Re = 500, where C Ps remains almost unchanged after d 0 /L is larger than 2. The greatest energy saving compared to an isolated plate that travels the same distance is seen at d 0 = 0.3L and Re = 1000, with a reduction of 39% in C Ps . In  addition to the initial position, the Reynolds number is also a significant factor in determining the swimming energy consumption. For example, 183% more energy is needed to swim the same distance at d 0 /L = 0.75 when Re decreases from 2000 to 500. Generally, it is more energy efficient to swim in the  The variation patterns of C Ps generally apply to the mean energy cost coefficient C Ph during the stationholding after reaching the target. The closer to the cylinder the target is, the less energy is required by the plate to stay there. A larger Re also contributes to the reduction of energy consumption. Not surprisingly, the curves of the converged force adjustment factor α h show similar patterns as d 0 and Re are varied to that of C Ps and C Ph whose calculations are dependent on instantaneous force amplitude.
The instantaneous swimming distance relative to the initial position, the swimming velocity, and the force amplitude factor α produced by the controller, at d 0 = 0.3L and d 0 = 6.0L, are presented in figures 9-11, respectively. It can be found that the swimming distance and velocity curves when d 0 = 6.0L are close to that of an isolated plate, implying the hydrodynamic benefit from the cylinder upstream is rather weak at this remote initial position. The plate is firstly drifted downstream at the beginning of the swimming (before t = 11T) when the plate departs from a position far from the cylinder. This is related to the small force amplitude (α) in the early phase as presented in figure 11, under which the generated thrust is not large enough to overcome the drag. But this downstream drift is not seen at d 0 = 0.3L and is significantly reduced at a higher Reynold number, 2000. Under the closed-loop feedback control, the plate experiences remarkable   acceleration and reaches the maximum swimming speed at around t = 17T when the maximum actuation force is applied as shown in figures 10 and 11. After this peak swimming velocity, the exerted actuation force is reduced under the feedback control as the plate approaches the target. As a result, the plate decelerates and swims slowly towards the target to avoid significant overshoot.
The swimming speed during the holding station at the target is relatively small and fluctuates upstream and downstream symmetrically, as plotted in figure 10. We find that the converged swimming velocity at the larger Re (2000) during the position-holding phase seems more stable in contrast to a lower Re (500) in terms of smaller undulation amplitude (±0.01U versus ±0.02U). Furthermore, a smaller d 0 also yields a smaller undulation amplitude of the swimming speed. In general, the implemented feedback controller is shown to be effective in the motion control of the considered flexible plate involving dynamic FSI, evidenced by the well-converged plate position, swimming speed, and actuation force amplitude factor, as shown in figures 9-11.
In addition to the above swimming performance, we also study the flexible structural deformation due to the dynamic interplay of the structure and the fluid with the actuation force being applied. The tip displacement of the trailing edge of the plate is depicted in figure 12. The maximum tip displacement of the plate at Re = 500 is greater than that at Re = 2000 attributed to a larger actuation force applied on the surface during the station holding. A closer initial location, and thus a closer target position relative to the cylinder, seems to have little effect on the maximum amplitude of the displacement. Nevertheless, a more noticeable oscillation of the displacement amplitude is presented at d 0 /L = 0.3 and 1.0, but not at d 0 /L = 6.0 for both Reynolds numbers, which is marked in figure 12. It should be mentioned that the normalized displacement oscillation period T oscillation U/L, defined as the time interval of the maximum trailing-edge tip displacements of the plate as shown in figure 12, at d 0 /L = 0.3 is 4.79 and 4.27 for Re = 500 and 2000, respectively, which are close to the vortex shedding period from the cylinder, i.e. 4.69 and 4.39, respectively. Therefore, we conjuncture that this displacement oscillation is mainly caused by the nearby downstream periodic vortex shedding from the cylinder surface. This demonstrates the interaction between the flexible deformation and nearby vortices, which is not presented if the plate is far from the cylinder.
For further insight into the flexible deformation of the plate under actuation force, we plot the envelope of the midline in one actuation period in figure 13. The periodic body deformation is generally symmetric about the horizontal balanced position (y-axis), which shows similarity to the body deformation patterns of a real larval zebrafish (Müller and van Leeuwen 2004) under the present force actuation mechanism in terms of the arch-shaped midline at the posterior body (see figure 2 in their work). The observed tip deformation at Re = 500 is larger than that at Re = 2000 due to a larger actuation force as shown in figure 11. At t = 0.5T and t = T when the instantaneous force amplitude is zero, a more noticeable arch near the middle of the body length is seen at Re = 500. In comparison, the body at Re = 2000 presents a more flat profile at this instant.
The instantaneous Z-vorticity contour around the plate during station-holding at two different d 0 is depicted in figure 14. A pair of vortices, each with clockwise and anti-clockwise direction, is shed from the trailing edge of the plate, roughly at t = 0.25T and t =0.75T, when the plate flaps reversely from one extreme side to another side, respectively, during one period. The effect of the presence of the fixed cylinder is mainly reflected by the morphologies of the wake vortices. Namely, the vortices behind the plate are more compressed towards the negative x-direction at d 0 /L = 0.3, which presents a 'V' style layout from the plate tail to the downstream till the cylinder. In contrast, at d 0 /L = 6.0, when the plate is far from the cylinder, the vortices show more regular patterns which have been seen behind a flapping foil or plate in static fluid and uniform flow , Zhu et al 2014. This difference is caused by the high pressure on the left-side surface of the cylinder, as shown in figure 15. Due to this high pressure, the low-pressure zone near the tail of the plate is reduced in size at d 0 /L = 0.3 compared with that at d 0 /L = 6.0, which can be found in a comparison of figures 15(a) and (b). Apart from this noticeable reduction of low-pressure, the other pressure distributions between these two are generally identical. This would help the plate gain a drag reduction, and thus, less energy is required to hold the station for d 0 /L = 0.3 as presented in figures 8(c) and (d). The vortices patterns and the pressure distribution at Re = 2000 are qualitatively similar to the counterparts at Re = 500. Therefore, they are not shown here.

The effect of force frequency f * and phase shift angle ϕ
In this section, we study the effect of f * and ϕ on the swimming performance and station-holding capability with Re = 500 and k p = 15. We now use a larger k p than k p = 7 in section 4.2 to yield a faster initial acceleration of the plate. This aims at avoiding the Figure 16. The overall energy expenditure coefficient C Ps (a) and time T s (b) required by the plate to reach the target by traveling the same distance L, the mean energy expenditure coefficient C Ph (c), and the converged adjustment factor α h of actuation force amplitude (d) when the plate holds station near the target at different initial distances relative to the cylinder d 0 at different f * with Re = 500 and ϕ = 0. numerical divergence due to the clashes of the cylinder grid cluster and the plate grid cluster for some cases considered in this section. This clash is caused by the downstream drift of the plate by the incoming flow when the produced thrust is overwhelmed by the drag force at the initial swimming phase if a small k p is used.
The results of C Ps , T s , C Ph and α h at different force frequencies are summarized in figure 16. As can be seen, the actuation frequency has a significant effect on all four parameters. Specifically, much more energy and swimming time are required at f * = 2.0 when d 0 /L is beyond 1 for the plate to travel the same distance and hold the station at the target compared with other f * . The most energy savings and time savings are yielded at f * = 2.5 for d 0 /L > 1. A larger force amplitude is needed to swim forward and hold the station at f * = 2.0, which even approaches the maximum allowed force amplitude (α ≈ 1) at a large d 0 , e.g. 6.0L. Our simulation also shows that an isolated plate in the absence of the cylinder is unable to swim forward or hold the position, and it is drifted downstream, at f * = 2.0 under the current control scheme even with the maximum force amplitude C ef0 = 2.0 (α = 1).
A close inspection on figure 17(a) reveals that the reduced C Ps at f * = 1.5, compared to f * = 2.0, may be attributed to a faster arrival at the target position, despite that a larger maximum swimming speed is seen at t = 10T for f * = 2.0. By comparing figures 17(c) and (d), we also find that a larger force amplitude and higher actuation frequency at f * = 2.0 do not lead to greater tip amplitude than that at f * = 1.5, as demonstrated in Dai et al (2012), implying a complex interaction between the flexible deformation and actuation parameters. Figure 18 depicts the instantaneous Z-vorticity contours around the plate and cylinder during station holding near the target at d 0 /L = 0.3 under the two different actuation force frequencies. The vortices distribution around the left-side surface of the cylinder is more disturbed by the shedding vortices from the plate tail at f * = 1.5. Thus, the shedding vortices by the cylinder have a less regular configuration compared with that behind the cylinder at f * = 2.0 and 2.5, as shown in figures 14 and 18. This may be attributed to the greater oscillation amplitude of the plate's tip at f * = 1.5. The maximum tip displacements for the three frequencies are 0.147L, 0.134L, and 0.127L at f * = 1.5, 2.0 and 2.5, respectively. Larger flapping motion at the trailing edge enhances the capability of the plate to communicate momentum to the surrounding flow (Michelin and Llewellyn Smith 2009), thus leading to a stronger influence on the vortices shed from the cylinder at f * = 1.5.      We next explore the effect of the phase shift ϕ of the actuation force in equation (2) at Re = 500, f * = 2.5, and d 0 /L = 1.0. The results of C Ps , T s , C Ph and α h are presented in figure 19. The phase shift has a similar effect on C Ps and T s as ϕ increases. Under the parameters studied, a maximum and minimum value of C Ps and T s are seen at ϕ = 60 and 120 deg. In comparison, the energy-saving and swimming timesaving are roughly 16% from ϕ = 60 deg to 120 deg. Nevertheless, the phase shift ϕ has a negligible effect on C Ph and α h in which the maximum difference caused by ϕ is below 1%. Therefore, the phase shift generally only affects the swimming phase from the initial location to the target, while the station holding is far less sensitive to ϕ. Likewise, the influence of the phase shift on the instantaneous variations of swimming distance x s , swimming speed u b , and force adjustment factor α is mainly noticeable during the swimming phase, as shown in figure 20.

Conclusions
The swimming and station-holding performance of a flexible plate in front of a fixed cylinder in a uniform flow is numerically studied by using a high-fidelity FSI solver. The deformation of the two-dimensional plate is actuated by the imposed actuation force that mimics the muscle force of the fish body and fin. The amplitude of this actuation force is dynamically tuned by a feedback PD controller to achieve the control objective, i.e. swimming from the initial position to a designated target point and holding station there. The simulation results show that the proposed control scheme to adjust the actuation force amplitude, based on swimming speed and the instantaneous gap between the plate location and target, is effective and robust when coupled with the FSI solver. The plate successfully achieves the control objectives of reaching the target and holding the station stably despite the drag force due to the incoming flow at different flow environments (Reynold numbers, and initial position relative to the cylinder) and different actuation force formations (force frequencies and phase shift angles).
Nevertheless, these flow conditions and force actuation formulas indeed result in performance changes, in terms of swimming time T s and energy consumption C Ps traveling over the same distance (one body length) from the initial position to the target, and the energy expenditure C Ph required to hold the station at the target. We find that energy consumption is reduced for swimming and station holding at a higher Re compared with a smaller one. The plate which departs from an initial position d 0 closer to the cylinder can also benefit energy savings reflected by smaller C Ps and C Ph , compared to that with larger d 0 at the same Reynolds number, force actuation formula, and control scheme. This may be attributed to the reduction of the low-pressure zone in size near the trailing edge of the plate (see figure 15) due to the downstream cylinder, which provides additional thrust boost (or drag reduction) to help save energy. Meanwhile, the flexible deformation of the plate is also affected by the nearby cylinder flow field, reflected by the frequency resonance of the oscillation of the maximum tip displacement with the shedding vortices from the cylinder (see figure 12). These benefits and influences are reduced when the plate is far away from the cylinder. If d 0 is large enough, the swimming time and energy consumption are almost the same as that in a uniform flow in the absence of the cylinder.
The force frequency has a significant effect on swimming performance. Under the parameters studied, much more swimming time and energy consumption are required at f * = 2.0 compared with other frequencies when all the other variables are the same. Even at f * = 2.0, an isolated plate in a uniform flow fails to produce sufficient thrust to overcome the drag force to swim forward, despite the applied largest force amplitude allowed in this work. In comparison, the effect of the phase shift ϕ of the actuation force is less significant. To be more specific, the difference in swimming energy consumption, C Ps , from the initial to the target position attributed to varied ϕ is below 16%, while it has a negligible effect on station-holding performance after the plate reaches the target.
In this study, a simple controller is used for motion control, which proves to be effective for the considered control objectives under the simulation environment. Refined PD tuning or even more advanced control strategies are anticipated for more complex performance specifications. It is also noted that the plate is restricted to move only in the axial direction in this work. Our numerical tests show that there may be a non-negligible lateral movement if it is considered due to the phase shift between the applied force and the actual bending of the plate under the present actuation scheme. Specific measures may be needed when the lateral movement is considered to reduce the deviation, e.g. by introducing an additional controller to tune the force direction, which would be an interesting topic for future study.
(www.archie-west.ac.uk) based at the University of Strathclyde. The first author thanks China Scholarship Council (CSC) for financial support during his study in the UK.

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

Appendix. Validation
This validation involving a self-propelled flexible plate in the wake of two tandem cylinders described in Wang et al (2019), as shown in figure 21. In this simulation, a prescribed motion at the leading edge of the plate is imposed as H (t) = H 0 sin(2πft) with H 0 being the maximum heave amplitude and f denoting the heave frequency, while the trailing edge of the plate is set free. The distance between the two cylinders D x is 4.5D with D being the diameter of the two cylinders. According to the studies in Wang et al (2019), if the plate is able to hold stationary in the wake, the flapping frequency needs to be identical to that of the vortex shedding. Therefore, we firstly perform a fluid-only simulation, i.e. the flow over the two tandem cylinders without the flexible plate in the wake of them, and obtain the vortex shedding Strouhal number, St = f v L/U 0 = 0.176 (with f v being the vortex shedding frequency, L being the length of the plate and U 0 denoting the incoming flow velocity), which is close to the value 0.185 reported in Wang et al (2019). Then the dimensionless heave motion frequency is defined as f * = fL/U 0 = 0.176. In this validation, the governing parameters matching with the counterparts in Wang et al (2019) are the Reynolds number Re = U 0 L/ν = 200, the mass ratio σ = ρ s h/ρ f L = 1, and the bending rigidity K = EI/ρ f U 2 D 3 = 0.5.
There are three typical modes of the plate in the wake of the cylinders, i.e. drift upstream (DU), DD, and holding stationary (HD) depending on the flapping amplitudes H 0 and initial horizontal positions of the plate G 0 , which are summarized by Wang et al (2019) as shown in figure 22. For the HS mode, we quantitatively compare the simulated results of the horizontal positions of the leading edge of the plate when it is holding station at H 0 = 0.3L with that reported in their work in table 2. The typical trajectories of the three motion modes in our simulations are shown in figure 23. Therefore, it can be seen that our simulation results present good agreement with the counterparts in Wang et al (2019).