Numerical Investigation of a Two-Element Wingsail for Ship Auxiliary Propulsion

: The rigid wingsail is a new type of propulsion equipment which greatly improves the performance of the sailboat under the conditions of upwind and downwind. However, such sail-assisted devices are not common in large ships because the multi-element wingsail is sensitive to changes in upstream ﬂow, making them di ﬃ cult to operate. This problem shows the need for aerodynamic study of wingsails. A model of two-element wingsail is established and simulated by the steady and unsteady RANS approach with the k- ω SST turbulence model and compared with the known experimental data to ensure the accuracy of the numerical simulation. Then, some key design and structural parameters (camber, the rotating axis position of the ﬂap, angle of attack, ﬂap thickness) are used to characterize the aerodynamic characteristics of the wingsail. The results show that the position of the rotating shaft of the ﬂap has little inﬂuence on the lift coe ﬃ cient at low camber. When stall occurs, the lift coe ﬃ cient ﬁrst increases and then decreases as the ﬂap axis moves backward, which also delays the stall angle at a low camber. At the high camber of AOA = 6 ◦ , the lift coe ﬃ cient always increases with the increase of the rotating axis position of the ﬂap; especially between 85% and 95%, the lift coe ﬃ cient increases suddenly, which is caused by the disappearance of large-scale ﬂow separation on the suction surface of the ﬂap. It reﬂects the nonlinear coupling e ﬀ ect between camber of wingsail and the rotating axis position of the ﬂap


Introduction
The rigid wingsail is a common auxiliary propulsion device. Due to its characteristics of environmental protection and energy saving, the rigid wingsail has been applied to various types of ships. In 1980, two rectangular rigid wingsails with a total sail area of 194.4 m 2 on "Shin Aitoku Maru" were installed in Japan, which was the world's first modern sail-assisted commercial tanker. After four years of actual sailing, oil tankers can save 8.5% a year compared with conventional ships [1]. Although the international oil price fluctuates greatly, the research on sail-assisted vessels varies from country to country. In 2018, China Shipbuilding Group delivered the world's first Very Large Crude Carrier (VLCC), "Kai Li", with wingsails ( Figure 1). The trial results show that the energy-saving efficiency of the VLCC is obvious [2]. At the same time, the wingsails have also been rapidly developed in the field of unpowered navigation. Particularly in 2010, BMW Oracle was powered by a 60-meterhigh, multi-element wingsail, which won the 33rd America's Cup. The flexible size and excellent performance of the wingsail has rekindled the interest of the shipping industry [3]. However, due to the phenomenon of flow separation or stall on the surface of the wingsail, the propulsion performance of the wingsail will be deteriorated, which will seriously affect the stability rapidly developed in the field of unpowered navigation. Particularly in 2010, BMW Oracle was powered by a 60-meterhigh, multi-element wingsail, which won the 33rd America's Cup. The flexible size and excellent performance of the wingsail has rekindled the interest of the shipping industry [3]. However, due to the phenomenon of flow separation or stall on the surface of the wingsail, the propulsion performance of the wingsail will be deteriorated, which will seriously affect the stability of the ship. The New Zealand Chiefs[4] suffered a shipwreck in the 35th America's Cup ( Figure 2). Therefore, it is very important to find an effective method to control flow separation or delay stall.  In order to improve the stall characteristics, different flow control methods are used, which can be divided into active control and passive control. The active control method has been applied to the controllable loop sail [5], fluid injection [6,7], turbine sail [8], Magnus sail [9], trailing flap [10,11], leading slat [11,12], etc. The passive method is also applied to different types of technologies, such as Walker sails [13], deformed flaps [14] and leading-edge tubercles [15], etc. powered by a 60-meterhigh, multi-element wingsail, which won the 33rd America's Cup. The flexible size and excellent performance of the wingsail has rekindled the interest of the shipping industry [3]. However, due to the phenomenon of flow separation or stall on the surface of the wingsail, the propulsion performance of the wingsail will be deteriorated, which will seriously affect the stability of the ship. The New Zealand Chiefs[4] suffered a shipwreck in the 35th America's Cup ( Figure 2). Therefore, it is very important to find an effective method to control flow separation or delay stall.  In order to improve the stall characteristics, different flow control methods are used, which can be divided into active control and passive control. The active control method has been applied to the controllable loop sail [5], fluid injection [6,7], turbine sail [8], Magnus sail [9], trailing flap [10,11], leading slat [11,12], etc. The passive method is also applied to different types of technologies, such as Walker sails [13], deformed flaps [14] and leading-edge tubercles [15], etc. In order to improve the stall characteristics, different flow control methods are used, which can be divided into active control and passive control. The active control method has been applied to the controllable loop sail [5], fluid injection [6,7], turbine sail [8], Magnus sail [9], trailing flap [10,11], leading slat [11,12], etc. The passive method is also applied to different types of technologies, such as Walker sails [13], deformed flaps [14] and leading-edge tubercles [15], etc.
In recent years, the flap setting of the wingsail is considered to be a very feasible active control method to control the flow separation. This idea has been applied in the American Cup Sailing Competition with great success. The rigid wingsail consists of two or three symmetrical wings. There is a gap between them to control the wingsail camber on starboard tack and port tack (Figure 3), so as to improve the propulsion performance and delay stall. Many scholars have also carried out research on the aerodynamic characteristics.
The aim of this paper is to understand the influence of geometric parameters on aerodynamic performance. The parametric design of the wingsail is selected and described in Figure 4. In order to In 1996, Daniel [11] designed a high-performance, three-element wingsail. The experimental results show that the maximum thrust coefficient of the three-element wingsail is increased by 68%, the stall angle is delayed between 4 • and 6 • , and the thrust in the whole operation area is also improved. At the most efficient wind direction angle of the wingsail, the thrust is increased by 83%. This proves the superiority of multi-element wingsail in propulsion characteristics. In 2003, Fujiwara [16,17] switched from a triangular sail to a rectangular one. Its propulsion performance is superior to that of Nojiri's hybrid dynamic sail. In 2015, he cooperated with Qiao Li [14] to modify the hybrid sail, replacing the flap with a rigid plate and controlling the rotation of the wing and the plate. Therefore, the sail was also called a variable camber sail(VCS). Through simulation and test, the aerodynamic performance of the VCS was better than that of naca0021 sail and the plate sail. Vincent [18] carried out the simulating research of the two-element wingsail with some key design parameters (camber, slot width, angle of attack, flap thickness, etc.) The results show that the transition phenomenon of flap boundary layer occurs due to the existence of laminar separated bubbles and the interaction between the wing and flap boundary layer in the slot region. The stall is related to slot leakage flow with nonlinear coupling and the leakage flow is affected by flap deflection, slot width and flap thickness. Through PIV measurement and numerical simulation, Alessandro [19] explained the flow phenomena of rigid wingsails at low and high cambers in detail, especially the flow behavior near the slot between two element wings. He pointed out that the flow field in the gap is the key factor affecting the performance of the wingsail. However, the internal relationship between the change of slot size and the development of wingsail aerodynamic characteristics is still unclear. In order to better understand the aerodynamic characteristics and better control of the rigid wingsail, some numerical studies are proposed in this paper.

The Geometry of the Wingsail
The geometry is a simplified configuration of a two-element wingsail with 1/20th model-scale (as Figure 3 and Table 1). The nondimensional number for characterizing the flow of fluid Reynolds number is defined by Equation (1): where µ is the viscosity coefficient of air. The results of naca0018 airfoil with Re = 5 × 10 5 (the wind speed is assumed to be 20 m/s) are compared with the existing experiments to ensure the accuracy of the numerical results. The general view of the geometry and its characteristic parameters are as follows: The aim of this paper is to understand the influence of geometric parameters on aerodynamic performance. The parametric design of the wingsail is selected and described in Figure 4. In order to simplify the general problem, it has been decided to focus on the flap deflection angle d, the position of flap rotation axis in the direction of the wing chord X r , and flap thickness e 2 /c 2 . The angle of attack (α) of the wing represents the angle of attack (AOA) of the wingsail, as seen in Figure 4.
Based on the preliminary simulation and previous research [11,18], the initial wingsail configuration has chord ratio c 1 /c 2 = 3:2, wing thickness e 1 /c 1 = 18%, flap thickness e 2 /c 2 = 15%, one flap rotation axis X r /c 1 = 85%, slot width g/c 1 = 2.4%. Its name will be r1.5t1815 X85g2.4. Adding the angle of attack (α) and flap deflection angle (d) after this name, and the full configuration name of the wingsail will be r1.5t1815 X85g2.4α6 d15, for α = 6 • , d = 15 • .  Mesh Details. The quality of grid determines the accuracy of numerical results. Particular attention has been paid to refine the mesh in the slot, being a region of interaction between the wake of the wing, flow from the slot and the flap boundary layer. Mesh resolution close to the slot has been encrypted as shown in Figure 6. In addition, the mesh size of the flap leading-edge curve is set as 0.4545% (in the case of mesh number 9.86M). In order to ensure the accurate simulation of the flow field near the wall, the height of the first grid cell closest to the wall should be suitable. The nondimensional wall distance for a wall-bounded flow y + is defined by Equation (2): where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10 −5 c on the wingsail surface. In this case, the profile of the dimensionless wall distance y + on the surface of the wingsail with AOA =6° is shown in Figure 7. It can be seen that the value of y + on the sail surface is about 1, ensuring the accuracy of numerical calculation results. The growth ratio of the prism layer is 1.2.The total number of nodes in the grid is about 9.86 million.

The Grid Structure
Computational Domain. Before mesh generation, the calculation domain of geometric model should be determined. The computational domain of wingsail is large enough (32c × 30c × 10c) to consider the side effect as negligible ( Figure 5). AOA is adjusted by rotating the direction of the airflow. The chord of the wingsail is 0.35 m and the span 0.7 m. In order to simplify the model, the atmospheric boundary layer is not considered. The boundary conditions used are velocity inlets on three sides (inlet, starboard boundary, port boundary) and a static pressure outlet equal to the far-field pressure. The velocity at the inlet is uniform, and the value is the same as the free flow velocity. The wingsail surface and the bottom surface of the computational domain are defined as anti-slip walls.   Mesh Details. The quality of grid determines the accuracy of numerical results. Particular attention has been paid to refine the mesh in the slot, being a region of interaction between the wake of the wing, flow from the slot and the flap boundary layer. Mesh resolution close to the slot has been encrypted as shown in Figure 6. In addition, the mesh size of the flap leading-edge curve is set as 0.4545% (in the case of mesh number 9.86M). In order to ensure the accurate simulation of the flow field near the wall, the height of the first grid cell closest to the wall should be suitable. The nondimensional wall distance for a wall-bounded flow y + is defined by Equation (2): where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10 −5 c on the wingsail surface. In this case, the profile of the dimensionless wall distance y + on the surface of the wingsail with AOA =6° is shown in Figure 7. It can be seen that the value of y + on the sail surface is about 1, ensuring the accuracy of numerical calculation results. The growth ratio of the prism layer is 1.2.The total number of nodes in the grid is about 9.86 million.
(a) (b)  Mesh Details. The quality of grid determines the accuracy of numerical results. Particular attention has been paid to refine the mesh in the slot, being a region of interaction between the wake of the wing, flow from the slot and the flap boundary layer. Mesh resolution close to the slot has been encrypted as shown in Figure 6. In addition, the mesh size of the flap leading-edge curve is set as 0.4545% (in the case of mesh number 9.86M). In order to ensure the accurate simulation of the flow field near the wall, the height of the first grid cell closest to the wall should be suitable. The nondimensional wall distance for a wall-bounded flow y + is defined by Equation (2): where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as 2.871 × 10 −5 c on the wingsail surface. In this case, the profile of the dimensionless wall distance y + on the surface of the wingsail with AOA = 6 • is shown in Figure 7. It can be seen that the value of y + on the sail surface is about 1, ensuring the accuracy of numerical calculation results. The growth ratio of the prism layer is 1.2. The total number of nodes in the grid is about 9.86 million. where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10 −5 c on the wingsail surface. In this case, the profile of the dimensionless wall distance y + on the surface of the wingsail with AOA =6° is shown in Figure 7. It can be seen that the value of y + on the sail surface is about 1, ensuring the accuracy of numerical calculation results. The growth ratio of the prism layer is 1.2.The total number of nodes in the grid is about 9.86 million.
(a) (b)  Mesh Independence. It is considered a necessary step to analyze the influence of grid characteristics. When Re =5×10 5 , four different grid numbers (including 2.57, 4.22, 9.86 and 14.4M) are used to estimate the grid sensitivity. In the case of AOA = 6°, the grid independence analysis is carried out with steady Reynolds Average N-S equation and AOA=15° with URANS. Figure 8 reports the lift and drag coefficients of the unmodified sail model. As shown in Figure 8, the number of grids has a great influence on the lift coefficient, and the error of lift coefficient is less than 0.3% between 9.86 and 14.4M. Based on these calculations, the increase in the number of grids (relative to the case of 9.86M) results in a small change in lift. And the drag coefficient, whose change can be considered as acceptable. Therefore, the grid number of 9.86M is used in this study. In order to further verify the reliability of the grid with URANS, the influence on the flow field with different grid sizes is also analyzed when Re=5×10 5 and AOA =15°. For four different grids (including 2.57, 4.22, 9.86 and 14.4M), the velocity vector on the x-y plane is shown in Figure 9. It can Mesh Independence. It is considered a necessary step to analyze the influence of grid characteristics. When Re = 5 × 10 5 , four different grid numbers (including 2.57, 4.22, 9.86 and 14.4M) are used to estimate the grid sensitivity. In the case of AOA = 6 • , the grid independence analysis is carried out with steady Reynolds Average N-S equation and AOA = 15 • with URANS. Figure 8 reports the lift and drag coefficients of the unmodified sail model. As shown in Figure 8, the number of grids has a great influence on the lift coefficient, and the error of lift coefficient is less than 0.3% between 9.86 and 14.4M. Based on these calculations, the increase in the number of grids (relative to the case of 9.86M) results in a small change in lift. And the drag coefficient, whose change can be considered as acceptable. Therefore, the grid number of 9.86M is used in this study. Mesh Independence. It is considered a necessary step to analyze the influence of grid characteristics. When Re =5×10 5 , four different grid numbers (including 2.57, 4.22, 9.86 and 14.4M) are used to estimate the grid sensitivity. In the case of AOA = 6°, the grid independence analysis is carried out with steady Reynolds Average N-S equation and AOA=15° with URANS. Figure 8 reports the lift and drag coefficients of the unmodified sail model. As shown in Figure 8, the number of grids has a great influence on the lift coefficient, and the error of lift coefficient is less than 0.3% between 9.86 and 14.4M. Based on these calculations, the increase in the number of grids (relative to the case of 9.86M) results in a small change in lift. And the drag coefficient, whose change can be considered as acceptable. Therefore, the grid number of 9.86M is used in this study. In order to further verify the reliability of the grid with URANS, the influence on the flow field with different grid sizes is also analyzed when Re=5×10 5 and AOA =15°. For four different grids (including 2.57, 4.22, 9.86 and 14.4M), the velocity vector on the x-y plane is shown in Figure 9. It can be observed that there is a large separation vortex on the suction surface of the flap and a small  In order to further verify the reliability of the grid with URANS, the influence on the flow field with different grid sizes is also analyzed when Re = 5 × 10 5 and AOA = 15 • . For four different grids (including 2.57, 4. 22, 9.86 and 14.4M), the velocity vector on the x-y plane is shown in Figure 9. It can be observed that there is a large separation vortex on the suction surface of the flap and a small vortex on the wake of the wing, which is generated by grids of 9.86 and 14.4M, respectively, but there is no obvious flow separation for grids of 2.57 and 4.22M. There was no significant difference in the flow profile between 9.86 and 14.4M. Therefore, the grid number of 9.86M is suitable for the research of wingsail 3D.

Computational Approach
Numericalmethod. Menter [20] established the k-ω turbulence model of shear stress transport (SST). Hassan [21] compared the numerical predicted results of standard, RNG and realizable K-ε, standard and SST K-ω models with the experimental data, and found that the K -ω SST model is the most accurate prediction model. So the steady Reynolds Averaged Navier-Stokes(RANS) equations are solved thanks to ANSYS Fluent finite-volume solver with low angle of attack before stall occurs; then, the unsteady Reynolds Average N-S equation is used based on the SST K-ω models when physical instabilities exist [22,23]. As the low Reynolds number is based on the chord of the wingsail (Re=5×10 5 ), the fluid is considered as incompressible (the Mach number of the inlet is about 0.062). These models are assigned solvers which control the number of the time step. The selected time step size is computed in the case of CFL=V△t/△x=1. The impact of time models on the aerodynamic performance of the wingsail is investigated at Re=5×10 5 (i.e., V=20.87 m/s), and the mesh size of the leading-edge curve of the flap is taken as the length interval (i.e., △x=5×10 −4 m). Therefore, the magnitude of the time step △t is set at a value of 2.5×10 −5 s in this simulation investigation. The

Computational Approach
Numericalmethod. Menter [20] established the k-ω turbulence model of shear stress transport (SST). Hassan [21] compared the numerical predicted results of standard, RNG and realizable K-ε, standard and SST K-ω models with the experimental data, and found that the K-ω SST model is the most accurate prediction model. So the steady Reynolds Averaged Navier-Stokes(RANS) equations are solved thanks to ANSYS Fluent finite-volume solver with low angle of attack before stall occurs; then, the unsteady Reynolds Average N-S equation is used based on the SST K-ω models when physical instabilities exist [22,23]. As the low Reynolds number is based on the chord of the wingsail (Re = 5 × 10 5 ), the fluid is considered as incompressible (the Mach number of the inlet is about 0.062). These models are 7 of 15 assigned solvers which control the number of the time step. The selected time step size is computed in the case of CFL = V∆t/∆x = 1. The impact of time models on the aerodynamic performance of the wingsail is investigated at Re = 5 × 10 5 (i.e., V = 20.87 m/s), and the mesh size of the leading-edge curve of the flap is taken as the length interval (i.e., ∆x = 5 × 10 −4 m). Therefore, the magnitude of the time step ∆t is set at a value of 2.5 × 10 −5 s in this simulation investigation. The turbulence intensity at inlet is 1%. The discrete format is quick. Simple algorithm is adopted for the pressure velocity coupling scheme [24].
Model validation. In order to verify the numerical results, the lift coefficient of the two-dimensional wingsail model in a free flow environment is compared with the experimental results, as shown in Figure 10. Daniel [11] had carried out experimental research on the lift/drag characteristics of naca0018 wing. In the case of Re=5×10 5 , the RANS method is verified. Before stall (AOA < 15°), the numerical calculation results of the lift coefficient and drag coefficient are close to the experimental values, and the estimation error is less than 5%.At the same time, the prediction of the position of the stall angle is accurate. Although the results of numerical prediction are different from the experimental values after stall, it is acceptable that the 2D numerical calculations used for the qualitative study of the lift/drag characteristics of the wingsail.

Two-Dimensional Wingsail Configurations Study
The lift force and drag force are converted into dimensionless lift coefficient CL and drag coefficient CD, respectively. They are defined as follows: where L and D are the lift force and drag force of the wingsail, respectively.AR is the aspect area of wingsail, include the wing and flap. The lift coefficient and drag coefficient characteristics of wingsail with different geometric parameters are compared.
Camber and the rotation axis of the flap effect: as the key parameter, the rotating axis position of the flap in the direction of the wing chord Xr is selected as 75% and 85%; five cambers (d=5°,10°,15°,20°,25°) have been selected as Figure 11.
When Xr is 85%, the stall does not occur advance with the increase of flap deflection angle at low camber as Xr is 75%. It can be illustrated that the slot width with Xr=75%exceeds that with Xr=85%due to the coupling between the camber of the wingsail and the slot width (as Figure 12). The fluid through the slot is enough to supplement the loss of the wing wake, then delay the stall. When the flap deflection angle is more than 15°, the stall angle begins to reduce, which is caused by flow  In the case of Re = 5 × 10 5 , the RANS method is verified. Before stall (AOA < 15 • ), the numerical calculation results of the lift coefficient and drag coefficient are close to the experimental values, and the estimation error is less than 5%. At the same time, the prediction of the position of the stall angle is accurate. Although the results of numerical prediction are different from the experimental values after stall, it is acceptable that the 2D numerical calculations used for the qualitative study of the lift/drag characteristics of the wingsail.

Two-Dimensional Wingsail Configurations Study
The lift force and drag force are converted into dimensionless lift coefficient C L and drag coefficient C D , respectively. They are defined as follows:   Camber and flap thickness effect: Nine selected simulations are summarized in Table 2 to illustrate the coupling between the camber (flap deflection) and the thickness of the flap. It can be seen that the flap thickness (e2/c2) has little effect on the lift and drag coefficients of d= 5° and 15° (low camber), but has significant effect on d= 25° (high camber) which shows the nonlinear coupling between flap thickness and flap deflection angle.
For high camber, the thickening flap leads to an increase of the leading-edge radius which decreases the pressure coefficient suction peak and has postponed the stall of the wingsail. It emphasizes the complexity of the slot flow of a two-element wingsail through the strong coupling between flap deflection, thickness and slot width. The slot effect may be useful to recall some famous papers with Gentry [25] and Smith [26] which illustrate the complex effect among the wing wake, the leakage flow from the slot and the flap boundary layer as may be seen in Figure 13. Therefore, the slot of the wingsail includes a strong design constraint which include the rotating axis position of the flap, flap deflection, flap thickness and slot width.  Figure 11. Convergence of lift coefficients as a function of AOA and flap deflection angle with (a) X r = 85% and (b) X r = 75%.
When X r is 85%, the stall does not occur advance with the increase of flap deflection angle at low camber as X r is 75%. It can be illustrated that the slot width with X r = 75% exceeds that with X r = 85% due to the coupling between the camber of the wingsail and the slot width (as Figure 12). The fluid through the slot is enough to supplement the loss of the wing wake, then delay the stall. When the flap deflection angle is more than 15 • , the stall angle begins to reduce, which is caused by flow separation of the wing wake, then the flap stalls. The maximum lift coefficient of the wingsail is 2.29 when the flap deflection angle is 20 • with AOA = 12 • . However, due to the reduction of stall angle, the comprehensive performance of the wingsail with d = 20 • is not as good as that with d = 15 • . When the flap deflection angle adds to 25 • , the lift coefficient decreases in the range of full angles of attack, which may be caused by the stalls. Therefore, it is necessary to consider the rotating axis position of the flap, flap deflection angle and slot width when choosing flap geometry parameters.  Camber and flap thickness effect: Nine selected simulations are summarized in Table 2 to illustrate the coupling between the camber (flap deflection) and the thickness of the flap. It can be seen that the flap thickness (e2/c2) has little effect on the lift and drag coefficients of d= 5° and 15° (low camber), but has significant effect on d= 25° (high camber) which shows the nonlinear coupling between flap thickness and flap deflection angle.
For high camber, the thickening flap leads to an increase of the leading-edge radius which decreases the pressure coefficient suction peak and has postponed the stall of the wingsail. It emphasizes the complexity of the slot flow of a two-element wingsail through the strong coupling between flap deflection, thickness and slot width. The slot effect may be useful to recall some famous papers with Gentry [25] and Smith [26] which illustrate the complex effect among the wing wake, the leakage flow from the slot and the flap boundary layer as may be seen in Figure 13. Therefore, the slot of the wingsail includes a strong design constraint which include the rotating axis position of the flap, flap deflection, flap thickness and slot width.  Figure 12. The slot width with X r = 85% and 75%.
Camber and flap thickness effect: Nine selected simulations are summarized in Table 2 to illustrate the coupling between the camber (flap deflection) and the thickness of the flap. It can be seen that the flap thickness (e 2 /c 2 ) has little effect on the lift and drag coefficients of d = 5 • and 15 • (low camber), but has significant effect on d = 25 • (high camber) which shows the nonlinear coupling between flap thickness and flap deflection angle.
For high camber, the thickening flap leads to an increase of the leading-edge radius which decreases the pressure coefficient suction peak and has postponed the stall of the wingsail. It emphasizes the complexity of the slot flow of a two-element wingsail through the strong coupling between flap deflection, thickness and slot width. The slot effect may be useful to recall some famous papers with Gentry [25] and Smith [26] which illustrate the complex effect among the wing wake, the leakage flow from the slot and the flap boundary layer as may be seen in Figure 13. Therefore, the slot of the wingsail includes a strong design constraint which include the rotating axis position of the flap, flap deflection, flap thickness and slot width.

Aerodynamic Performance
The three-dimensional flow around the wingsail at low camber (d= 15°) and high camber (d= 25°) is studied in detail. In order to better study the effect of the rotating axis position of the flap on the aerodynamic characteristics and delayed stall of the wingsail, we analyzed the lift characteristics of the wingsail when the rotating axis of the flap is located at different positions of the wing chord (Xr=75%, 80%, 85%, 90% and 95%) atα=6° and 15°。 Figure 14 shows the lift characteristics for different rotating axis positions of the flap. The position of flap rotating axis has little effect on lift coefficient at α=6° (stall has not yet occurred) with low camper. The lift coefficient first increases and then reduces with the position of flap rotation axis moving backward; especially from 80% to 85%, it drops suddenly at α=15°(stall has occurred). It can be illustrated that the change is caused by the flow separation of the wing wake and the flap suction surface from Figure 16. At α=6° with high camber, the lift coefficient increases with the position of the flap rotating axis moving backward, especially from 85% to 95%. It can be explained from Figure  17b that the flow separation on the suction surface of the flap has disappeared.

Aerodynamic Performance
The three-dimensional flow around the wingsail at low camber (d = 15 • ) and high camber (d = 25 • ) is studied in detail. In order to better study the effect of the rotating axis position of the flap on the aerodynamic characteristics and delayed stall of the wingsail, we analyzed the lift characteristics of the wingsail when the rotating axis of the flap is located at different positions of the wing chord (X r = 75%, 80%, 85%, 90% and 95%) at α = 6 • and 15 • . Figure 14 shows the lift characteristics for different rotating axis positions of the flap. The position of flap rotating axis has little effect on lift coefficient at α = 6 • (stall has not yet occurred) with low camper. The lift coefficient first increases and then reduces with the position of flap rotation axis moving backward; especially from 80% to 85%, it drops suddenly at α = 15 • (stall has occurred). It can be illustrated that the change is caused by the flow separation of the wing wake and the flap suction surface from Figure 16. At α = 6 • with high camber, the lift coefficient increases with the position of the flap rotating axis moving backward, especially from 85% to 95%. It can be explained from Figure  17b

Aerodynamic Performance
The three-dimensional flow around the wingsail at low camber (d= 15°) and high camber (d= 25°) is studied in detail. In order to better study the effect of the rotating axis position of the flap on the aerodynamic characteristics and delayed stall of the wingsail, we analyzed the lift characteristics of the wingsail when the rotating axis of the flap is located at different positions of the wing chord (Xr=75%, 80%, 85%, 90% and 95%) atα=6° and 15°。 Figure 14 shows the lift characteristics for different rotating axis positions of the flap. The position of flap rotating axis has little effect on lift coefficient at α=6° (stall has not yet occurred) with low camper. The lift coefficient first increases and then reduces with the position of flap rotation axis moving backward; especially from 80% to 85%, it drops suddenly at α=15°(stall has occurred). It can be illustrated that the change is caused by the flow separation of the wing wake and the flap suction surface from Figure 16. At α=6° with high camber, the lift coefficient increases with the position of the flap rotating axis moving backward, especially from 85% to 95%. It can be explained from Figure  17b that the flow separation on the suction surface of the flap has disappeared.

Streamlines
As for the wingsails with different rotating axis positions of the flap at α = 15 • and d = 15 • , the streamlines on the mid-span of wingsail are depicted in Figure 15. It can be found that a small separation vortex appears in the wake of the wing and a large separation vortex appears in the regions over the suction surface of the flap at X r = 90%.  An obvious separation helix appears on the suction surface of the wing at Xr=95%, which indicates that the vortex has been formed. When the rotating axis position of the flap is moves backward from 90% to 95%, the flow separation on the suction surface of flap disappears. We guess the fluid flowing through the smaller gap does not supplement the wake of the wing, but flows along the flap boundary layer.

Velocity Magnitude Contours
It can be seen from Figure 17 that the suction surface of the wing has a large flow separation at α=15° with low camber. At Xr=85%, the fluid flowing through the slot complements the flow Figure 16. Limiting streamline and static pressure contours on suction surface for the two-element wingsail at (a) X r = 75% (b) X r = 80% (c) X r = 85% (d) X r = 90% and (e) X r = 95% with α = 15 • , d = 15 • .

Velocity Magnitude Contours
It can be seen from Figure 17 that the suction surface of the wing has a large flow separation at α = 15 • with low camber. At X r = 85%, the fluid flowing through the slot complements the flow separation of the wing wake and the flap suction surface. There is no vortex in the mid-span of the wingsail. The slot jet can be observed clearly at X r = 90%. As a result of the reduced slot width due to the rotating axis position of the flap backward, the slot jet only divides the vortex of the wing wake and there is a large-scale flow separation on the suction surface of the flap, which causes deep stall of the flap. This phenomenon has also been observed and described in Biber [27] and Chapin [18]. Because of the smaller slot width at X r = 95%, the slot jet only flows along the boundary layer of the flap, which has little effect on the separation flow of the wing wake. If the slot is further reduced or not set, the flap setting may aggravate the separation of the wing wake, such as the research of the hybrid sail designed by Qiao Li [14]. However, the forward movement of the rotating axis position of the flap is limited by the flap deflection angle. As shown in Figure 18, when the rotating axis position of the flap moves forward from 90% to 85% at α= 6 • , the large-scale flow separation occurs on the suction surface of the flap, as the phenomenon seen by Fiumara [19] in the two-element sail experiment in 2016. separation of the wing wake and the flap suction surface. There is no vortex in the mid-span of the wingsail. The slot jet can be observed clearly at Xr = 90%. As a result of the reduced slot width due to the rotating axis position of the flap backward, the slot jet only divides the vortex of the wing wake and there is a large-scale flow separation on the suction surface of the flap, which causes deep stall of the flap. This phenomenon has also been observed and described in Biber [27] and Chapin [18]. Because of the smaller slot width at Xr = 95%, the slot jet only flows along the boundary layer of the flap, which has little effect on the separation flow of the wing wake. If the slot is further reduced or not set, the flap setting may aggravate the separation of the wing wake, such as the research of the hybrid sail designed by Qiao Li [14]. However, the forward movement of the rotating axis position of the flap is limited by the flap deflection angle. As shown in Figure 18, when the rotating axis position of the flap moves forward from 90% to 85% at α= 6°, the large-scale flow separation occurs on the suction surface of the flap, as the phenomenon seen by Fiumara [19] in the two-element sail experiment in 2016.
(a) X r = 85% (b) X r = 90% (c) X r = 95% separation of the wing wake and the flap suction surface. There is no vortex in the mid-span of the wingsail. The slot jet can be observed clearly at Xr = 90%. As a result of the reduced slot width due to the rotating axis position of the flap backward, the slot jet only divides the vortex of the wing wake and there is a large-scale flow separation on the suction surface of the flap, which causes deep stall of the flap. This phenomenon has also been observed and described in Biber [27] and Chapin [18]. Because of the smaller slot width at Xr = 95%, the slot jet only flows along the boundary layer of the flap, which has little effect on the separation flow of the wing wake. If the slot is further reduced or not set, the flap setting may aggravate the separation of the wing wake, such as the research of the hybrid sail designed by Qiao Li [14]. However, the forward movement of the rotating axis position of the flap is limited by the flap deflection angle. As shown in Figure 18, when the rotating axis position of the flap moves forward from 90% to 85% at α= 6°, the large-scale flow separation occurs on the suction surface of the flap, as the phenomenon seen by Fiumara [19] in the two-element sail experiment in 2016.

Conclusions
By studying the influence of flap geometric parameters on the aerodynamic characteristics of two-element wingsail under steady and unsteady conditions, two-dimensional and three-dimensional related parameters were simulated using the same Reynolds number (Re =5×10 5 ).
The 2D simulation results show that, when the rotating axis position of the flap is located at 85% of the wing chord, the thickening flap leads to an increase of the leading-edge radius which decreases the pressure coefficient suction peak and has postponed the stall of the wingsail for high camber. It reflects the nonlinear coupling effect between wingsail camber and flap thickness. When the rotating axis of the flap is located at the 75% of the wing chord, the stall angle is delayed with the increase of the flap deflection angle at low camber. When selecting the geometric parameters of the flap, factors such as the position of the flap rotation axis, the flap deflection angle, and the flap thickness need to be considered comprehensively.
The 3D simulation mainly studies the influence of the flap rotating axis position and the flap deflection angle on the stall characteristics of the wingsail. When stall has not yet occurred with low camber, the rotating axis position of the flap has little effect on the lift coefficient, while stall has occurred, the lift coefficient increases first and then reduces with the rotating axis position of the flap moving backward. The flow separation of the suction surface of the wing expands from the root to the top and the return area becomes larger, especially from the 80% to the 85%, the lift coefficient drops suddenly. This is caused by the flow separation between the wing wake and the flap suction surface. At high camber with AOA=6°, the lift coefficient always increases with the position of the flap rotation axis, especially from 85% to 95%, the lift coefficient suddenly rises, which is caused by the disappearance of large-scale flow separation of the flap suction surface. Therefore, the slot width is an important factor affecting the flow separation of the wing wake and the suction surface of the flap, where size is affected by the rotating axis position of the flap and the flap deflection angle. When the flap deflection angle is adjusted to obtain a large lift coefficient, the restriction of the rotating axis position of the flap must be considered to ensure a reasonable stall angle range.

Conclusions
By studying the influence of flap geometric parameters on the aerodynamic characteristics of two-element wingsail under steady and unsteady conditions, two-dimensional and three-dimensional related parameters were simulated using the same Reynolds number (Re = 5 × 10 5 ).
The 2D simulation results show that, when the rotating axis position of the flap is located at 85% of the wing chord, the thickening flap leads to an increase of the leading-edge radius which decreases the pressure coefficient suction peak and has postponed the stall of the wingsail for high camber. It reflects the nonlinear coupling effect between wingsail camber and flap thickness. When the rotating axis of the flap is located at the 75% of the wing chord, the stall angle is delayed with the increase of the flap deflection angle at low camber. When selecting the geometric parameters of the flap, factors such as the position of the flap rotation axis, the flap deflection angle, and the flap thickness need to be considered comprehensively.
The 3D simulation mainly studies the influence of the flap rotating axis position and the flap deflection angle on the stall characteristics of the wingsail. When stall has not yet occurred with low camber, the rotating axis position of the flap has little effect on the lift coefficient, while stall has occurred, the lift coefficient increases first and then reduces with the rotating axis position of the flap moving backward. The flow separation of the suction surface of the wing expands from the root to the top and the return area becomes larger, especially from the 80% to the 85%, the lift coefficient drops suddenly. This is caused by the flow separation between the wing wake and the flap suction surface. At high camber with AOA = 6 • , the lift coefficient always increases with the position of the flap rotation axis, especially from 85% to 95%, the lift coefficient suddenly rises, which is caused by the disappearance of large-scale flow separation of the flap suction surface. Therefore, the slot width is an important factor affecting the flow separation of the wing wake and the suction surface of the flap, where size is affected by the rotating axis position of the flap and the flap deflection angle. When the flap deflection angle is adjusted to obtain a large lift coefficient, the restriction of the rotating axis position of the flap must be considered to ensure a reasonable stall angle range.   The position of flap rotation axis in the direction of the wing chord [-]