Research on the Undulatory Motion Mechanism of Seahorse Based on Dynamic Mesh

The seahorse relies on the undulatory motion of the dorsal fin to generate thrust, which makes it possess quite high maneuverability and efficiency, and due to its low volume of the dorsal fin, it is conducive to the study of miniaturization of the driving mechanism. This paper carried out a study on the undulatory motion mechanism of the seahorse's dorsal fin and proposed a dynamic model of the interaction between the seahorse's dorsal fin and seawater based on the hydrodynamic properties of seawater and the theory of fluid-structure coupling. A simulation model was established using the Fluent software, and the 3D fluid dynamic mesh was used to study the undulatory motion mechanism of the seahorse's dorsal fin. The effect of the swing frequency, amplitude, and wavelength of the seahorse's dorsal fin on its propulsion performance was studied. On this basis, an optimized design method was used to design a bionic seahorse's dorsal fin undulatory motion mechanism. The paper has important guiding significance for the research and miniaturization of new underwater vehicles.


Introduction
With the increasing demand for natural resources in modern society, the speed of exploitation of terrestrial resources is difficult to meet people's needs for material life. The ocean, which accounts for 71% of the entire surface of the earth, has become a treasure trove of resources for all countries. There are not only abundant fishery resources and mineral resources but also sufficient energy resources, such as large oil fields and combustible ice. Whether it is economic or military, the treasure of the ocean is attractive enough, and the rapid development of underwater vehicles will become inevitable. The seahorse relies on the undulatory motion of the dorsal fin to generate thrust, which makes it possess quite high maneuverability and efficiency, and due to its low volume of the dorsal fin, it is conducive to the study of miniaturization of the driving mechanism. This article has carried out research on the undulatory motion mechanism of the seahorse's dorsal fin.
The current research on the seahorse mainly focuses on the following aspects: some focus on the distribution and breeding of seahorse populations in global seas [1][2][3]; some focus on the evolution of the seahorse, as well as basic biological characteristics and living habits [4,5]; others focus on the kinematics and dynamics of the seahorse's dorsal and caudal fins [6][7][8][9][10]. In addition, a high-speed camera system is used to study the undulatory motion of the seahorse's dorsal fin [11]; and also, some focus on the study of the physiological structure and the mechanical properties of the muscle of the seahorse's dorsal fin [12].
So far, few studies have been conducted on the undulatory motion mechanism of the seahorse's dorsal fin. The reason is mainly due to the uncertainty of living seahorse's movement, which leads to considerable difficulties in the setup of the experimental device, and this uncontrollable movement will also make it difficult for high-speed cameras to obtain sufficient illumination and focus. In addition, differences between seahorse's species, between different genders, and between different individuals will hinder systematic and reproducible research. There are so many difficulties in living animal experiments, so many scholars usually manufacture bionic prototypes for mechanism research. This method can solve the above-mentioned problems, but the prototype dorsal fin is difficult to achieve the swing amplitude and frequency like a living seahorse, and the difficulty in manufacturing and control causes the experimental results to be inaccurate. In addition to physical experiments, simulation can also be used to study the mechanism of seahorse's movement. However, due to insufficient computer computing power and other reasons, early simulations were mainly two-dimensional plane simulations, which were difficult to directly explore the movement mechanism. Nowadays, the computing power is greatly improved. Compared with physical experiments, fluid-structure coupling simulation can improve the efficiency of experiments and reduce the cost of making prototypes. At the same time, it can reduce the impact of the uncontrollable motion of the living seahorse on the experiment, and it is beneficial to the systematic and repeatable research [13].
Based on the theory of fluid-structure interaction, this paper uses Fluent software to construct a dynamic model of the interaction between the seahorse dorsal fin and seawater. The influence of the different swing frequency, wavelength, and amplitude of the dorsal fin on the undulatory motion of the seahorse's dorsal fin was analyzed, which provide important support for the research on the undulatory motion mechanism of the seahorse's dorsal fin. On this basis, an optimal design method was used to design a seahorse-like dorsal fin wave motion mechanism, which has important guiding significance for the development of new underwater vehicle research and the miniaturization of the vehicle.

Physical
Model. The seahorse's dorsal fin is composed of fin rays and fin membrane. The length of adult seahorse's dorsal fin is generally between 3 and 25 mm, and the number of fin rays is between 10 and 30. Some types of seahorse are shown in Figure 1 [14], and the number of fins and dorsal fin length of some types of seahorse are shown in Table 1.
The common perpendicular of the fins is defined as the z -axis, and this direction is called the chord direction. The y -axis is perpendicular to the z-axis and straight down, and this direction is called the span direction. Finally, according to the right-hand spiral, determine the x-axis direction. The   2 Applied Bionics and Biomechanics coordinate system is shown in Figure 2. Points on the dorsal fin of the seahorse approximately rotate around the z-axis.
For the real seahorse's dorsal fin swinging, its swing amplitude should change with the z-axis position. In order to simplify the model, it is considered that the swing amplitude does not change with the z-axis position. The kinematics equation of the dorsal fin is where l is the fin length, L is the total length of the dorsal fin, A is the swing amplitude, f is the swing frequency, λ is the swing wavelength, and t is time.
First, calculate the coordinates of a large number of points on the dorsal fin in excel. Then, import them to Solidworks to generate a point cloud, as shown in Figure 3(a). Through surface treatment, the physical model of the dorsal fin is obtained as shown in Figure 3(b).

Fluid-Structure Coupling Dynamics
Modelling. In order to obtain the dynamic model of the seahorse's dorsal fin, the dynamic equation of the seahorse's dorsal fin needs to be established first. First, establish two coordinate systems, one of which is fixedly connected to the seahorse's dorsal fin and is called the coordinate system O ′ . And the other is an inertial coordinate system which is called the coordinate system O. The two coordinate systems coincide at the initial moment.
The position of the origin of the coordinate system O ′ is given by: where a, b, c is the position of the origin of the coordinate system O′ in the coordinate system O. a 0 , b 0 , c 0 is the initial position of the origin of the coordinate system O′ in the coordinate system O. v x , v y , v z is the speed of origin of the coordinate system O′ along the x-, y-, and z-axes in the coordinate system O.

Applied Bionics and Biomechanics
The relationship between Euler angle change rate and angular velocity is given by: Then, we can get the Euler angle by: where α, ϕ, θ is the Euler angle of the coordinate system O ′ . α 0 , ϕ 0 , θ 0 is the initial Euler angle of the coordinate system O′. ω x , ω y , ω z is the angular velocity of the coordinate system O ′ .
The speed in formula (2) is obtained by: where v ! 0 is the initial velocity of the origin of the coordinate system O ′ in the coordinate system O. The angular velocity in formula (4) is solved by where ω ! 0 is the initial angular velocity of the origin of the coordinate system O ′ .
The position of any point on the dorsal fin in the coordinate system O is given by: where HðtÞ is the coordinate transformation matrix. r where f x , f y , and f z are the mass force components of unit mass fluid in x, y, and z directions. v x , v y , and v z are the velocity components of the fluid in the x, y, and z directions. p is relative pressure. v is the kinematic viscosity of the fluid.
The continuity equation for viscous fluid is given by: There are three existing turbulence numerical simulation methods: Direct Numerical Simulation (DNS), Reynolds Average Navier-Stokes(RANS), and Large Eddy Simulation (LES). RANS is the application of statistical theory of turbulence, which is the simulation method commonly used in engineering. Usually based on Boussinesq's eddy viscosity hypothesis, the zero equation, one equation, or two equations are introduced to close the equation. The zeroequation model has a common shortcoming, that is, the turbulence viscosity coefficient only depends on the local flow  Applied Bionics and Biomechanics parameters, and has nothing to do with the flow elsewhere, which is inconsistent with experimental observations. On this basis, a one-equation model and two-equation model were developed. Among the two-equation model, the SST model has certain accuracy and consumes limited computing resources and has higher calculation accuracy for nearwall surfaces compared to the other two-equation models. Therefore, the SST model is used in this project.

Simulation Model.
Use the grid processing software ICEM to mesh the computing space. The area near the fin surface needs to be focused, so the density of mesh nodes near the fin surface increases. The height of the first cell perpendicular to the body surface is 0.05 mm. This height is chosen to make the y + of most of the cells in contact with the body surface fall within the effective range of the standard wall function. At the same time, considering the calculation speed, take the outer division step of the model as 0.1 mm, and the model is shown in Figure 4. Import the ICEM file into Fluent, and set the seawater parameters according to the seahorse's living environment, and then, you can simulate the undulatory motion process of the seahorse's dorsal fin under different conditions. Set the model boundary parameters as shown in Table 2, and the Fluent simulation preset parameters as shown in Table 3.
The undulatory motion of the seahorse's dorsal fin is controlled by using Fluent UDF. In the numerical calculation process, the instantaneous force acting on the dorsal fin is obtained by integrating the fin surface pressure and shear stress.

Results and Discussion
3.1. Simulation Process. Import the Fluent calculation results into CFD POST for postprocessing, and then, we can get the flow field distribution at any time. Take 1 Hz frequency, 200 mm wavelength, and π/5 swing amplitude for qualitative explanation. The pressure distribution on the surface of the dorsal fin in 0.25 s, 0.5 s, 0.75 s, and 1 s is shown in Figure 5.
It can be clearly seen from Figure 5 that one side of the surface of the dorsal fin that pushes the water flow is a high-pressure area, and the other side is a low-pressure area.
Take a section as shown in Figure 6 to obtain the flow field pressure distribution of the section, as shown in Figure 7.
It can be seen more clearly from Figure 7 that one side of the fin along the wave propagation direction is a high-pressure area, and the other side is a low-pressure area. The pressure difference between the two sides results in the generation of forces in the x and z directions. Since the dorsal fin of this example is composed of two complete sine waves, the forces generated in the x direction cancel each other out. Take another section as shown in Figure 8 to obtain the flow field pressure distribution of the section, as shown in Figure 9.
It can be seen from Figure 9 that when the dorsal fin ray swings, due to the pressure difference between the upper and lower fins, a force in the negative direction of the y-axis is generated.

The Effect of Dorsal Fin Swing Frequency on Propulsion.
Use Fluent to simulate the five swing frequencies of 1 Hz, 10 Hz, 35 Hz, 50 Hz, and 100 Hz. The wavelengths are all 200 mm, and the swing amplitudes are all π/5. Record the force of the dorsal fin in x, y, and z directions, respectively, and perform curve fitting in matlab. Since the force scatter diagram has not stabilized in the first few periods, the third        Applied Bionics and Biomechanics and fourth period scatter diagrams are selected for fitting.
Since the minimum frequency of 1 Hz and the maximum frequency of 100 Hz are too far apart, the abscissa is set to time multiplied by the frequency of the corresponding working condition, in order to show the difference of different working conditions more clearly and intuitively. By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing frequencies in the x direction, the following results are obtained.
It can be clearly seen from Figure 10(a) that as the frequency increases, the average force d and the fluctuation amplitude a in the x direction both increase. It can be seen from Table 4 that the approximate sinusoidal frequency of the force in the x direction is basically the same as the swing frequency of the dorsal fin, and the phase is also basically the same. From Figure 10(b), the fluctuation amplitude a and frequency f can be better fitted with a quadratic function, and the relationship between the amplitude and frequency of the force fluctuation in the x direction of the dorsal fin can be obtained by: However, in Figure 10(c), the average force d does not change much after 35 Hz.
By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing   It can be clearly seen from Figure 11(a) that as the frequency increases, the average force d and the fluctuation amplitude a in the y direction both increase significantly. It can be seen from Table 5 that the approximate sinusoidal frequency of the force in the y direction is twice the swing frequency of the dorsal fin, and the phase difference is approximately π/2. From Figure 11(b), the amplitude of fluctuation a and frequency f can be better fitted with a quadratic function, and the relationship between the amplitude of the force fluctuation in the y direction of the dorsal fin and the frequency can be obtained by: From Figure 11(c), the average force d and frequency f can also be better fitted with a quadratic function, and the relationship between the average force in the y direction of the dorsal fin and the frequency can be obtained by:    Applied Bionics and Biomechanics By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing frequencies in the z direction, the following results are obtained.
It can be clearly seen from Figure 12(a) that as the frequency increases, the average force d and the fluctuation amplitude a in the z direction both increase significantly. It can be seen from Table 6 that the approximate sinusoidal frequency of the force in the z direction is twice the swing frequency of the dorsal fin, and the phase is basically the same. From Figure 12(b), the fluctuation amplitude a and frequency f can be better fitted with a quadratic function, and the relationship between the amplitude of the force fluctuation in the z direction of the dorsal fin and the frequency can be obtained by From Figure 12(c), the mean force d and frequency f can also be better fitted with a quadratic function, and the relationship between the mean force in the z direction of the dorsal fin and the frequency can be obtained by: In summary, the average force and fluctuation amplitude in each direction increase with the increase of frequency, and the average force and fluctuation amplitude in the y and z directions have a quadratic relationship with frequency. The amplitude of the force fluctuation in the x direction also has a quadratic relationship with the frequency, but the   [15], the authors presented an experimental investigation of flexible panels actuated with heave oscillations at their leading edge. Results were presented from kinematic video analysis, particle image velocimetry, and direct force measurements. They draw the conclusion "The magnitudes of both signals (net thrust and power) increase with heaving frequency, as expected." And they gave the net thrust and power curves of different frequencies as Figure 13. It can be seen that as the frequency increases, the net thrust gradually increases, and it becomes an approximate sine function in the swing period, which verifies the correctness of the simulation in this study.

The Influence of the Swing Wavelength of the Dorsal Fin on the Propulsion Force.
Since the total length of the dorsal fin is 400 mm, in order to minimize the force fluctuations in the x direction, an integer number of traveling waves should be selected for the entire dorsal fin. Therefore, Fluent is used to simulate the three swing frequencies of 400 mm, 200 mm, and 133 mm, respectively. The frequencies are all 10 Hz, and the swing amplitudes are all π/5. Record the force of the dorsal fin in x, y, and z directions, respectively, and perform curve fitting in matlab. Since the force scatter diagram has not stabilized in the first few periods, the third and fourth period scatter diagrams are selected for fitting.
By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing wavelengths in the x direction, the following results are obtained.
It can be clearly seen from Figure 14 and Table 7 that as the wavelength increases, the average force d in the x  10 Applied Bionics and Biomechanics direction increases. When the swing wavelength is 400 mm, reducing the wavelength can significantly suppress the fluctuation amplitude a of the force in the x direction. But after the wavelength of 200 mm, the influence of decreasing the wavelength on the amplitude a of the force fluctuation in the x direction is greatly reduced. The approximate sinusoidal frequency of the force in the x direction is basically the same as the swing frequency of the dorsal fin, and the phase changes with wavelength. By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing wavelengths in the y direction, the following results are obtained.
It can be seen from Figure 15 and Table 8 that the average force d in the y direction increases significantly with the increase of the swing wavelength, but the fluctuation amplitude a hardly changes.
By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing  It can be seen from Figure 16 and Table 9 that the mean value d of the force in the z direction increases significantly with the increase of the swing wavelength, but the fluctuation amplitude a hardly changes.
Based on the analysis of the force in the above three directions, it can be seen that when the wavelength is between 400 mm and 200 mm, as the wavelength decreases, the x direction fluctuation of the dorsal fin is significantly suppressed while below 200 mm, the impact is small. At the same time, as the wavelength increases, the mean value of the force in the y and z directions increases significantly, but the fluctuation range is almost unchanged.
3.4. The Influence of the Swing Amplitude of the Dorsal Fin on the Propulsion Force. Use Fluent to simulate the three swing amplitudes, the frequency is 10 Hz, and the wavelength is 200 mm. Record the force of the dorsal fin in x, y, and z directions, respectively, and perform curve fitting in matlab. Since the force scatter diagram has not stabilized in the first few periods, the third and fourth period scatter diagrams are selected for fitting.
By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing amplitudes in the x direction, the following results are obtained.
From Figure 17 and Table 10, it can be seen that the swing amplitude has little effect on the force in the x direction. No matter how the swing amplitude changes, the average force d in the x direction fluctuates near the 0 line, and the fluctuation amplitude a is relatively small. The approximate sinusoidal frequency of the force in the x direction is basically the same as the swing frequency of the dorsal fin, and the phase is also basically the same.
By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing amplitudes in the y direction, the following results are obtained.
It can be seen from Figure 18 that as the swing amplitude increases, the average force d and the fluctuation amplitude a in the y direction both increase. It can be seen from Table 11 that the approximate sinusoidal frequency of the force in the y direction is twice the swing frequency of the dorsal fin, and the phase difference is π/2. By using a custom function f ðxÞ = a sin ð2πbx + cÞ + d to curve-fit the force data of the dorsal fin at different swing amplitude in the z direction, the following results are obtained.  It can be seen from Figure 19 that as the swing amplitude increases, the average force d in the z direction and the fluctuation amplitude a both increase. It can be seen from Table 12 that the approximate sinusoidal frequency of the force in the z direction is twice the swing frequency of the dorsal fin, and the phase is the same.
In summary, the average value of the force in each direction increases with the increase of the swing amplitude, but the influence on the force in the x direction is negligible. At the same time, the increase of the swing will cause the fluctuation of the force to increase.

Design of Undulatory Motion Mechanism Imitating
Seahorse's Dorsal Fin. In order to realize that adjacent fins oscillate with a fixed phase difference, a crank-rocker mechanism driven by a crankshaft is designed, as shown in Figure 20.
A crankshaft is composed of a left part, a right part, and a plurality of middle parts, as shown in Figure 21. The left part has two protruding ends, the left of which is in the cen-ter while the right one with eccentricity. There is a shaft flat position on the right-side extension to match the middle part. An eccentric hole on the left side of the middle part is matched with the left part, and an eccentric shaft on the right side is matched with the next middle part. There is also a flat shaft position, and there is a 45°phase difference between the left hole and the right shaft, which makes adjacent fins produce a fixed phase difference. The right part also has an eccentric hole to fit with the middle part, and a concentric shaft on the right is connected to the motor. When making crankshaft parts, consider hollowing out the middle disc to reduce the moment of inertia.
The comprehensive problem of the function mechanism means that the functional relationship between the input and output angles corresponding to the rocker and the crank is required to be as close as possible to the given function ψ = f ðφÞ. The comprehensive theory of the mechanism proves that the kinematics of the crank-rocker mechanism has nothing to do with the actual length of the rod, but only with the shape of the mechanism, that is, the relative length between the rods. Let us set the relative length of the crankshaft as l 1 = 1. The relative length of the connecting rod is l 2 . The relative length of the pendulum is l 3 . The relative length of the frame is l 4 . And the initial angle of the crankshaft and the swing lever is φ 0 , ψ 0 , so there are 5 variables. Therefore, at most five sets of corresponding angular positions can be accurately met. If the organization is required to best approximate the expected function in more positions, the optimal synthesis method can be used. In the  13 Applied Bionics and Biomechanics optimization design of the kinematics of the planar four-bar linkage, the objective function is generally established according to the kinematics parameters of the mechanism. For example, the movement realized by a four-bar linkage mechanism is an input-output angular function derived from the geometric relationship of the mechanism's movement, and it is required to have the smallest deviation from a given function within a certain range of motion. In the optimization design of linkage mechanism dynamics, it is relatively simple to use the pressure angle and transmission angle in the mechanism as important indicators for the motion analysis and dynamic analysis of the mechanism. In order to obtain good transmission performance and increase the reliability of the mechanism, it is necessary to select the best mechanism dynamics parameters so that the maximum pressure angle is the smallest or the minimum transmission angle is the largest during the movement of the mechanism. In this project, the crank angle is required to be at any position in a circle, that is, when φ = φ 0~ð φ 0 + 2πÞ, the output angle of the rocker is as close to the function ψ = ψ 0 + π/5ðsin ðφ − φ 0 − π/2Þ + 1Þ as possible. Assuming that the initial position angles of the crank and the joystick φ 0 , ψ 0 correspond to the input and output angles when the joystick is in the right extreme position. In this way, φ 0 , ψ 0 is no longer an independent variable, so there are three relative lever length variables left. Since the three-dimensional search is more complicated and time-consuming, it is assumed that the relative length of the rack is l 4 = 15. Take the relative length of the connecting rod l 2 and the relative length of the rocker l 3 as design variables.
The crank-rocker mechanism is in accordance with the corresponding relationship between the input and output angles between the driving crank and the driven rocker. The independent parameters include the relative length of the rod l 2 /l 1 , l 3 /l 1 , and l 4 /l 1 and the initial angular positions of the crank and the rocker φ 0 and ψ 0 . As shown in Figure 22, assume that the acute angle between the crank and the rocker and the frame when the pendulum reaches the right limit position is taken as the initial position angle φ 0 and ψ 0 . Due to the geometric relationship of the right limit position, the crank and the connecting rod are collinear, so the two initial angles can be determined according to the geometric relationship of the initial position.
Therefore, the two initial angles can be expressed by other rod lengths and are no longer independent parameters. According to the above assumption, the crank length is the unit length l 1 = 1. Frame length is l 4 = 15. Therefore, the length of the connecting rod l 2 and the length of the rocker l 3 are selected as design variables X = l 2 . This turns into a two-dimensional optimization design problem. Taking the least square deviation of the output angle of the mechanism as the design goal, the given function and the actual function are discretized, and the discrete deviation function is obtained. The sum of the discrete deviation functions is used as the objective function, as shown in where ψ i is the expected output angle and ψ si is the actual output angle.
According to the given functional relationship and the corresponding relationship between the two initial angles, the output angle expression can be obtained by where s is the uniform number of discrete points of the crank angle φ in the interval φ 0~ð φ 0 + 2πÞ and i is the serial number of each discrete point. The actual output angle expression can be determined according to the geometric relationship of the movement of the mechanism, as shown in Figure 23.
Among them, in ΔBDC and ΔABD, applying the law of cosines, we can get In order to make the transmission performance of the mechanism better, the minimum transmission angle of the mechanism γ min ≥ 45°and the maximum transmission angle of the mechanism γ max ≤ 135°. When the crank and the frame are collinear, the mechanism has the minimum or maximum transmission angle, as shown in Figure 24. When the mechanism is in these two positions, the law of cosines can be used to obtain by:  Figure 19: z-direction force on the dorsal fin at different swing amplitude.
It can be seen that the output angle function of the mechanism fits well with the objective function.

Conclusions
In this paper, a theoretical model of the interaction between the seahorse's dorsal fin and seawater is firstly carried out, and a fluid-structure coupling model suitable for this subject is derived, and then, a simulation model of the seahorse's dorsal fin and seawater is established based on the fluidstructure coupling theory and the physical properties of seawater. And refer to the Fluent database for parameter matching. Taking a certain working condition as an example, the interaction process between the seahorse's dorsal fin and the sea is analyzed, and the pressure distribution on the surface of the seahorse's dorsal fin and the flow field around it during the undulatory motion of the seahorse's dorsal fin is obtained through data postprocessing. Then, using the controlled variable method, keeping other variables the same, the influence of the swing frequency, wavelength, and amplitude of the seahorse's dorsal fin on the instantaneous force and average force in various directions of the seahorse's dorsal fin during the undulatory motion is studied, and the swing frequency, wavelength, and amplitude of the seahorse's dorsal fin are summarized. For the influence of the instantaneous force and the average force on the seahorse's dorsal fin in various directions during undulatory motion, finally, an optimal design method is used to design a seahorse-like dorsal fin undulatory motion mechanism, which has important guiding significance for the development of new underwater vehicle research and the miniaturization of the vehicle.

Data Availability
The data used to support the findings of this study are included within the article.