Numerical Simulation of a Multi-Body System Mimicking Coupled Active and Passive Movements of Fish Swimming

: A multi-body system model is proposed for the mimicking of swimming fish with coupled active and passive movements. The relevant algorithms of the kinematics and dynamics of the multi-body system and coupled fluid solver are developed and fully validated. A simplified three-body model is applied for the investigation of the hydrodynamic performance of both an active pitch motion and passive movement. In general, there is an optimal stiffness, under which the model swims with the fastest velocity. The effect of the damper can be drawn only when the stiffness is small. Comparing with the rigid tail, the flexible tail leads to a faster speed when the stiffness and damping coefficients are in a suitable range.


Introduction
The propulsion and maneuvering abilities of fish [1,2] have attracted a lot of attention in the past, as they have various applications to industry related products, such as robotic fish [3][4][5] and autonomous underwater vehicles (AUV) [6][7][8]. Multidisciplinary knowledge is required to fully explore the inherent mechanism which involves fluid mechanics for the estimation unsteady hydrodynamics force, structure analysis for the evaluation of fish body deformation and biology concepts for the determination of various signal transformations via the neurological system. Consequently, a completed and comprehensive replication of a fish's movement system must consider all the above-mentioned aspects, and particularly, couple the fish's internal muscle forces with its external hydrodynamics force and the flexible fish body deformation controlled by the nervous system [9,10].
Due to the complexity in the mechanism of fish swimming, existing investigations on this subject are limited and could be roughly divided into following two steps, in the order of increasing physical complexity and computational or experimental novelty gradually: (i) a prescribed locomotion of swimming fish [11][12][13] and (ii) a passive motion with the lower order degrees of freedom (DOF) [14][15][16]. The latter one generally considers interactions between external hydrodynamic forces with fish body movements. However, as a real fish in nature, both active locomotion and passive deformation are fully coupled. Therefore, a comprehensive model is required which includes both active and passive movements of fish swimming as well as their interactions with external flow.
On the other hand, from the point view of robot design, a whole body can be divided into several segment as a multi-body system with hinges (links) connecting those segments either actively or passively. Recently, a series of relevant work has been reported in the literature. Eldredge and co-workers [17,18] developed an articulated three-body model with either prescribed angular movements of the two hinges or passive hinges to mimic the undulatory motion of a slender body, and its interaction to incoming flow was investigated by the viscous vortex particle method. Kajtar and Monaghan [19] studied the motion of three linked ellipses moving through a viscous fluid in two dimensions coupled with their smoothed particle hydrodynamics solver and the angles between the ellipses change with time in a specified manner. An advanced algorithm for a multi-body system coupled with slender body theory for estimation of hydrodynamic forces was developed by Boyer and co-workers [20][21][22], and various control methods, e.g., active, passive, etc. and various morphologies of multi-body, e.g., serial for fish-like robot, tree for bird/insectlike robot, etc. were considered in their model. Though practically fish motion is propelled by muscle stimulation instead of an electric motor, it is complicated to produce a manmade muscle to fully mimic the real animal muscle. In this situation, the multi-body system with participations by internal torque, body stiffness, and phase delay, etc., is of importance for providing a clue on advanced biomimetic investigation. The undulation motion can be mimicked accurately by adjusting the ways of control on linked hinges, and the more segments and hinges applied, the more precise undulation motion is achieved. Moreover, the coupling procedure with fluid solver is straightforward, and hydrodynamic performance can be obtained by the mature CFD techniques.
Inspired by the work mentioned above, a multi-body model mimicking fish swimming is proposed in current work, in which coupled active and passive movements are considered. In addition, a Navier-Stokes equations solver is also applied to fully resolve the hydrodynamic force acting on the multi-body model. On one hand, as both active and passive movements are considered, the obtained results may provide a comprehensive understanding of the inherent mechanism of fish swimming. On the other hand, as a prototype of robotic fish the obtained results may also provide sufficient information for robotic fish design and relevant control strategies as well. In the following, the proposed multi-body model is first introduced in Section 2 and the simulation methods are also briefed. In Section 3, two typical validation cases are given and the results are compared with those in the literatures. Results and discussion of the kinematics and hydrodynamic performance of current multi-body model are presented in Section 4, and the conclusions are drawn in Section 5.

Coupled Active and Passive Model
In general, the real fish body and its movement can be well represented by the multibody system with a serial of rigid bodied connected by active and passive hinge joints. The preliminary model is shown in Figure 1a, in which the active hinges can be used for the prescription of active movement of fish body and the passive hinges can be used for the representation of internal flexibility (elasticity) of fish body and tail.
For simplicity, a three-body model as shown in Figure 1b is employed in current study. Each rigid body has an elliptic shape with aspect ratio of 0.1 and chord length C. The left two rigid bodied are connected by a passive hinge and the angle between centerlines of the two bodies are r1, whereas the right two bodies are connected by an active hinge and the evolution of the angle r2 is prescribed.

Fluid Solver
The governing equations for fluid flow around the multi-body model are the twodimensional (2D) continuity and momentum equations for incompressible viscous fluid, which are: where u is the velocity vector of flow, p the pressure.  and μ are the density and viscosity of the fluid, respectively. The governing equations are solved by commercial software ANSYS Fluent. The features of dynamic mesh and user defined function are activated for calculating and capturing body motions. The forces and moments that exerting on moving surfaces by surrounding fluid are expressed in terms of the integrated pressure and viscous stress, and both of them are taken into account at every time step in this study. The size of computational domain is 42C × 23C. The boundary condition on the rigid segment surface is set as no slip wall boundary, an inlet velocity with magnitude of zero is defined on the left side of computational domain, and a pressure outlet is applied on the right boundary as shown in Figure 2. Regarding to the dynamic mesh method, the remeshing and smoothing parameters are both chosen carefully. The re-meshing is accomplished by the local cell method, and the smoothing process is done with a diffusion function. The parameters applied are all well tested by a mesh density independent test. A 2D pressure-based transient fluid solver is selected, and the fluid field is set as laminar viscous model. A fractional-step method (FSM) scheme is activated under the pressure-velocity coupling panel. The spatial discretization of both the pressure and momentum are with the second order upwind accuracy.

Dynamics of Multi-Body System
In order to capture and derive the kinematics of the multi-body system, a basic vector, Xstate is employed to characterize essential motion variables of the system in global coordinate, and it includes information of the position of reference body, linear and angular velocities of reference body, and angular velocity of all the hinge joints. Meanwhile the information of input joint conditions is represented by Pcontrol, and it includes the pitching acceleration for active joint and the torque applied on the passive joint. By Newton's second law and also the information of external forces fext (hydrodynamic forces from fluid solver), the kinematics of the multi-body system can be updated consequently. For simplicity, the kinematic evolution of the multi-body system can be summarized by a function as: (2) where Noutput indicates the variables to be obtained including information of torque of actuating the active joint, and angular acceleration of passive joint.
However, concerning the kinematics of each segment in the system, the forces on the hinges connecting the rigid body are still unknown. Following the algorithm proposed by Boyer and Porez [20], three recursive loops are conducted which includes the first forward recursion on the kinematics, the backward recursion on the external forward dynamics and the second forward recursion on the internal dynamics. The detailed illustration of the implementation of the algorithm can be found from Reference [20][21][22].
Moreover, the coupling process between multi-body system algorithm and fluid solver is made through an interactive data transferring between two solvers with fluid forces and moments (fext) of fluid solver and statement vector (Xstate) of multi-body system. Basically, the simulation iteration loop starts from updating the imposed position which is available from the last time step, and fluid domain is solved for obtaining hydrodynamic forces and moments, and then with the fluid forces and moments as input conditions, the multi-body dynamic solver calculate updated position for next time step. Especially, before the next time step starts, a fourth-order predictor and fifth-order corrector time-explicit discretization method of Equation (3) is utilized to achieve an accurate solution. (3) where t is the time step and the symbol ( ) represents time derivative.

Validations
To assess the reliability of our coupling method with both the multi-body system algorithm and fluid solver, two validation cases, with either active joint or passive joint models, are carried out. To sum up, the results agree well with the previous numerical work.
The multi-body system algorithm with active hinge joints is first validated with the case presented by Eldredge [17], in which a three-body rigid system (similar to the model shown in Figure 1b) with active hinges controlled by prescribed angles in quiescent fluid is simulated. In particular, the model is built as a massless articulated system and each rigid body has an identical elliptical section area of aspect ratio 0.1 with chord length c, and the distance from tip to hinge is set to 0.1c. Each hinge between the pair of the rigid bodies is independently controlled, with the prescribed angle as a function of time (t), The entire linked system has 3DoF in x, y and pitching directions, i.e., it is free to move in all directions in two dimensions under propulsion of linked segments controlled by prescribed angles. The undulation Reynolds number is fixed at 200. The results are well matched with Eldredge [17] in translational velocity components as shown in Figure 3. To validate the multi-body system algorithm with passive joints, the model with two rigid segments linked by a torsional spring is simulated and validated with the results by Toomey and Eldredge [18]. The upper rigid body has been imposed with sway and pitch motion, and the motion of the lower body is governed by fluid forces and constraint of the linked hinge. Our results of induced angle and fluid forces agree well with numerical simulations as shown in Figure 4.

Results and Discussions
In this section, the multi-body system is simulated with both active and passive joints, aiming to investigate tail flexibility effect on propulsion performance, where the tail flexibility is realized by applying different spring stiffness and damping coefficients at passive joint, and the analysis is made from the hydrodynamic points of view with considerations of fluid forces exerting on the model components, power consumption as well as vortex structures.
We assume that the left two rigid bodies mimic the fish tail with flexibility consisting of a passive joint in the middle, and the right rigid body is fish body with an active joint as the tail peduncle. There is a prescribed pitch angular motion between fish body and tail at the peduncle joint, mimicking fish tail that flaps to propel the entire body. This prescribed relative motion between head and tail is same as r2 in Section 3. The passive joint is a stiffness-damper spring, its torque (τ) is governed by following equation: where R* and K* are damping and stiffness, respectively. The damping and stiffness can be represented by the non-dimensional coefficients as: R  R * / ( fc 4 ), (6) where f is the flapping frequency. Referring to the stiffness and damping coefficients selected by Toomey and Eldredge [18], K and R vary between 5.1 to 51.4 and 0.2 to 0.7 respectively. In the current study, the parametric test is carried out with combinations of six stiffness coefficients (varying from 1 to 27) and two different damping coefficients (0 and 0.245). The smaller stiffness leads a softer tail. One more case with rigid tail is also included as a baseline for comparison. The rigid tail is designed by prescribing a zero-pitch angle on the passive joint all over the cycles.

Kinematics of the Multi-Body System
The trajectory of entire system and induced instantaneous pitch angle at the passive joint are the most desired kinematic quantities which can be directly measured from the current simulations. As shown in Figure 5a, the trajectory is tracked by monitoring the location of the hinge joint between articulated system body and tail, i.e., the position of the active joint. The system motion starts in quiescent fluid condition, and accelerates gradually into a quasi-steady stage under undulatory propulsive motion. The induced motion of the system follows a zig-zag trace, and the moving direction is determined mainly by phase shift of the induced pitch at passive joint under different stiffness and damping coefficients. The orientation angle (α) of the motion trace line shown in Figure  5b, presenting the moving direction within global coordinate, is quantified by inverse tangent formula (α = atan(Y/X)) using the trajectory in Figure 5a. As plotted in Figure 5, the induced motion reaches the quasi-steady status after approximately 10 revolutions. The system undergoes a development procedure of balancing friction force and thrust force. At start-up stage, thrust force is bigger than friction force, and causes accelerated motion. Theoretically, in the Stokes regime, viscous force is proportional to moving velocity, so the viscous force and thrust force can balance with each other when the velocity increases, then leads to a stable quasi-steady stage. Although the studies on the start-up stage are of importance to understand the mechanisms of maneuverability and stability, the resulting analysis in the following sub-sections will mainly focus on the fully developed stage. The average values of the variables yielding global analysis can provide a general view of the induced locomotion under specific undulatory body posture. It is noted that velocities and forces in the following sub-sections are all transformed in semi-local coordinates with the axis along the travelling direction and perpendicular direction, which are decided by the average induced angle (average α in Figure  5b) in the quasi-steady stage.
The velocity components, Vx and Vy, are presented in the traveling direction and the perpendicular direction respectively. The instantaneous velocities in the quasi-steady stage with selected combinations of stiffness and damping coefficients are plotted in Figure 6. It is shown that all the amplitudes, mean values and phases of Vx vary with the parameters, while the amplitudes of Vy are slightly different with mean value remaining zero. In addition, the phase difference of Vx is a consequence of the phase shift of the induced pitch motion at the passive joint. The horizontal velocity component Vx is one of the most important variables that used to measure the propulsive performance. Therefore, the average Vx is calculated and plotted as in Figure 7. It can be seen that the mean travelling velocities increase dramatically when the stiffness coefficient is below 3.95, and decrease gradually after a mild rise to the peak point. The effect of the damping coefficient on the induced velocity is dependent on stiffness coefficient. The average velocities are larger with the damper applied when the stiffness coefficient is smaller than 11.05, while the differences disappear when the stiffness is bigger. It is interesting to observe that the articulated multi-body models in most cases travel faster than the one with rigid tail, and there is an optimal stiffness coefficient, leading to the fastest swimming velocity, as a result of appropriately induced angle and phase at passive joint. The increase of velocity with flexible tail agrees with findings in the previous work from Bergmann et al. [23], in which similar findings have been obtained through examining caudal fin elasticity effect by changing lumped parameters. The instantaneous induced angles at the passive joint (r1) under the selected parameters are shown in Figure 8. The angles are periodic in each revolution, so the amplitude (ramp) and phase (φ) are arranged in Figure 9 exploring the trends under different parameters. It can be seen that the stiffness plays an important role on both ramp and φ. It is reasonable that ramp is bigger when the joint is less stiff, though a spike occurs when the joint stiffness is further reduced. The magnitudes of all ramp are smaller than the imposed angle amplitude (57°) at the active joint. The phase (φ) decreases with bigger stiffness, and it shows that the induced motion turns to be more consistent with the active motion when the tail is stiffer. With the damper applied, the induced angle has smaller amplitude, and the impact of damping coefficient on ramp is dependent on the stiffness, that the difference of ramp with or without damper turns to be smaller when the joint stiffness becomes bigger. The damper causes a delayed action of the passive joint, and that leads to a bigger phase change.  The induced angle at passive joint is always a result coupling with external surrounding fluid and internal elements. It is noticed that the flexibility of the spring determines whether the induced pitch is dominated by internal or external variables at the passive joint. One example for internal properties domination is that when the spring stiffness is big enough, especially in the case of the rigid tail, there is no relative pitch motion induced between two connected rigid components. The external environment takes in charge of the induced motion when the spring stiffness is small, such that in a case with the passive joint as a fully revolute joint with no stiffness and damping, the soft part of the tail would just follow the fluid pattern generated by wake structure and bypass flow. This can explain the observation in Figure 9a that the induced angle turns smaller before and after the stiffness coefficient (K) of 1.97, as it transfers from external dominating condition to internal dominating conditions.

Hydrodynamic Performance
The instantaneous torque history is plotted in Figure 10, where τ1 is restoring torque induced by the interactions of the elements and fluid, and τ2 is input torque produced by electric actuator motor at active joint. τ1 is obtained from the stiffness-damper spring equation (Equation (4)), and τ2 is an output from the function of Equation (2). The torque curves are periodic and there are phase shifts between different cases, the average torque is approximately zero, and the amplitude mainly depends on the stiffness. It can be observed from Figure 11 that the torques of active joints are larger than those of passive joints, and all of them increase dramatically before reaching a stable level with increased spring stiffness. The torques of both active and passive joints are slightly bigger in the cases with flexible tails than those with rigid tail.  The power can be used to quantify energy consumption for propulsion. In the current multi-body system, the power applied at active joint is the only input energy resource, which contributes to the kinetic motion and can be consumed by the induced motion at passive joint only when damper is considered. The power, P, is calculated as the following: The results of the power history and average power are plotted in Figures 12 and 13, respectively, where P1 is power at passive joint and P2 is input power at active joint. It is seen from Figure 12 that the power curves at passive joint under different parameter selections have similar amplitudes, but phases of the curves shift. The power at the active joint is constantly positive whereas the amplitudes change. The positive and negative values denote the spring can either store power (negative P) or release power (positive P). The phase shift of power curve is a consequence of that of induced pitch angle. The average power at the passive joint without damper is approximately zero as shown in Figure  13a, indicating the spring only store or release power without power consumption. While with the damper, the average power turns negative, showing that power loss exists.  The flow structures (vorticity distribution/contour) around the multi-body system are illustrated in Figure 14. It is clearly shown that the typical reverse von Kármán vortex street is observed, although the undesired vortex leaks from the gap between adjacent elements. The strength of the vortex street can enhance the thrust force, and the wake structure is highly dependent on the undulatory profile, especially the flexibility of the tail which can be interpreted by the selections of input parameters. To sum up, the cases with faster swimming speeds show stronger vortex strength.

Conclusions
A multi-body system model is proposed for the mimicking of swimming fish with coupled active and passive movements. The relevant algorithms of the kinematics and dynamics of the multi-body system and coupled fluid solver are developed and fully validated. The test cases in current study are simplified which is a three-linked rigid segment system with an arbitrary sinusoidal pitch motion applied at active joint, and a linear spring model is applied at passive joint to mimic flexible tail. The flexibility mechanism can be examined from this schematic.
In general, there is an optimal stiffness, under which the model swims with the fastest velocity. The effect of the damper can be drawn only when the stiffness is small. The damper can shift the phase of the induced angle by delaying the response of the spring, and hence change the propulsion posture, which causes different swimming speeds. Comparing with the rigid tail, the flexible tail led to a faster speed when the stiffness and damping coefficients are in a suitable range. The properties of the spring stiffness and damper, the induced/prescribed pitch angle and torques are typical variables from the aspect of internal dynamics. The fast speed is a result of the increase of lateral forces produced by larger input torque, and also a consequence of less power loss due to a properly induced undulatory swimming pattern, which are all related to the subject of external dynamics. It can be seen from the test case that the variables from both the internal and external dynamics can be clearly illustrated.
The analysis of the internal variables effect on the external behavior should be more enhanced with further investigations. Future work would be to focus on activating advanced ways of control, and on implementing the method on a more practical physical model for further exploration of undulatory mechanisms.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.