*2.2. 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. **Figure 4.** Wingsail geometry and parameterization.

*yu y u y y* ρρτ **Figure 5.** Box domain for the freestream simulation. **Figure 5.** Box domain for the freestream simulation.

τ

+

( ) *wall*

0 *y y* ν μρ μ = <sup>∂</sup> == = <sup>∂</sup> (2) where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10−5c 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. *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): *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<sup>+</sup> is defined by Equation (2):

$$y^{+} = \frac{y\mu\_{\tau}}{\nu} = \frac{\rho y}{\mu} \sqrt{\frac{\tau\_{\text{null}}}{\rho}} = y \sqrt{\frac{\rho}{\mu} (\frac{\partial u}{\partial y})}\_{y=0} \tag{2}$$

(a) (b) **Figure 6.**The encrypted mesh in (**a**) mid-span and (**b**) bottom wallof the slot. where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10−5c 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 as 2.871 <sup>×</sup> <sup>10</sup>−<sup>5</sup> c on the wingsail surface. In this case, the profile of the dimensionless wall distance y<sup>+</sup> on the surface of the wingsail with AOA = 6 ◦ is shown in Figure 7. It can be seen that the value of y<sup>+</sup> 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.

**Figure 6.**The encrypted mesh in (**a**) mid-span and (**b**) bottom wallof the slot.

research of wingsail 3D.

research of wingsail 3D.

Cl [-]

0

= <sup>∂</sup> == = <sup>∂</sup> (2)

*y*

*y*

**Figure 4.** Wingsail geometry and parameterization.

**Figure 5.** Box domain for the freestream simulation.

*yu y u y y*

 μρ

τ

where y is the distance from the first mesh cell to the nearest wall. In this paper, y is specified as2.871×10−5c on the wingsail surface. In this case, the profile of the dimensionless wall distance y+ on

nondimensional wall distance for a wall-bounded flow y+ is defined by Equation (2):

τ ρ

ν

+

*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

( ) *wall*

 ρ

 μ

**Figure 6. Figure 6.** The encrypted mesh in ( The encrypted mesh in ( **a**) mid-span and ( **a b**) bottom wallof the slot. ) mid-span and (**b**) bottom wallof the slot.

*J. Mar. Sci. Eng.* **2019**, *7*, x FOR PEER REVIEW 6 of 16

**Figure 7.**Wall y+ contour on wingsail surface (AOA=6 deg). **Figure 7.** Wall y+ contour on wingsail surface (AOA = 6 deg). **Figure 7.**Wall y+ contour on wingsail surface (AOA=6 deg).

*Mesh Independence.* It is considered a necessary step to analyze the influence of grid characteristics. When *Re* =5×105, 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 *Mesh Independence.* It is considered a necessary step to analyze the influence of grid characteristics. When *Re* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> , 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×105, 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.

(**a**) (**b**) **Figure 8.** Convergence of (**a**) lift and (**b**) drag coefficients as a function of mesh number. **Figure 8.** Convergence of (**a**) lift and (**b**) drag coefficients as a function of mesh number.

(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

with different grid sizes is also analyzed when *Re*=5×105 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

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* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> and AOA <sup>=</sup> <sup>15</sup>◦ . 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. *J. Mar. Sci. Eng.* **2019**, *7*, x FOR PEER REVIEW 7 of 16

#### *2.3. Computational Approach 2.3. Computational Approach*

pressure velocity coupling scheme [24].

characteristics of naca0018 wing.

*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×105), the fluid is considered as incompressible (the Mach number of the inlet is about 0.062). *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* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> ), the fluid is considered as incompressible (the Mach number of the inlet is about 0.062). These models are

leading-edge curve of the flap is taken as the length interval (i.e., △x=5×10−4m). Therefore, the magnitude of the time step △t is set at a value of 2.5×10−5s in this simulation investigation. The turbulence intensity at inlet is 1%. The discrete format is quick. Simple algorithm is adopted for the

*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 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* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> (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., <sup>∆</sup><sup>x</sup> <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> m). Therefore, the magnitude of the time step <sup>∆</sup>t is set at a value of 2.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> 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. *J. Mar. Sci. Eng.* **2019**, *7*, x FOR PEER REVIEW 8 of 16

**Figure 10.** Comparison of lift coefficient (**a**) and drag coefficient(**b**) between test and CFD. **Figure 10.** Comparison of lift coefficient (**a**) and drag coefficient (**b**) between test and CFD.

In the case of *Re*=5×105, 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. In the case of *Re* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> , 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.

#### **3.Numerical Results 3. Numerical Results**

### *3.1. Two-Dimensional Wingsail Configurations Study 3.1. 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: The lift force and drag force are converted into dimensionless lift coefficient *C<sup>L</sup>* and drag coefficient *CD*, respectively. They are defined as follows:

$$C\_L = \frac{L}{0.5\rho v^2 A\_R} \tag{3}$$

$$C\_D = \frac{D}{0.5\rho v^2 A\_R} \tag{4}$$

where *L* and *D* are the lift force and drag force of the wingsail, respectively.*A*R is the aspect area of wingsail, include the wing and flap. The lift coefficient and drag coefficient characteristics of where *L* and *D* are the lift force and drag force of the wingsail, respectively. *A<sup>R</sup>* 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.

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. *Camber and the rotation axis of the flap e*ff*ect:* as the key parameter, the rotating axis position of the flap in the direction of the wing chord X<sup>r</sup> is selected as 75% and 85%; five cambers (*d* = 5 ◦ , 10◦ , 15◦ , 20◦ , 25◦ ) have been selected as Figure 11.

the flap, flap deflection angle and slot width when choosing flap geometry parameters.

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

When Xr is 85%, the stall does not occur advance with the increase of flap deflection angle at

1.6

**Figure 11.**Convergence of lift coefficients as a function of AOA and flap deflection angle with (**a**) Xr=85% and (b) Xr=75%. **Figure 11.** Convergence of lift coefficients as a function of AOA and flap deflection angle with (**a**) Xr = 85% and (**b**) Xr = 75%. 2 2

1.6

**Figure 12.** The slot width with Xr=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 (e2/c2) has little effect on the lift and drag coefficients of *d*= 5° and 15° (low camber), but When X<sup>r</sup> is 85%, the stall does not occur advance with the increase of flap deflection angle at low camber as X<sup>r</sup> is 75%. It can be illustrated that the slot width with X<sup>r</sup> = 75% exceeds that with X<sup>r</sup> = 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. (**a**) Xr=85% (**b**) Xr=75% **Figure 11.**Convergence of lift coefficients as a function of AOA and flap deflection angle with (**a**) Xr=85% and (b) Xr=75%. AOA [deg] Cl [-] 0 4 8 12 16 0.4 0.8 1.2 **d=5o Xr =85% d=10o Xr =85% d=15o Xr =85% d=20o Xr =85% d=25o Xr =85%** AOA [deg] Cl [-] 0 4 8 12 16 0.4 0.8 1.2 **d=5o Xr =75% d=10o Xr =75% d=15o Xr =75% d=20o Xr =75% d=25o Xr =75%**

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, **Figure 12.** The slot width with Xr=85% and 75%. **Figure 12.** The slot width with Xr = 85% and 75%.

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. **Table 2.**Results for α=6° with the SST transitional turbulence model. **Configurations e2/c2 d CL CD L/D**  *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 *Camber and flap thickness e*ff*ect:* 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.

r1.51815x0.85g2.4α6d5 15% 5° 1.005 0.0213 47.16 r1.51815x0.85g2.4α6d15 15% 15° 1.724 0.0338 50.94 r1.51815x0.85g2.4α6d25 15% 25° 1.951 0.06 32.53 r1.51812x0.85g2.4α6d5 12% 5° 1.002 0.0207 48.5 r1.51812x0.85g2.4α6d15 12% 15° 1.73 0.0336 51.47 r1.51812x0.85g2.4α6d25 12% 25° 1.606 0.1628 9.87 r1.51810x0.85g2.4α6d5 10% 5° 1.002 0.0208 48.22 r1.51810x0.85g2.4α6d15 10% 15° 1.74 0.0337 51.59 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. 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.

**Table 2.**Results for α=6° with the SST transitional turbulence model.

r1.51810x0.85g2.4α6d25 10% 25° 1.649 0.1649 10

**Configurations e2/c2 d CL CD L/D** 

r1.51815x0.85g2.4α6d15 15% 15° 1.724 0.0338 50.94 r1.51815x0.85g2.4α6d25 15% 25° 1.951 0.06 32.53 r1.51812x0.85g2.4α6d5 12% 5° 1.002 0.0207 48.5 r1.51812x0.85g2.4α6d15 12% 15° 1.73 0.0336 51.47 r1.51812x0.85g2.4α6d25 12% 25° 1.606 0.1628 9.87 r1.51810x0.85g2.4α6d5 10% 5° 1.002 0.0208 48.22 r1.51810x0.85g2.4α6d15 10% 15° 1.74 0.0337 51.59 r1.51810x0.85g2.4α6d25 10% 25° 1.649 0.1649 10


**Table 2.** Results for α = 6 ◦ with the SST transitional turbulence model.

**Figure 13.** The streamlines in the slot region (URANS). **Figure 13.** The streamlines in the slot region (URANS). **Figure 13.** The streamlines in the slot region (URANS).

#### *3.2. Three-Dimensional Study about Effect of Flap Rotation Axis Position 3.2. Three-Dimensional Study about E*ff*ect of Flap Rotation Axis Position 3.2. Three-Dimensional Study about Effect of Flap Rotation Axis Position*

#### 3.2.1. Aerodynamic Performance 3.2.1. Aerodynamic Performance

configuration.

regions over the suction surface of the flap at Xr=90%.

regions over the suction surface of the flap at Xr=90%.

3.2.2. Streamlines

3.2.2. Streamlines

3.2.1. 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 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 theaerodynamic 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<sup>r</sup> <sup>=</sup> 75%,80%, 85%, 90% and 95%) at <sup>α</sup> <sup>=</sup> <sup>6</sup> ◦ and 15◦ . 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°。

(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 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 axismoving backward; especially from 80% to 85%, it drops suddenly at <sup>α</sup> <sup>=</sup> <sup>15</sup>◦ (stall has occurred). It can be illustrated that the change is caused by the flow separation of the wing wake and the flap suctionsurface from Figure 16. At <sup>α</sup> <sup>=</sup> <sup>6</sup> ◦ 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 Figure17b that the flow separation on the suction surface of the flap has disappeared. 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.

**Figure 14.**Lift coefficient versus flap rotation axis position of the wingsail in low and high camber **Xr /c1 [-] Figure 14.**Lift coefficient versus flap rotation axis position of the wingsail in low and high camber configuration. **Figure 14.** Lift coefficient versus flap rotation axis position of the wingsail in low and high camber configuration.

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

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
