Thrust Force Generated by Heaving Motion of a Plate: The Role of Vortex-Induced Force

To understand the force acting on birds, insects, and fish, we take heaving motion as a simple example. This motion might deviate from the real one. However, since the mechanism of force generation is the vortex shedding due to the motion of an object, the heaving motion is important for understanding the force generated by unsteady motion. The vortices released from the object are closely related to the motion characteristics. To understand the force acting on an object, information about momentum change is necessary. However, in vortex systems, it is impossible to estimate the usual momentum. Instead of the momentum, the “ virtual momentum, ” or the impulse, is needed to generate the force. For calculating the virtual momentum, we traced all vortices over a whole period, which was carried out by using the vortex-element method. The force was then calculated based on the information on the vortices. We derived the thrust coefficient as a function of the ratio of the heaving to travelling velocity.


Introduction
Motion of insects or birds is inherently unsteady. The creatures utilise the unsteadiness efficiently. For example, a coherent structure called the leading edge vortex (LEV) plays an essential role in the generation of unsteady force. Many authors have published studies on the topic and hilighted its importance, experimentally and numerically. The magnitude of the unsteady force cannot be explained by a steady-state approach. In many cases, the unsteadiness generates greater forces more efficiently than that in the steady state [1,2]. Experiments have been conducted in three-dimensional space and numerical analyses have been carried out to understand the mechanism of force generation. These studies explained several aspects of unsteady phenomenon, but the role of vortices generated close to the object is still unclear. How does the behaviour of vortices affect the generation of force? In particular, how does momentum change depend on the force? We are not sure how to estimate the momentum of a vortex system, because the usual momentum has no definite value. Our aim is to establish a rule that governs the force generation by the momentum change. Characteristics such as the magnitude, the rotation direction, and the position are key to determining the momentum.
Unless we determine their properties, the evaluation of force cannot be made quantitatively.
When an object of a constant circulation Γ moves with a constant speed dx 0 =dt, a fluid force acts perpendicular to the direction of motion. The magnitude is known to be ρ dx 0 =dt ð ÞΓ. It should be noted that the magnitude is the derivative of the virtual momentum ρx 0 Γ with respect to time, see [3], Art.157. Here, ρ is the density of fluid. This is a simplest application of a well-known law that governs the conservation of virtual momentum. In other words, this is a typical example of the second law of motion in the vortex motion. In general, the virtual momentum plays an essential role in the generation of force instead of the normal momentum. As illustrated above, in unsteady flows, the virtual momentum is important for the generation of force. We would like to illustrate the role of the virtual momentum by applying it to a heaving motion of a thin plate.
A lot of attention has been paid to the dependence of parameters characterising the unsteadiness known as the reduced frequency or the Strouhal number of the propulsive motion of insects, fish and humans (for example, [4][5][6]). Here, we also discuss the dependence of the reduced frequency on the thrust.
The heaving motion of a thin plate is the simplest and most suitable example of the analysis of unsteady phenomena. In addition, the heaving motion is solved in the limit as the heaving amplitude becomes smaller. For investigating the unsteady phenomenon, the vortex motion is a key concept. The analytical tool used here is not specific and can be extended to wider problems.

Direct effect of a heaving plate
First, we have a look at the relation between the force acting on a body fixed in a stream and the free vortices flowing behind it. It is known that a drag acts on a still body set in the stream. We can see two vortex rows here, called the Kármán vortex street (see Figure 1(a)).
We can also notice another similar vortex street behind the flying birds and the swimming fish. However, the direction of rotation of the vortices is inverse. In the case of the Kármán street, a momentum defect is observed while the momentum seems to increase behind the birds and fish. In the latter case, a thrust acts on the object to move forward due to the increase in momentum. As an example, we show the vortex street appearing in heaving motion (see Figure 1(b)). In pitching motion, a similar street can be observed (see example, [7]). In general, those cases where backward momentum increases generate thrust acting against the flow. In the figure, the thick arrows denote the direction of the increased momentum.
To understand the mechanism of thrust generation we study the heaving motion of a thin plate in a uniform flow. We assume that the plate has a constant circulation Γ. Even in unsteady conditions, we assume that the fluid flows smoothly at the trailing edge according to Kutta's condition. The circulation Γ is determined by this smoothness condition. The velocity around the leading edge would diverge and hence the pressure may not be finite because the edge is a mathematical singular point.
To evaluate the force acting on an object, we usually integrate the pressure on the surface of the object. However, because a simple plate has two singular points at the leading and trailing edges. In particular, the estimation of the pressure at the leading edge is almost impossible when Kutta's condition is applied at the trailing edge. Instead of the integration of pressure, we apply Newton's second law of motion, which states that the force is a result of the momentum change. However, it is known that the estimation of momentum is almost impossible, and hence virtual momentum has to be used instead.

Effect of bound vortex
The coordinates system is shown in Figure 2. A thin aerofoil is located at z ¼ z 0 t ð Þ in the complex z-plane or at z 0 ¼ 0. The coordinates z and z 0 are related by the equation Consider a uniform flow whose velocity is U in the x-direction and a bound vortex of a constant circulation Γ around the plate and no free vortices. The circulation is positive when the fluid rotates in the anticlockwise direction, while the vorticity is positive for vortices rotating in the clockwise direction. The force X þ iY acting on the object located at z 0 ¼ x 0 þ iy 0 À Á is given by: where the dot denotes the derivative with respect to time t [8]. Here, the length of the plate is 4a (=L) and located parallel to the uniform flow (see Figure 2).
Confining ourselves to the oscillation only in the y-direction, or _ z 0 ¼ i_ y 0 t ð Þ, the force can be: For cases without any motion, the above equation is written simply as Y ¼ ÀρUΓ, which corresponds to the lift known as the Kutta-Joukowski theorem.
The second term on the right-hand side indicates the drag defined as where m 0 is called the virtual mass. The direction of this force is parallel to the direction of motion. Accordingly, this force which acts in the y-direction cannot contribute to the propulsion. The virtual mass for this thin plate is expressed as πρ L=2 ð Þ 2 ¼ m 0 ð Þ(see [9], Art. 9.222, [10][11][12] for the general discussion). This force acting only in the y-direction is independent of vortex formation and shedding. The force is not related to thrust, and hence we will not discuss this force any more. Finally, from Eq.(3) the force in the x-direction is This formula corresponds to the Kutta-Joukowski theorem. When the object with the circulation Γ is located at z ¼ z 1 , the virtual momentum is expressed as Àiρz 1 Γ.
Eq. (5) can be derived easily by considering the virtual momentum. For an object with a constant circulation Γ 1 located at the position z 1 , the momentum, or more precisely the virtual momentum, P, of the flow is expressed as Àiρz 1 Γ 1 . When the vortex moves at the speed _ z 1 , the force F acts on it as a result of momentum change, i.e.,

Effect of free vortex
Next, we proceed to discuss about the effect of free vortices on the force. The general rule for estimating the force, when the viscosity is negligible, is the Blasius formula, see [10]. Since the formula is valid only for steady flow conditions, it has to be extended to include the unsteady effect. The extended formula for the force X, Y ð Þ, as seen in, for example, [13], is given as where B denotes the path along the surface of an object in the anticlockwise direction. In the above equation, f z ð Þ is the complex potential defined by f z ð Þ ¼ ϕ x, y ð Þþiψ x, y ð Þ. Here ϕ and ψ are the velocity potential and the stream function, respectively. The bar denotes the complex conjugate. Thisin Figure 10 formula expresses two typical types of forces. One is the virtual momentum (VM) component, and the other the direct-interaction (DI) component. VM acts due to the change in momentum and DI is the direct interaction of the vortices with the body, which becomes important when the vortex is near the object. We denote the two forces F v for VM and F d for DI to distinguish between them. Before we discuss the general case, we consider a simple one where one free vortex κ 1 exists at z ¼ z 1 . The forces for VM and DI are expressed as where df =dz ½ c means the convection velocity at z ¼ z 1 by the vortex κ 1 . On the other hand, the force for DI is estimated from First, we consider Eq.(10). This force is dependent on the object form. To integrate it we map a plate in the z-plane to a circle of radius a in the ζ-plane as When a vortex is located at z ¼ z 1 , the integration can be carried out to give where and the convection velocity, It is easy to see that the right-hand side of Eq. (12) is pure imaginary, because the right-hand side expresses the sum of a complex and its complex conjugate. This means that the force has only a y-component. Therefore, the component F d is not related to the thrust. Hence, we will not discuss F d anymore. Only the VM would contribute to the thrust force.
From Eq.(10), we have where _ z 1 is the covection velocity of vortex κ 1 . The above equation is to Eq.(6), because 2πκ 1 =ÀΓ 1 . To determine the convection velocity _ z 1 , we apply the conformal mapping Eq.(13) and trace the vortex in the ζ-plane and then calculate the velocity in the z-plane. The moving speed of vortex κ 1 in the z-plane is already given by Eq. (13). Hence, we have Formulas (5) and (15) are the main targets for the calculation of thrust.

Determination of positions and velocities of a vortex
Now, we discuss how to generate a vortex under our boundary condition. What determines the vorticity and its position? Consider a flat plate set parallel to the flow (see Figure 2). Even in unsteady motion, the flow is subject to the condition that the fluid flows smoothly at the trailing edge. In other words, Kutta's condition at the edge must be satisfied at all times. We consider the heaving motion whose velocity, perpendicular to the plate is expressed as In the above equation, ν is the radian frequency of the heaving motion, and W T is the amplitude. Denoting the period of the oscillation as T, T ¼ 2π=ν.
Because the plate has a velocity in the y-direction at t ¼ 0, Kutta's condition is not satisfied. To satisfy the condition we set a new vortex at x ¼ 2a þ Δx, and we determine the vorticity κ 1 of the vortex so as to satisfy the condition. As for setting the initial position, [14] serves as a useful reference. The condition for the flow leaving the trailing edge smoothly determines κ 1 uniquely. Later at t ¼ Δt the vortex κ 1 moves away and hence the flow does not satisfy Kutta's condition again. To avoid the undesirable flow, we set a new vortex κ 2 at the same position as the initial position of κ 1 , i.e., at x ¼ 2a þ Δx. Kutta's condition fixes the value κ 2 uniquely. Similarly, the subsequent process determines sequentially κ i i ¼ 1, 2, … ð Þ . We proceed to the next step to discuss the problem of movement of vortices. A vortex moves by the other free vortices including the bound vortex and the uniform velocity. The induced velocity w ¼ u À iv at z by the vortex κ c located at z ¼ z c is written as: Actual calculations were done in the ζ-plane with respect to all the vortices including those of the mirror image. The calculation step was carried out at every time for the step Δt. See [14] for the suitable relation between Δx and Δt.

Calculation results
In the calculations, we determine the physical variables by choosing a ¼ 1, ρ ¼ 1 and U ¼ 1. In Figure 3, we show their positions and the direction of rotation for the case when ν ¼ 0:5 and W T ¼ 0:5 at t ¼ 19:8. The symbol + denotes the vortices of the clockwise rotation and those of the triangle (in red) the anticlockwise one, respectively. It is seen that the vortices rotating in clockwise direction gather at some places in the negative y-plane, while those rotating in anticlockwise direction gather in the positive y-plane. Figure 4 shows the positions at t=39.8. We can find three clusters of vortices of clockwise rotation at about x ¼7.5, 20, and 33, and three clusters of anticlockwise rotation at x ¼13, 27, and 37. The clusters of positive or negative vortices occur by the interaction of each vortex. At those positions, the vorticities concentrate and have a structure in a large scale. Nonlinearity is seen even for such low W T (=0.5). Three clusters of vortices rotating in the clockwise direction are in the area for y < 0, while three clusters rotating in the anticlockwise direction are in the area for y > 0. This array of two vortex streets would generate the downward flow, which suggests that the momentum is generated in the positive x-direction. Momentum generation in the positive x-direction means the generation of thrust force, as will be explained later. The deviation of arrays from the ordered ones is the result of nonlinearity. Figure 4 also shows the deviation of sinusoidal distribution of vortices. Next, we consider the positions of vortices at initial stages near t=0. Those vortices generated initially, which are distributed near x ¼ 40, fluctuate violently and move to the positive y-direction.
In Figure 5, the distribution of vortices κ i determined in the manner explained earlier is depicted. This plot shows the complex distribution of vortices based on the interactions among many vortices. This may explain the reason why the clusters are generated.

Calculation of force
In the following section, we describe calculations carried out when U=ρ=a=1 unless specified otherwise.

Direct force by movement of a plate with a circulation
According to Eq.(5), the movement in the y-direction of the plate with a circulation Γ 1 gives rise to the force X b normal to it, We investigate the thrust generation due to the movement of a thin flat plate in more detail. When the motion is subjected to Eq.(17), we consider the force in the y-direction at the initial stage t ≈ 0. Vortices rotating in the positive direction appear under the lower place near the trailing edge. Similarly, in the ζ-plane, the mirror images of the vortices rotating in the anticlockwise appear in a circle with radius a. In these, circumstances the circulation Γ around the circle is positive. In this case the force acts in the negative x-direction, because the sign of X b is negative from Eq.(18).
The case where two free vortices are outside the circle is shown in Figure 6. For more than a vortex in the flow field there must be mirror images whose sign is opposite to the free vortices. In general, at time t, many vortices of the same number of free vortices exist inside the circle.
When n vortices are released, the intensity of the vortices is expressed as a sum P n i¼1 κ i . At the same time, the sum of vortices within the circle of the radius a determines the circulation Γ of the bound vortex. The circulation of the bound vortex is expressed as Using this circulation, we try to evaluate the force generated by the heaving motion. In Eq.(18) by changing Γ 1 by Γ, we have the force, The variation of X b calculated by using the above equation is shown in Figure 7 as a function of nondimensional time τ ¼ t= 2a=U ð Þ. The variation of the heaving velocity w h τ ð Þ of the plate (Eq.(17) is also plotted there). Let us consider the initial stage when the plate moves upward. Vortices generated by the upward movement are rotating in the clockwise direction, as shown in Figures 4 or 5. At the initial stage, mirror images inside the circle of the radius a rotate in the anticlockwise direction. In other words, the plate has a positive circulation. Because the velocity w h is initially positive, the force X b is negative from Eq.(20). The force acts against the main flow, i.e., the plate is pulled by the fluid in the negative x-direction. Usually, the upward motion connects with positive circulation, and hence the force becomes almost negative. On the contrary, negative circulation occurs when the motion is downwards. As a result, the sign of X b has a negative value in the mean. The average value X b av ð Þ is À0.240. The index (av) stands for the mean over two periods of oscillation. It should be noted that when the absolute velocity of the plate |w h | is the maximum, the force becomes maximum. This means that during the generation of strong vortices, the pressure at the edge becomes large. Strictly speaking, slight time delay is also observed. This may be because of the effect of the convection due to other free vortices.
The  circulation must be taken into account. From this point of view, there is a room for reconsidering the results. As seen in Figure 7, the force X b varies with a period π=ν ¼ T=2 ð Þand has a negative value on an average. However, we did not take into account the variation of Γ. The circulation Γ around the plate changes with the same period π=ν. In Eq. (20) we took into account the differentiation of the vitual momentum partly, and it could not give the correct force induced by virtual momentum. The x-component of the real virtual momentum, ρΓ t ð Þy 0 t ð Þ, has two time-dependent variables, Γ and y 0 . To estimate the correct force X b , we should take into account the variation of the virtual momentum. The correct expression for the force: In the present situation, Γ and y 0 are both periodic function of time whose period is 2π=ν ¼ T ð Þ. The product of two periodic functions with the same period is also a periodic function. The differentiation with respect to t is also a periodic function whose average is zero. Finally, we conclude that the force X b gives no net force, or X b av ð Þ =0. Here the subscript (av) stands for the average over two periods of time, 2T.

Effect of moving vortices
In this subsection, we discuss the force resulting from the movement of free vortices. First, we show the result of the force in the y-direction. This problem was first discussed and the solution was analytically given by Kármán-Sears in the linear limit [15]. Their solution corresponds to the sum of the forces Y v and Y d . The force Y v has already been given in Eq. (15) only when one free vortex exists. For the present aim, however, the formula should be extended to include all the vortices. In the following, according to [15] the variation of force divided by 2aρU is shown. When W T =0.5 and ν=0.5 the variations are given in Figure 8. The variation of Y v in the VM, and that of the sum of Y v and Y d in the DI. The change of the sum Y v +Y d agrees well with the analytical result of [15]. In particular, the agreement becomes better for a lower W T . It is seen that the two components Y v and Y d have an importance of the same degree on the generation of force. At an initial stage, τ ≈ 0, the sum has a negative value, and the minimum value of the force occurs at the stage where the velocity in the y-direction becomes maximum, which corresponds to the initial instant τ ¼ 0. The force acts as a drag in this heaving motion. It is interesting to investigate the y-component of the force with respect to the virtual momentum. Because such a force in the y-direction is not related directly to the thrust force, hereafter we will not discuss it further. At the initial stage, it is seen from Figures 4 and 5 that vortices of positive vorticity appear. These vortices travel to the position near x ¼ 40 at τ ≈ 20.
Next we consider the thrust component of the force generated by the change of virtual momentum. From Eq.(15) the force component is expressed for a vortex κ 1 as Taking into account all the vortices existing in the flow field, we can get a complete set of the component for the present problem. The variation is shown in Figure 9. It seems to oscillate sinusoidally except for the initial stage and has a positive value in the mean. For this example, the mean value X v av ð Þ is calculated to be 0.148. This means that the force acts in the positive x-direction or the fluid pushes the plate to the direction of the flow. Similarly, the plate is adding the force to the flow as a reaction. In this sense, we may regard the positive X v av ð Þ as a thrust. Whether this force acts as a thrust or a drag depends on the combination of the sign of κ 1 and that of the velocity _ y 0 . Possible combinations are listed in Table 1. Behind the heaving plate there appear two vortex streets, as shown in  the vortices rotating in the positive direction move upward and those rotating in the negative direction moves downward. This tendency is pronounced for the vortices existing near the trailing edge. It is noted that the force has its peak when the heaving velocity w h is at the maximum or the minimum. The period of the force oscillation is T=2. When the plate passes through y ¼ 0, the vortex with a strong intensity is generated. At this instant, the force reaches the maximum. Table 1 suggests that the sign of the force X v is positive. In fact, it is seen from Figure 9 that the average of the force X v is calculated to be positive.
By comparing Figures 7 and 9 it is clear that the force component X v is small compared to X b . The heaving motion has an influence more effectively on the generation of force in the y-direction As far as the thrust force is concerned, however, the force X b produced directly by the heaving motion has no effect. Therefore, the force we should take into account is the force X v only as a thrust.

Effect of heaving amplitude on the force
It seems that the thrust force is generated due to the motion of the plate against the fluid. To understand the role of the heaving amplitude W T , we plotted X v av ð Þ as Force in the x-direction. a function of W T in Figure 10. As mentioned earlier, the subscript (av) means the average over two periods. It is easy to see that the thrust is proportional to W 2 T except when the W T value is lage. In this case, the proportional constant is estimated to be 0.592. In addition to X v av ð Þ , the variation of X b av ð Þ is also plotted for comparison.
Next, we show the variation of the thrust X v av ð Þ =2aρU as a function of U in Figure 11. The curve seems to be inversely proportional to U except for large values of U. This means that the thrust X v does not depend on the velocity U. When W T =0.5, the constant of proportionality is estimated as 0.148.

Effect of heaving frequency on the force
From the previous subsection, it can be seen that X v is proportional to W 2 T and does not depend on U. To confirm it, we have plotted the nondimensional variable X v av ð Þ = L=2 ð Þ=ρU 2 as a function of W T =U ð Þ 2 in Figure 12 for three different W T 's, i.e., 0.3, 0.5 and 0.7. The coefficient X v av ð Þ =2aρU 2 is called the thrust coefficient denoted as C L . The coefficient has almost linear relation to the velocity ratio W T =U. The constant of proportionality must be nondimensional. In such unsteady locomotions, the most important dimensionless parameter is the reduced frequency k ¼ 2aρν=U ð Þor the Froude number. However, Figure 12 gives no defined dependence of the reduced frequency on the coefficient C T . The plotted data include various values of k between 0:2 ≤ k ≤ 2:5. We can draw our conclusion that the reduced frequency k does not affect the thrust coefficient in this heaving motion.
We summarise the thrust coefficient in the nondimensional form, where L is the chord length being equal to 4a.

Concluding remarks
Thrust force can be generated by a simple heaving motion of a plate. The force is perpendicular to the direction of oscillation. A pair of rows of vortices plays an important role in the generation of the force. The two vortex streets give rise to an increase in momentum in the direction normal to the direction of oscillation. The  word "momentum" here does not mean the usual momentum but the virtual one, because the usual momentum cannot be determined in such a vortex system. The direct integration of the pressure around the surface of a body is not a correct way to know the thrust generation. Application of the virtual momentum to the generation of force made the estimation of the force possible.
In general, the most important parameter characterising the unsteady flow is the reduced frequency k or Froude number. How the parameter plays a part in the generation of force has been a main concern of many people. Many researchers have tried to address this problem experimentally. However, the task was difficult to address, and only few researchers have answered analytically.
Our result is for the coefficient of thrust C T , The proportional constant is nondimensional and does not depend on the parameter k expressing the unsteadiness of flow. The thrust force X v av ð Þ is independent of uniform velocity U, and therefore the coefficient C T is proportional to W T =U ð Þ 2 . Although our analysis is confined to the heaving motion of a thin plate, we summarise that the force due to the vortex movement can be expressed as a function of nondimensional quantity in a simple form. It is expected that our analysis could apply to more complex movement of an aerofoil.