108, 108′ is moving in positive y-direction (the leading side) and negative while the foil 108, 108′ is moving in negative y-direction (the trailing side). A smooth square wave ssw is proposed for the base shape of the target angle of attack α_target:

α

t
⁢
a
⁢
r
⁢
g
⁢
e
⁢
t

=

A
⁡
(

1
+

B
A

⁢
s
⁢
s
⁢
w

)

⁢
s
⁢
s
⁢
w

,

(
1
)

where ssw is

ssw
=

cos
⁢

(

θ
-

1
2

⁢

sin
⁢

(

2
⁢
θ

)

)

.

This profile with exemplary amplitudes A=14° and B=1° is presented in FIG. 3 using the dashed curve. The use of this asymmetric smooth square wave for side force adjustment may be elaborated in an embodiment (see FIG. 6 and the description related to it). The continuous line in FIG. 3 refers to an example of a symmetric ssw, A=15°, B=0°.

FIG. 4 illustrates an example of the foil 108 rotating at an angular speed ω at an arbitrary foil location expressed as function of the rotation angle θ of the foil wheel 106. The data on the pitch angle γ(θ) of a foil 108, 108′ may be formed such that a target angle of attack α (with chosen amplitude) is attained as closely as possible. Mathematically, the situation can be characterized in a following manner, for example:

Vr=Vin−ωR

Vr,x=Vx−ωRx=Vx+ωR*sin(θ)

Vr,y=Vy−ωRy=Vy−ωR*cos(θ)

B(θ)=atan(Vr,y/Vr,x)

γ(θ)=β(θ)−α(θ),

where atan( ) is the inverse tangent function (known also as an arcus tangent function), α(θ) is an angle of attack, β(θ) is a relative velocity angle, γ(θ) is the pitch angle of a foil, Vr is a relative velocity, Vin is a spatially variable inflow velocity of water, i.e. the wake field W, ωR is a foil velocity, ω is an angular speed of a foil wheel and R is a radius of a foil wheel. The inflow velocity relative to the foil wheel Vin and the rotational velocity of the foil around the foils wheel Loll comprise the relative velocity Vr=Vin−ωR towards the foil. Note that both Vin and hence Vr may be assumed to vary with the rotation angle of the foil wheel θ due to external wake field W from the marine vehicle and/or the wake field W induced by the foils themselves. This relative velocity forms the angle β(θ) with the x-axis. While the pitch angle of a foil is γ(θ), the angle of flow relative to the foil, i.e., the angle of attack becomes α(θ)=β(θ)−γ(θ) which is equal to the equation γ(θ)=β(θ)−α(θ).

Writing out the relative velocity flow angle β(θ) and the (symmetric) target angle of attack α(θ), the pitch angle γ(θ) of a foil 108, 108′ may be written as

γ
⁡
(
θ
)

=

a
⁢

tan
⁢

(

V
y
+

-

cos
⁢

θ

V
x
+

+

sin
⁢

θ

)

-

α
⁡
(
θ
)

=

tan

-
1

(

V
y

-

ω
⁢
R
⁢

cos
⁢

θ

V
x

+

ω
⁢
R
⁢

sin
⁢

θ

)

-

A
*
cos
⁢

(

θ
-

1
2

⁢

sin
⁢

(

2
⁢
θ

)

)

(

2
⁢
A

)

or more conveniently by scaling the relative velocity components by the foil wheel rotational speed ωR (scaling is optional):

γ
⁡
(
θ
)

=

a
⁢

tan
⁢

(

V
y
+

-

cos
⁢

θ

V
x
+

+

sin
⁢

θ

)

-

A
*
cos
⁢

(

θ
-

1
2

⁢

sin
⁢

(

2
⁢
θ

)

)

,

(

2
⁢
B

)

where V_x⁺ and V_y⁺ are the scaled contributions or components of the inflow velocity Vin of the wake field Win the directions of the coordinate system x and y. The wake field W contributions may be obtained via particle image velocimetry (PIV) from the actual device, for instance.

Although the measurement with at least on sensor 122 may be possible, the fact that not only is the wake field W required as a function θ but also exactly at the same instant an angle θ of the foil wheel 106, 106′ is required, may make it challenging. In an embodiment, the wake field W contributions may be recorded from a CFD simulation. This may be accomplished by placing a computational probe 500 adjacent to a foil 108, 108′ such that it follows the rotation of the foil wheel 106, 106′. In an embodiment, the computational probe does not the foil pitching and may remain static relative to a pivot point of the foil it follows. In this embodiment, the vector pointing from the pivot point to the probe is always parallel to the x-axis. An example of positioning the computational probes 500 at different foil positions is shown in FIG. 5.

A distance DD between a probe 500 and a corresponding foil 108 adjacent to it may be such that the probe 500 is not be too close to the foil 108 such that the foil itself disturbs the flow adjacent to it too much and not too far such that the correlation between the velocity of the water movement at the probe 500 and at the foil adjacent to it remains high. In an embodiment, a suitable distance DD may be about ¾ chord length in front of the pivot point. More generally, the range of the distance DD may be from about 0.5 to about 1 chord length, for example. In principle, the wake field W should be representative over the foil's span (average in z-direction) but for simplicity it may be assumed that a single point at about the middle is good enough. In the case the wake field W is measured, the at least one sensor 122 may be placed in a corresponding location as the computational probe 500.

The wake field W contributions may have quite complex shapes.

In an embodiment, contributions of the wake field W may be continually iterated within the CFD simulation. The control system, particularly a PID (Proportional Integral Derivative) control system requires a continuous analytical function with at least one continuous derivative. If the wake field W is collected at Δθ=1° intervals, i.e. at every degree of a 360° rotation of the foil wheel 106, 106′, from a CFD simulation or from measurements performed by the at least one sensor 122, a discrete Fourier transform may be utilized to form a good model based on continuous trigonometric functions with a total of 360 (full rotation) terms for each wake contribution.

However, feeding such an amount of data on the wake field W to a real machine's control system may be impractical. Furthermore, the real wake field W may not be similar to the simulation in all details. Hence in an embodiment, only a constant average term and/or a few lowest frequency sine and cosine terms from a Fourier series representation of the wake field contributions may be retained. For example three lowest frequency sine and cosine terms may be enough. These truncated functions may then be processed according to Equation (2A) or (2B) in the controller 112, for example. Based on this kind of an algorithm, the controller 112 may form the data on the pitch angle of a foil in a new CFD simulation. Since the wake field W also depends on the pitch angle of a foil, this may be done iteratively in an embodiment. The iterative process may be performed as follows:

- 1) Initialize the pitch angle of a foil in Equation (2B) with

V
x
+

=

-

V

x
,

foil
⁢

wheel

ω
⁢
R

⁢

and
⁢

⁢

V
y
+

=
0

- 2) Obtain converged CFD solution and record actual V_x⁺ and V_y⁺
- 3) Form Fourier series representations by retaining the constant and only a few sine and cosine terms
- 4) Feed a truncated Fourier series representations V_x,trunc⁺ and V_y,trunc⁺ to Equation (1)
- 5) Obtain a new converged CFD solution and record new actual V_x⁺ and V_y⁺
- 6) Form a new truncated Fourier series representations V_x,trunc⁺ and V_y,trunc⁺
- 7) Under-relax to avoid overshoot by averaging with the previous round V_x,trunc⁺ and V_y,trunc⁺
- 8) Return to step 4).

In step 7), too large a change i.e. overshoot should be avoided because a new iterated pitch function also changes the measured V_x⁺ and V_y⁺ and large overshoot will slow down or even prevent the convergence of an iterative process. In an embodiment, only two iterative rounds of the steps 1) to 8) may be enough, because beyond that the changes typically become marginal. A comparison of the contributions of a truncated wake field W to actual ones are presented in FIG. 6 and its explanation.

In an embodiment, while the angle of attack α(θ) is kept constant for a maximized length round the rotation of the foil wheel 106, 106′, and the data on the pitch angle γ(θ) is iterated to a single design point (velocity and RPM), the contributions V_x⁺ and V_y⁺ of the wake field W may be fixed in Equation (1), which determines the pitch angle γ(θ) of a foil 108, 108′. Then, the pitch trajectory of a foil is fixed regardless of the speed and RPM (Rotation Per Minute), and the foil wheel 106, 106′ acts similarly to a fixed pitch screw propeller (or a trochoidal foil wheel with constant eccentricity). If it is assumed that the RPM is constant, at slow speeds the actual angle of attack α(θ) is large, and a high thrust with reduced efficiency is obtained. As the speed increases towards the designed speed, the actual angle of attack α(θ) approaches the target behaviour, and a high efficiency is recovered.

However, this behaviour may also be improved in an embodiment. Decompose now a scaled wake field V_xy⁺ into a sum

V_xy⁺=V_xy,0⁺+V_xy,ind⁺ (3)

where the scaled wake field V_xy⁺ is expressed as a sum of an undisturbed contribution V_xy,0⁺ without a foil wheel 106, 106′, which is with the marine vehicle 100 in self-propulsion situation, real or simulated, and an extra contribution V_xy,ind⁺ induced by the foils 108, 108′.

The undisturbed contribution V_xy,0⁺ of the wake field is directly available since undisturbed speed and RPM are easily obtained. In an embodiment, it may be assumed that velocities of the induced dimensional contribution V_xy,indof the wake field W are locally (at any given θ) directly proportional to the wheel rotational speed ωR. It may be reasonable to assume that the induced velocity is indeed proportional to the velocity inducing it. With this assumption, the scaled induced contributions V_xy,ind⁺ of the wake field W are locally constant but still variable with respect to the angle θ of a rotation of the foil wheel 106:

V_xy,ind⁺=(V_xy⁺−V_xy,0⁺)_design=locally_constant. (4)

This means that the scaled induced contributions V_xy,ind⁺ of the wake field W are cyclically the same for every rotation of around the rotation axis of the foil wheel 106. Hence, by determining the scaled wake field V_xy⁺ at the design speed and RPM, it is possible to use Equations (3) and (4) to reconstruct a reasonable estimate of the wake field W at any other operating point. Equation (4) is approximate but overall the reconstructed wake field W will give a close match to the target angle of attack α(θ) at an arbitrary operating point. With the reconstructed wake field W, a high efficiency can be obtained over a wide range of operating points.

While the angle of attack α(θ) is constant along a maximized length of the rotation of the foil wheel 106, 106′, the pitch angle of a foil may reduce the side force. In an embodiment, it is possible to modify the target angle of attack α(θ) in order to adjust the side force to zero or a large value. In an embodiment, the side force may be adjusted to large value in order to steer the marine vehicle 100 to turn. From the steering perspective, it may however be desirable that the propulsion system 102 has zero or only a weak side force in a straight ahead condition. Both the instantaneous thrust and the side force depend on the corresponding angle of attack α(θ). Well below the stall angle, a larger angle of attack α(θ) increases both.

For example, with a constant ±15° angle of attack, the net side force is positive and by increasing the leading side (around θ=0/360°) angle of attack or decreasing the trailing side (around θ=180°) angle of attack, the positive net side force will decrease. For this, an asymmetric constant angle of attack, Equation (1), may be inserted to Equation (2A) or (2B) in place of the symmetric constant angle of attack, and setting A=14° and B=1°, for example. An asymmetric smooth square wave can be obtained with approximately constant value +15° on the leading side and −13° on the trailing side (see FIG. 3). To compensate for the loss of thrust due to a smaller angle of attack on the trailing side, the rotation rate may be increased. With these settings, the iterative process to obtain an approximate wake field W is repeated. An example of the resulting actual angle of attack α(θ) with the target angle of attack α(θ)_targetand the instantaneous foil efficiency are presented in FIG. 6. Again, the actual angle of attack α(θ) follows the target angle of attack α(θ)_targetquite well, and high efficiency levels are maintained. The performance results behave as expected. The efficiency is further increased and the side force is reduced compared with the prior art in this optimization example, but the thrust is slightly reduced. The reduced thrust may be compensated with a further increase in rotation rate or by adjusting the asymmetry a bit down, say to A=14.1 and B=0.9°. The central performance figures of this +15°/−13° asymmetric constant angle of attack case are listed in Table 1. For reference, the corresponding performance figures of the symmetric +15°/−15° constant angle of attack case are listed in Table 2.

TABLE 1