Propulsive performance of a rigid wingsail with crescent-shaped profiles

.


Introduction
Transportation accounts for a considerable proportion of greenhouse gas (GHG) emissions.According to 2019 data (Pathak et al., 2022), transportation produces 15% of the total GHG emissions worldwide, with shipping representing 9% of transportation emissions.According to 2022 data (UNCTAD, 2022), at least 80% of global trade by volume and more than 70% by value is carried by shipping.For instance, a pure car and truck carrier (PCTC) may consume 30 to 60 tons of fossil fuel per day, depending on how it is operated (Bialystocki and Konovessis, 2016).In the late 2010s, ocean-going ships consumed an average of 330 million tons of fuel per year and emitted approximately 1056 million tons of CO 2 (IMO, 2020).A previous study (Gadonneix et al., 2011) projected that the world's transportation needs could double by 2050.To address this problem, the International Maritime Organization (IMO) has agreed that shipping must become more energy-efficient and reduce GHG emissions to 50% of 2008 levels by 2050 (IMO, 2018).
Wind-assisted ship propulsion for large commercial ships is considered one of the most promising ways to reduce shipping's dependence on fossil fuels.Various innovative sail technologies have been proposed, such as rotor sails, vertical airfoils (also termed Ventifoils or suction wings), kites, wind turbines, and various wingsails (Khan et al., 2021).
Several of these technologies are already used on passenger and merchant vessels, while some are still being subjected to optimization or full-scale testing in research projects (Cairns et al., 2021).Lu and Ringsberg (2020) compared three sail technologies (the Flettner rotor, the DynaRig, and a classical wingsail with airfoil profiles) with respect to the fuel savings achieved for a specific ship sailing on specific voyage routes.Their study showed that wind-assisted ship propulsion technologies reduced fuel consumption by several percentages (e.g., 5.6% ∼ 8.9%).The amount of fuel savings depended on many factors for each of the three technologies.One of the crucial factors was the sail's performance for the wide range of angles of attack (α) related to the ship's heading direction.A sail's performance also depends on the aerodynamic interactions between the sails on the ship if more than one sail is installed.
Several numerical analyses have been conducted to evaluate the performance of wind-assisted ship propulsion systems based on rigid sails.Ouchi et al. (2011) performed full-scale computational fluid dynamics (CFD) simulations to evaluate the propulsive performance of a nine-wingsail system and carried out a case study for evaluation.Viola et al. (2015) developed a numerical optimization procedure for a rigid wingsail using the Reynolds-averaged Navier-Stokes (RANS) equations solver, which yielded an efficient parametric sail aerodynamic analysis method.Persson et al. (2019) presented simplified approaches to modeling wind-assisted ship propulsion systems using a limited number of CFD simulation results to extrapolate propulsive performance under various conditions.Tillig and Ringsberg (2020) presented a novel approach to analyzing aero-and hydrodynamic interaction effects on wind-propelled ships.The low-aspect-ratio wing theory was applied and modified to predict the lift and drag forces of hulls sailing at various drift angles.The sails' aerodynamic interaction effects were captured by numerically solving the Navier-Stokes equations for incompressible creeping flow.Malmek et al. (2020) developed two cost-effective aerodynamic methods to predict the performance of large-scale wingsails.One was based on the lifting line theory of potential flow in combination with precalculated two-dimensional RANS CFD data, and the other is a vortex lattice method (VLM).Zhu (2020) and Blount and Portell (2021) performed detached eddy simulations (DES) to study the performance of wingsails with a NACA 0015 profile under downwind conditions.
The studies mentioned above were mainly based on wingsails with conventional airfoil profiles, such as the NACA series.It is worth noting that the profiles have no cambers to adapt to wind directions in practice.However, including a camber in the profile geometry can substantially increase the lift coefficient.This mechanism is also valid for sails.Atkinson (2019) performed three-dimensional CFD simulations for studying a segmented rigid sail.Nikmanesh (2021) proposed a cambered crescent-shaped profile and predicted its propulsive performance by conducting unsteady RANS (uRANS) simulations based on the Spalart-Allmaras turbulence model.The results showed that this type of profile provides a higher C L , resulting in greater thrust.Chen et al. (2022) studied the aerodynamic characteristics of a set of arc-shaped wingsails by performing two-dimensional simulations.
This study introduced the conceptual design of a telescopic wingsail with a new crescent-shaped profile.Two-and three-dimensional CFD simulations were performed using uRANS with the k-ω shear stress transport (SST) turbulence model to evaluate its performance.One objective of this study was to examine suitable sail configurations with advantageous lift force coefficients (C L ) to maximize the thrust force (F T ) for large sailing merchant ships and their ship operation profiles.Two rigid wingsails with different sectional profiles (NACA 0015 and crescent-shaped profiles) were compared in terms of their propulsive performance.The analysis also explored the flow field properties, including the flow separation points and the tip vortices.Using the shipmodeling platform ShipCLEAN (Tillig et al., 2019), a case study was conducted in which the crescent-shaped wingsail was applied to a ship to evaluate the wingsail's propulsive performance.

Conceptual design
The propulsive performance of a wind-assisted ship propulsion concept based on a telescopic rigid wingsail with a crescent-shaped sectional profile was analyzed in this study.The rig is designed to have a telescopic function so that it can be operated under various conditions.For example, under weak-wind conditions, the wingsails can be fully expanded to capture maximum thrust, while they can be retracted under strong-wind conditions to protect the structures.Some existing conceptual wingsails, such as the Oceanbird (Workinn, 2021), were designed to have a larger section area at the bottom and a smaller section area close to the tip.In the proposed concept, in contrast, because of the telescopic function, the wingsail is designed to be uniformly extruded, which means that the shape of the sectional profile remains the same through the spanwise direction, as shown in Fig. 1.

Wind triangle
For ships using rigid wingsails, the wind load on the wingsails produces thrust that propels the ships.The thrust force is usually transferred from the wingsails to the ship through the mast.The external loads on the wingsail depend on the apparent wind speed, which is the wind that the ship and the sails experience.
Fig. 2 shows the wind triangle.The apparent wind speed (V AW ), i.e., the wind speed relative to the ship, and the apparent wind angle (θ AW ) can be calculated using Equations ( 1) and (2), respectively (Kimball, 2009).
The external loads on the sail include the forces and the moments, as Fig. 2 shows.From an aerodynamics perspective, the component of the total force on the sail that is parallel to V AW is the drag force (F D ), while  H. Zhu et al. that perpendicular to V AW is the lift force (F L ).In practical terms, the component that is parallel to the ship's speed (V S ) is the thrust force (F T ), which can be calculated using Equation (3), while that perpendicular to V S is the side force (F S ).The magnitude of F T represents the propulsive performance of the wind-assisted propulsion system, so one of the most important objectives of a wind-assisted propulsion concept is to generate as large an F T as possible.F S does not contribute to propulsion and causes heeling, rolling, drift, and induced resistance.
In this study, CFD simulations were conducted to predict the aerodynamic forces F L and F D .The nondimensional force coefficients C L and   4) and ( 5), respectively.The moment coefficient is calculated as shown in Equation ( 6). (4) (5) Similarly, the thrust force coefficient can be obtained using Equation (7).

Crescent-shaped sectional profile
For a sail to operate equally well on both starboard and port side tack, it can either be symmetrical with respect to the chord or symmetrical with respect to the normal of the chord.According to a rough estimation based on the thin-foil theory (Houghton and Carruthers, 1982), a sectional profile with a substantial camber is expected to provide a higher C L than an aerodynamically symmetrical profile, e.g., NACA 0015.Although the propulsive performance of the NACA 0015 profile can be improved by introducing a flap at the trailing edge, it would make the telescopic function hard to achieve.
Compared with traditional airfoil profiles, such as NACA 0015, a cambered profile that is symmetric at both edges would be operated differently when changing the tack.The traditional airfoil profile will always have the same leading edge, but the pressure side, i.e., the highpressure side, will become the suction side, i.e., the low-pressure side, when changing from one tack to the other, as illustrated in Fig. 3(a).The profiles will then operate like a modern Bermuda-type sail.The crescentshaped profile will swap the leading and trailing edge when changing tack, but the same side will always be the pressure side, and the opposite side will always be the suction side (see Fig. 3(b)).
The geometry of the crescent-shaped profile used in this study is presented in Fig. 4. The geometry is symmetric with respect to the center line that is normal to the profile chord.The chord length of the profile is 14 m.The maximum thickness of 2 m is at the middle of the chord.Given that the surfaces on the pressure and suction sides are curved, a camber is formed along the chord.
The propulsive performance of the profile varies depending on the chosen design parameters.The optimal parameters depend on both the flow and structural conditions.Based on initial estimates (Nikmanesh, 2021), in which a series of two-dimensional CFD simulations were performed, the chosen profile with the best thrust coefficient was expected to be reasonably close to the final design.

CFD simulation 2.2.1. Physical condition
CFD simulations were performed to predict the external load on the sail.Mesh generators and solvers in the commercial software STAR-CCM+ (Siemens PLM Software, 2021, 2021) were used in the present study.
A global coordinate system was introduced for the simulation cases considered in the study.The X axis was in the direction of the inlet flow.The Y axis pointed from the pressure side to the suction side.The Z axis pointed vertically upward from the bottom to the top.The origin was located at the bottom surface.For the crescent-shaped profile, the origin was at the center of the mean camber line (Fig. 5(a)), whereas for NACA 0015, it was located at the leading edge (Fig. 5(b)).
To simulate the critical condition, i.e., a condition with high wind speed, a uniformly distributed inlet flow velocity of 25 m/s was set for the wind speed.The Reynolds number (Re), calculated based on the chord length of the sail sectional profile, was 2.3 × 10 7 .Table 1 lists the properties of the fluid (for air at 25 • C).

Domain and boundary conditions
Two-dimensional simulations were first performed to identify the critical angle of attack (α c ) and provide a preliminary view of the flow field.The results showed that a peak value of C L was obtained when α was approximately 20 • .Three-dimensional simulations with two sets of boundary conditions were then conducted to study the threedimensional flow characteristics and the influence of tip vortices.
The first set of boundary conditions included periodic boundary conditions applied to the top and bottom boundaries of the computational domain in the spanwise direction.The second set, which represented the real conditions more closely, included symmetric boundary conditions imposed at the bottom boundary (where the sail is attached) and the top surface, as Fig. 6 shows.In both sets of boundary conditions, the side boundaries in the crossflow direction were set to the pressure outlets to make the simulations represent the real conditions more closely.In the two-dimensional simulations, the arrangement of the inlet and outlet boundaries followed Fig. 6(a).A non-slip boundary condition was specified for the wall of the wingsail.The upstream boundary of the domain was assigned as the velocity inlet.A pressure outlet boundary condition with a zero-pressure loss coefficient was imposed on the downstream boundary of the domain.These two sets of boundary conditions and two-dimensional simulations were compared to assess how different three-dimensional physical characteristics influence the predicted wind loads, as discussed in Sections 3.1 and 3.2.
The computational domains were rectangular.Full-scale geometries were simulated in this study.Fig. 7 presents the size and arrangement of the three-dimensional simulations with the freestream tip and symmetric bottom boundary conditions.The total height of the sail, i.e., the spanwise length, was 72 m in the simulation cases conducted with both sets of boundary conditions.For simulations with periodic top and bottom boundary conditions, the spanwise size of the domain was the same as the spanwise length of the sail.For the two-dimensional simulations, the size of the domain followed the size of the bottom boundary shown in Fig. 7.

Turbulence modeling
Turbulence was simulated using the uRANS equations with the k-ω SST model (Menter, 1993).This two-equation turbulence model computes the eddy viscosity turbulence based on the turbulence kinetic energy and the specific turbulence dissipation rate.The convection term of the model is discretized with a second-order upwind scheme.The k-ω SST model interprets the standard k-ω model within the inner region of the boundary layer, and as in the freestream, it switches to the k-ε model to avoid the problem that the k-ω model is sensitive to the inlet freestream turbulence.A blended wall treatment approach was applied to the RANS equations (Wilcox, 1989).This approach has the advantage of treating complex geometries with local flow characteristics.Because the velocity over complex walls varies over a wide range and the geometry of the wingsail profile has a curvature, it is difficult to ensure that y + in all cells adjacent to the walls are either above a high value or below a small value, which is needed to apply a conventional wall treatment model.In contrast, the blended wall treatment is regarded as a function of the local y + .Blended wall laws are employed to model smooth variable changes in the buffer layer between the viscous sublayer and the logarithmic region.The Reichardt law (Reichardt, 1951) is utilized for the momentum equations.Zeng et al. (2023) studied the Re sensitivity of a wingsail based on a camber profile and found that the force coefficients are not sensitive to Re for Re values greater than 1 × 10 5 , which indicates that transition does not affect the wind loads for high-Re conditions.In this study, the gamma transition model (Menter et al., 2015), which solves for turbulence intermittency to predict the onset of the transition in the turbulence boundary layer, was introduced to include the influence of the transition flow.For deep-stall conditions, attached flow only exists in a very small area near the edges, so the influence of the transition under large-α conditions is not studied.For the two-dimensional simulation case with α = 20 • as an example, the turbulence on the suction side develops very close to the leading edge, as shown in Fig. 8, where high-k t Fig. 12. Two-dimensional mesh independence study based on different refinement strategies.
H. Zhu et al. areas can be found near the leading edge.The transition model mainly influences the turbulence on the pressure side.Comparing the time-averaged pressure and wall shear stress shows that simulations without the transition model may underestimate C L by 0.3% and C D by 3.5%.However, because F D is relatively much lower (normally less than 10%) than F L for small-α conditions, the influence of the transition flow on the propulsive performance can be ignored.Furthermore, the transition model substantially reduces the simulation speed because an extra equation needs to be solved, so the transition model was not applied in this study.

Solver and discretization schemes
A finite volume method was utilized to discretize the governing equations.The method employs a segregated flow solver that uses the semi-implicit method for the pressure-linked equations (SIMPLE) algorithm (Patankar, 1980).The flow was assumed to be incompressible in this study because of the low freestream Mach number.
The convection fluxes on the cell faces were discretized using a hybrid second-order upwind and bounded-central scheme.The diffusion fluxes on both the internal and boundary cell faces were discretized using a second-order scheme.The second-order hybrid Gauss-LSQ method was used in the gradient computation.This method involves the reconstruction of the field values in a cell face, such as the secondary gradients of the diffusion fluxes and the pressure gradients, as well as the rate-of-strain tensors used in turbulence models.A second-order implicit method was applied to discretize the time derivative.The time-marching procedure adopts iterations within each time step.

Numerical mesh 2.3.1. Mesh typology
Unstructured meshes with a trimmed cell topology were mainly used for the simulations.Fig. 9 shows the mesh of the three-dimensional simulation with a freestream tip and symmetric bottom, together with typical cell sizes.The cells had a uniform size in each region at each refinement level.The region near the foil and the wake region were both refined, which can be seen in the section plane at Z = 0.5H in Fig. 9(b).The mesh in the wake region was refined by two parameters: the length and the separate angle of the wake refinement.A cylindrical volumetric mesh refinement with a length of 1.1H was introduced to refine the mesh near the foil.The diameter of the cylinder was 20 m.Flow separation points were expected to be distributed around the two edges, so similarly, a more refined set of meshing was applied to the region near the edges of the foil to capture the flow separation phenomena, as Fig. 9(c) shows.The refinement, except for the prism layer, was based on the base size, l base , indicated in Fig. 9.
Prism layers were generated near the wall of the foil to resolve the flow in the boundary layers, as shown in Fig. 9(d).The absolute total thickness of the prism layer was 0.5 m, and the number of prism layers was 55 for all of the simulation cases.Because the crescent-shaped profile has a quite large camber, a strong flow separation is expected.Thus, the y + values of the first layer of cells near the wall should be less   than 1 to achieve a sufficiently detailed and accurate assessment of the boundary layer flow.In the simulation cases, the order of magnitude of y + is approximately 10 − 1 on most areas of the wall.To obtain this low y + value, the near-wall thickness of the prism layer was set to an absolute value of 1 × 10 − 5 m, which did not change during global mesh refinement.
Polyhedral meshes were also used to cross-compare with the trimmed meshes described above.By following the same refinement strategy, which means having consistently refined regions, the polyhedral mesh shown in Fig. 10 was generated.Compared with the trimmed mesh, the cell sizes develop smoothly in the polyhedral mesh.The impact of mesh typology is discussed in Section 2.3.3.
As shown in Fig. 3, when θ AW is greater than 90 • , the wingsail operates with a high α.When θ AW is approximately 180 • , i.e., the downwind condition, F D is mainly used as the thrust for propulsion.Two-dimensional and three-dimensional simulations were performed for deep-stall conditions, such as when α = 90 • .Under deep-stall conditions, von Kármán vortex streets are expected to spread for quite a long distance in downstream areas, so an extra downstream refinement for the downstream areas was introduced when α ≥ 35 • (Fig. 11).

Two-dimensional mesh independence
Two-dimensional mesh independence studies were carried out by considering two aspects: the refinement strategy and the size of the cells.
With respect to the refinement strategy, three factors (the existence of the near-foil refinement, the length, and the separate angle of the wake refinement) were studied to ensure that the mesh quality was sufficient to obtain converged CFD results.The time-averaged values, as well as the oscillating amplitudes of C L and C D , were regarded as indicators of the mesh resolution independence.The base size (l base in Fig. 9) for these cases was 0.5 m.Fig. 12(a) shows the results based on the mesh without near-foil refinement with different separate angles of wake refinement.The length of wake refinement was fixed at 60 m.When the separate angle is larger than 0.3 rad, the values of the force coefficients converge.Similarly, by fixing the separate angle of wake refinement at 0.3 rad, a few sets of mesh with different lengths of wake refinement can be studied, as shown in Fig. 12(b), indicating that the length of wake refinement should be greater than 60 m to exclude the influence of the mesh.Fig. 12(c) shows the effects of near-foil refinement based on a series of mesh sets having a 0.3 rad separate angle of wake refinement.The near-foil refinement does not greatly increase the total number of cells.Although the effects on C L and C D are not obvious, this local refinement was retained in subsequent analyses to reduce the need to change the mesh when varying α.In summary, the mesh follows a strategy of wake refinement with a separate angle of 0.3 rad and a length of 60 m, as well as near-foil refinement.
Two-dimensional simulations with different base sizes were carried out by following a certain refinement strategy.The variation of the base size affects the entire mesh, except for the prism layer mesh in the normal direction.As Fig. 13 shows, when the base size is less than 0.35 m, the difference in the force coefficients from the two most refined cases is less than 1% for C L and 4% for C D .Because C D is only approximately 8% of C L , the thrust force is mainly based on the lift.Because of the limitation of two-dimensional simulations that the flow separation cannot be well resolved, the two-dimensional mesh independence study did not go deeper.

Three-dimensional mesh independence
A series of three-dimensional meshes were generated based on the same refinement strategy.Five meshes were generated to study the three-dimensional mesh independence.For convenience, the meshes from coarsest to finance were numbered from 1 to 5. The force coefficient results for meshes 3, 4, and 5 were similar (Fig. 14).The difference in C L between meshes 3 and 5 was 0.37%, and the difference in C D is 0.36%.This suggests that when the mesh is more refined than mesh 3, the simulation results are irrelevant to the mesh.The number of cells in mesh 3 was 23,735,358.The following three-dimensional simulation results are all based on mesh 3.
As mentioned in Section 2.3.1, another three-dimensional simulation based on a polyhedral mesh with a similar total number of cells was also performed to exclude the influence of mesh typology.The differences between the force coefficients of the two types of mesh were less than 1%, so it can be concluded that the mesh typology did not affect the simulation results.

Performance evaluation
Predicting the performance of a wind-assisted ship propulsion system requires consideration of the full system, i.e., the sails and the ship.The F S and M yaw introduced by wind-assisted ship propulsion systems must be compensated for by a drift in the ship and rudder angle, both of which introduce added resistance, causing both a lower net thrust, i.e., the thrust of the sail minus the added resistance caused by the sail, and the need for sail trim optimization to achieve the best performance under any given constraints, e.g., the effects of the rudder or heel angle.Thus, a model with at least four degrees of freedom (4DOFs), i.e., surge, drift, yaw, and heel, must be used.The performance of a wind-assisted ship propulsion system also depends on the ship on which it is installed, which means that any performance or comparison study must involve a case study ship.Here, a tanker with a deadweight of approximately 100, 000 tons was used.The ship's dimensions are given in Table 2.
ShipCLEAN was used as a performance prediction tool in this study.ShipCLEAN is a generic model developed to provide accurate predictions with very little input data, such as those elements listed in Table 2, respecting 4DOFs.The remaining dimensions, such as L pp , depth, and superstructure dimensions, are estimated by the model using empirical formulas (Tillig et al., 2017).To predict the performance under real operation conditions, V TW , θ TW (see Section 3.3.3),currents, and water temperature are specified.Then wave heights and directions are specified or evaluated to match the wind speed, a specified fetch, and the wind direction (Tillig and Ringsberg, 2019).Detailed information about how the hydrodynamic forces and moments, including the hydrodynamic compensation of the side forces created by the wingsail, has been documented in previous publications (Tillig and Ringsberg, 2019;Tillig et al., 2017Tillig et al., , 2018)).In general, steady-state conditions are evaluated by the model, so the sum of all forces and moments in the 4DOFs must be zero, as shown in Equation ( 8).The model solves Equation ( 8) iteratively because forces and moments in different directions are H.Zhu et al. dependent on each other.
The achievable accuracy of the power prediction was demonstrated by Tillig et al. (2018).The development, applicability, and achievable accuracy of ShipCLEAN have been documented by Tillig andRingsberg (2019, 2020).The uncertainties in the prediction of sail forces and aeroand hydrodynamic interactions have been discussed by Thies and Ringsberg (2022).

Force coefficients
According to wind tunnel tests (Sheldahl and Klimas, 1981), when the Reynolds number is 1 × 10 7 , the critical angle of attack, α c , of the NACA 0015 foil is approximately 16 • , and the maximum C L is approximately 1.42.A previous study (Nikmanesh, 2021) indicated that, for the newly introduced crescent-shaped profile, α c is around 20 • .The following comparison between two-and three-dimensional simulations, as well as the comparison between different boundary conditions, are for the crescent-shaped foil.
For the crescent-shaped foil, with respect to force coefficients, the differences between two-and three-dimensional simulations with periodic top and bottom boundary conditions are 7.5% for C L and 58.1% for C D , as Fig. 15 shows.This suggests that two-dimensional simulations substantially overestimate the force coefficients even under low-α conditions, and C D is more sensitive than C L .In two-dimensional simulations, because of the limitation of the spanwise flow, the vortices are constrained and not well-developed (Park et al., 2017).Overestimation of force coefficients, especially C D , has also been observed in high-Re CFD simulations of airfoils with strong flow separation (Zhu, 2020).Comparing the results for different boundary conditions shows that these two cases produce force coefficients with obvious differences.When there is a freestream tip, the lift force on the foil decreases by 5.7% , and the drag force increases by 35.0%.
As Fig. 15 shows, the crescent-shaped profile has substantially higher force coefficients and is expected to provide better propulsive performance.Furthermore, because of tip vortices, the reduction in C L for the wingsail with the crescent-shaped profile is 5.7%, while the value is 10.7% for the sail with the NACA 0015 profile.Therefore, the effects of the tip vortices on the crescent-shaped profile are not as substantial as those on the NACA 0015 wingsail.On the other hand, because of flow separation, as discussed in Section 3.2.1, the crescent-shaped wingsail shows remarkable oscillation amplitude of the force coefficients.Therefore, the crescent-shaped sail is believed to suffer from more serious flutter, which increases the requirements of the sail structure.
To analyze the propulsive performance of the wingsail under various wind conditions, i.e., a large range of wind directions, for the proposed crescent-shaped foil, two-dimensional simulations with α values from − 2 • to 95 • were performed to determine α c , and 12 three-dimensional simulations were performed with different α values.The plots of the force coefficients are shown in Fig. 16, where the blue lines and symbols represent the two-dimensional results, and the red lines and symbols represent the three-dimensional results.Two-dimensional simulations predict much higher C D values than three-Fig.18. Time-averaged pressure coefficient distribution over the crescent-shaped profile.The inlet flow velocity is in the positive X direction.
H. Zhu et al. dimensional simulations, as Najjar and Vanka (1995) also observed.According to the two-dimensional simulations, the highest C D is approximately 3.7 when α = 80 • .For high-Re conditions, C D for a flat plate is approximately 1.98, and for a semicircle opening upstream, it is approximately 2.30 (Hoerner, 1976).The C D of the crescent-shaped profile should probably be between these two values.Furthermore, because of the tip vortices, C D is expected to be even lower, so the predicted C D is unreasonably high.In addition, C D suddenly decreases when α increases from 80 • to 90 • in the two-dimensional simulations.
These unreasonable results suggest that, because of the strong flow separation, two-dimensional simulations cannot provide reasonable predictions of the force coefficients when α > 20 • , as explained in Sections 3.1.2and 3.2.1.However, another assumption is still held that two-dimensional simulations correctly reveal the trend of the force coefficients.Based on this assumption, a hybrid two-and three-dimensional simulation method is proposed.First, a series of two-dimensional simulation cases with various α values were performed to determine the time-averaged force coefficients.Next, a limited number of three-dimensional simulations with distinctive α values, such as those that correspond to peak or valley values of force coefficients, were performed.Then, the ratios of C L and C D between the two-and three-dimensional simulations was calculated.Finally, C L and C D were rescaled based on the ratios.In this way, the computational capacity was reduced, because threedimensional simulations are much more computationally demanding than two-dimensional simulations.In Fig. 16, the purple dashed line

Pressure distributions
Fig. 17 presents the time-averaged pressure coefficient distribution along the surface of the profile.Because of the large camber of the crescent-shaped profile, the pressure difference between the pressure side and the suction side is obvious even when α = 0 • , as shown in Fig. 17(a).
When α ≤ 20 • , two-dimensional and three-dimensional simulations provide similar predictions of the pressure distribution.According to the three-dimensional simulations, the highest C L is obtained when α = 20 • .As Fig. 17(b) shows, the leading part of the profile, i.e., the upstream half, mainly contributes to generating the lift force because C p is quite low, which results in a large pressure difference.In addition, around the trailing edge (the chord-wise location is at approximately 7 m), the wall pressure at the pressure side becomes lower than that at the suction side.This inversion of wall pressure probably negatively affects the lift force coefficient, and the shape of the edges is expected to be optimized in the future.
However, when α is large, e.g., α = 35 • , as shown in Fig. 17(c), the two-dimensional and three-dimensional results are quite different.Based on two-dimensional simulations, the pressure on the suction side is unreasonably low, leading to unreasonably high force coefficients.
The same is observed when α is 60 • (Fig. 17(d)) and 85 • (Fig. 17(e)).One exception is that when α = 90 • , i.e., when the chord line is perpendicular to the inlet flow velocity, the suction-side pressure predicted by the two-dimensional simulations is even higher than that based on the threedimensional simulations.In the two-dimensional results, positive C p is widely distributed on the suction side.This strange phenomenon is regarded as the reason for the sudden decrease in C D in Fig. 16(a).The crescent-shaped profile has a sharp leading edge and a blunt trailing edge compared with conventional airfoils.Thus, C p is very high at the leading edge and negative at the pressure side when approaching the trailing edge.To improve the propulsive performance of wingsails with crescent-shaped profiles, the shape of the edges should be optimized.The proposed crescent-shaped profile is symmetric about the mid-chord to facilitate the operation of the wingsail, so multi-objective optimization of the aerodynamic performance of the edges should be conducted to give appropriate consideration to both edges.according to the two-dimensional and three-dimensional results, respectively.The areas in red represent the region where V X > V AW , i.e., where the streamwise velocity is higher than the inlet flow velocity, whereas those in blue represent the region where V X < V Aw .One of the most important characteristics of the flow field is the remarkable flow separation.For NACA 0015, there is usually no or only one flow separation point, depending on α.Flow separation does not on the suction side close to the trailing edge and also on the pressure side close to the leading edge.A large region with low velocity or even reversed flow appears attached to the pressure side.In Fig. 19(a), where α = 10 • , the positive and negative of wall shear stress converse, which indicates that flow separation occurs.When α = 20 • , a high-velocity area occurs near the suction side, causing low pressure and leading to multiple flow separation points on the downstream half of the profile, as shown in Fig. 19(b).When flow separation occurs in the boundary layer, a shear layer and a separated flow region between the shear layer and surface are formed, which modifies the pressure distribution.Flow separation is thought to increase pressure drag.

Flow separation
When α ≤ 20 • , the flow field predicted by the two-dimensional and three-dimensional simulations show some similarities.However, for large-α conditions, two-dimensional and three-dimensional simulations provide totally different results.Take the pair of simulation cases with α = 35 • as an example.In the two-dimensional results, the vortices are strong, but in the three-dimensional results, the vortices are difficult to recognize.In Fig. 20(d) and (e), where α are separately 25 • and 35 • for the two-dimensional simulations, a high-velocity reversed-flow region can be seen along the suction side.This high-velocity region causes low pressure, leading to a strong lift force.This is regarded as the main mechanism for the peak of C L in Fig. 16(b).However, this phenomenon is not seen in the three-dimensional results (Fig. 21(d) and (e)), in which V X still shows positive values in the area along the suction side, which explains why the peak of C L at α = 35 • is not as notable as in the twodimensional results.20(g)).These differences cause the uneven distribution of C p on the suction side in Figs. 17 and 18 (the blue lines and shadows represent the two-dimensional results).
However, three-dimensional simulations yield different results.Take the condition of α = 85 • (Figs.20(g) and 21(g)) as an example.In the three-dimensional results, the flow close to the suction side is almost quiescent at a very low velocity.Similar flow characteristics have been observed in similar studies.For instance, Castelli et al. (2012) used numerical simulations to study a flat plate with α = 90 • and found large low-velocity areas in the wake.This finding explains why the C p on the suction side predicted by the three-dimensional simulations is much higher than that predicted by two-dimensional simulations.In addition, in Figs.20(g) and 21(h), the downstream and upstream fields are both influenced by the wingsail profile.V X is always lower than V AW around the crescent-shaped profile.Hence, it can be inferred that when the wingsail is operated under α ≈ 90 • , i.e., under downwind conditions, the interaction among multiple wingsails might be considerable, which can be analyzed in future studies.From the iso-surface plot in Fig. 22, where colorful contours represent the streamwise vorticity, circular patterns of rotating air left behind the tip of the foil, which are the tip vortices, can be easily recognized (see Fig. 22(a) and (c)).Tip vortices cause a reduction in C L , as discussed in Section 3.2.2.Wingsails with different sectional profiles have different wake characteristics.For the crescent-shaped profile, vortex shedding is much more substantial.
Numerous vortex tubes can be seen in the wake region (Fig. 22(a) and (b)).However, for the NACA 0015 profile (see Fig. 22(c) and (d)), only limited vortex tubes can be seen developing on the suction side.Vortex shedding is regarded as the main reason for the oscillation of force coefficients.Thus, the oscillation of force coefficients of the crescent-shaped wingsail is much stronger than that of the wingsail with the NACA 0015 profile.
For the crescent-shaped profile, vortex shedding is obvious in the downstream region of the suction side and the trailing edge for simulations with both types of boundary conditions.However, in the downstream flow field, the structure of the vortex tubes is much more complex when the top and bottom boundaries are periodic.In Fig. 22(a), excluding the tip vortices, the vertical vortex tubes are uniform and easily recognized.Nevertheless, when the top and bottom boundaries are periodic (Fig. 22(b)), because of the constraints of the top and bottom boundaries, many horizontal vortex tubes can be seen.A similar phenomenon is also found in the NACA 0015 profile.This phenomenon explains the larger oscillation amplitude of the force coefficients.

Tip vortices
Another important characteristic of the flow field is the phenomenon of tip vortices, which is believed to be the main reason for the lift reduction when the boundary conditions are changed from a periodic top and bottom to a free tip and symmetric bottom.Fig. 23 presents the time-averaged wall pressure distribution at various section planes.When the top and bottom have periodic boundary conditions (Fig. 23  (b)), the pressure distribution along the foil is similar at different Z positions.For instance, the wall pressure distributions at the section planes of Z = 0.25H and Z = 0.95H do not show obvious differences.However, if the boundary conditions are a free tip and symmetric bottom (Fig. 23(a)), the pressure difference between the pressure side and the suction side decreases toward the tip of the wingsail, leading to a lift reduction.
In the streamline plot in Fig. 24, in which the color represents the spanwise velocity, flow is evident around the tip from the pressure side to the suction side.The vortices over the top of the wingsail are evident.The tip vortices affect the flow field around the tip over a very large area.The thickness of the wingsail at the mid-chord is 2 m, so, according to  For both sectional profiles, at the leading edge, two tip vortices develop independently at the pressure side and the suction side (Fig. 25).The tip vortex that develops at the pressure side, i.e., the tip leakage vortex, is considerably stronger than the other vortex.Hence, at around the mid-chord, these two tip vortices begin to fuse together.The tip vortices generated by the wingsail with the crescent-shaped profile are stronger and dissipate more slowly than those induced by the NACA 0015 wingsail.For the wingsail with the crescent-shaped profile, the structure of the vortex becomes complex at the tip closer to the trailing edge, as illustrated by the sectional plane X = 0.41L c in Fig. 25(a).In the wake region of the wingsail, at the sectional plane X = 0.63L c in Fig. 25 (a), a main large tip vortex is finally formed.However, for the wingsail with the NACA 0015 section, the main large vortex forms early around the sectional plane X = 0.91L c in Fig. 25 The tip vortices influence the flow field over a very large area, so proposing an effective solution to minimize the lift reduction caused by tip vortices is challenging.Some attempts have been made to achieve this, e.g., adding a wing flap or disc on top to prevent the round-tip flow illustrated in Fig. 24, but the results have not been ideal.Other tip design alternatives to optimize tip flow may be developed in future studies to increase C L .

Deep-stall conditions
When θ AW is close to 90 • , the wingsail may be operated with a large α, i.e., under deep-stall conditions.Fig. 26 shows that the distribution of ω X on the iso-surface of Q = 5 s − 2 .Tip vortices are still notable when α = 45 • (Fig. 26 More than one wingsail is usually installed on a ship to capture more wind power.In related studies (Malmek et al., 2020;Ouchi et al., 2011), the distance between the wingsails, i.e., the distance of the rotational axis, is usually less than two times the chord length.As Figs. 22 and 26 show, for wingsails with crescent-shaped profiles, the influence on the wake flow spreads more than four times the chord length in the downstream region.Therefore, it can be inferred that the interactions among multiple wingsails with crescent-shaped profiles are more substantial than those based on the conventional NACA series.A numerical study of the interaction of multiple wingsails may be conducted in the future.

Thrust for various wind directions
Based on the force coefficients from the CFD simulations with the freestream tip setup, the propulsive performance was evaluated, with the assumption that the force coefficients remain the same when the  H. Zhu et al. apparent wind speed changes following the direction of navigation, because the force coefficients are believed to be not sensitive to Re.
According to Lu and Ringsberg (2020), the best propulsive performance, i.e., the maximum C T , of a rigid wingsail is obtained when the point of the sail is a beam reach, which means that the apparent wind angle is approximately 90 • or 270 • .At this point of sail, the wingsail is always operated with the α that provides the highest C L .Thus, to compare the propulsive performance of a wingsail with the studied crescent-shaped profile and NACA 0015, the condition under which F L mainly contributes to the thrust is selected.Fig. 27 presents plots of C T vs. θ AW for the crescent-shaped profile and the NACA 0015, with C T calculated using Equations ( 3) and (7).For both profiles, α is fixed, which means that α = 20 • for the crescent-shaped profile and α = 16 • for NACA 0015.The crescent-shaped profile generates a noticeably greater thrust than the symmetrical NACA 0015 profile because of the higher C L .
It should also be noted that, usually, when the point of sail is luffing, close-hauled, or beam reach, i.e., θ AW is from 30 • to approximately 90 • , F L is the main source of thrust.However, under other conditions, the wingsail may be operated in another way to use F D .Therefore, to predict the propulsive performance for all apparent wind directions, an enumeration method is used for C T with θ AW in the range of 0   H. Zhu et al. and α in the range of 0 • to 90 • .The results for the propulsion (the highest C T at different θ AW ) and how the wingsail is operated (the α value that is applied to obtain the highest C T ) are plotted in Fig. 28.Because a polar diagram is always symmetric, only half of the polar plot, i.e., 0 • ≤ θ AW ≤ 180 • , is presented.
When the ship navigates against the wind, i.e., the luffing point of sail, C T is low or even less than 0. During that time, the wingsail is operated with a very low α to reduce the extra resistance.As θ AW increases, e.g., when θ AW = 60 • , C T increases, and the wingsail is operated with α c having the maximum C L .For θ AW in a wide range from 60 • to 180 • , the wingsail can provide appreciable propulsion.However, the wingsail is not operated in the same way.For instance, when 30 • < θ AW < 120 • , α should be approximately 20 • to achieve the maximum C L , and F L is the main source of thrust.When 120 • < θ AW < 150 • , i.e., the point of sail is board reach, the optimal α is approximately 40 • with C L ≈ C D , and the wingsail uses both F D and F L for propulsion.When the point of sail is running, i.e., when θ AW is approximately 180 • , F D is mainly used, and the wingsail is operated with α ≈ 80 • .
The rescaled two-dimensional results predict propulsive performance similar to the three-dimensional simulations.Therefore, the hybrid two-and three-dimensional method, which involves performing two-dimensional simulations together with a limited number of threedimensional simulations and then rescaling the force coefficients from the two-dimensional results to fit the three-simulation results, yields a  H. Zhu et al. reliable prediction of the propulsive performance while reducing the computational demand.

Expected fuel savings
The first step in calculating the expected fuel savings is to create polar plots of the fuel savings under different wind strength conditions over several true wind angles.Fig. 29 presents a polar plot of the fuel savings for the case study tanker with one crescent sail in 10 kn and 20 kn of wind.The sail is positioned 5 m behind the forward perpendicular at the centerline of the ship.
The polar plot shows that maximum fuel savings of approximately 9% at θ TW = 90 • in V TW = 10 kn and 25% in θ TW = 90 • in V TW = 20 kn can be expected.However, the fuel savings vary over the true wind angle.Thus, the performance must be predicted using actual routes and realistic weather conditions, as described in the following section.

Prediction of long-term fuel savings
To predict the long-term fuel savings, the AIS data of the case study ship were used to derive the position and speed during the year 2018.The environmental conditions were retrieved from the Copernicus Marine Environment Monitoring Service (CMEMS) and were updated every 3 h.The route is illustrated in Fig. 30.The wind conditions (true wind angle and true wind speed) are illustrated by the wind scatter plots in Fig. 31.Because this case is symmetrical, all apparent wind angles are between 0 and 180 • .
Total savings of 9.5% were achieved with the crescent sail.For comparison purposes, the simulations were repeated with a Flettner rotor (5 m in diameter and 30 m in height) positioned similarly.The Flettner rotor would have resulted in 9.8% savings on the same route.To better understand the bandwidth of potential savings during a year of operation, Fig. 32 presents a histogram of fuel savings respecting each waypoint, i.e., 1 point per 3 hours, for the full year.Additional drag was created on only a very few occasions.Approximately 34% of the time, the fuel savings were greater than 5%.The maximum additional fuel consumption was less than 1%.The maximum fuel saving was 97.5%.

Conclusions
This study focused on the analysis of a new crescent-shaped sail profile based on high-fidelity numerical simulations.The computational method used was an unsteady RANS analysis with the k-ω SST model.The numerical settings, including the boundary conditions, were investigated to understand their effects on the prediction accuracy of the sail aerodynamics.
It was found that the two-dimensional simulations generally overestimated the force coefficients compared to three-dimensional simulations.This might be because vortex evolution and coherence in the spanwise direction are excluded from two-dimensional simulations.Another reason might be that the tip vortices due to the freestream tip cause a lift reduction.However, both two-and three-dimensional simulations showed similar trends in the force changes with respect to the angles of attack, so by rescaling the two-dimensional results to a limited number of three-dimensional simulation cases, the propulsive performance was predicted with reasonable accuracy.
To understand the tip vortex effects on aerodynamics, such as F D and F L , the computational domain and boundary conditions were investigated using two setups.One setup used periodic boundary conditions at the spanwise side boundaries, whereas the other used symmetric boundary conditions at both side boundaries.In the second case, the topside boundary was positioned far away from the top-side edge of the sail to eliminate the influence of the boundary on the flow.The sail configuration was attached to the bottom side boundary, representing the water-free surface, and was thus imposed with the symmetry boundary condition.The second case reproduced the vortices induced   H. Zhu et al.

Fig. 6 .
Fig. 6.Two sets of boundary conditions and computational domains for the three-dimensional simulations.

Fig. 7 .
Fig. 7.The computational domain of the three-dimensional simulations.

Fig. 13 .
Fig. 13.Two-dimensional mesh independence study based on different base sizes.

Fig. 14 .
Fig. 14.Three-dimensional mesh independence study based on different base sizes.

Fig. 15 .
Fig. 15.Force coefficients of NACA 0015 and crescent-shaped foils from two-and three-dimensional simulations with different boundary conditions.The error bars represent the amplitudes of oscillation.

Fig. 16 .
Fig. 16.Time-averaged C L and C D curves.

Fig. 17 .
Fig. 17.Time-averaged C p distribution along the wall with different α.Three-dimensional results are based on the condition with a freestream tip.

Fig. 16
Fig. 16(a) reveals that there is no clear α c for this type of profile.Two peaks of C L occur when α = 20 • and α = 35 • .Based on two-dimensional simulations, the highest C L is approximately 2.7 when α = 35 • , whereas based on three-dimensional simulations, the highest C L is approximately 1.9 when α = 20 • .Both two-dimensional and three-dimensional simulations show that when α ≤ 80 • , C D increases as α increases (Fig. 16(b)).

Fig. 19 .
Fig. 19.Wall shear stress coefficient distributions along the wall of the wingsail.For each subfigure, left: entire profile; right: near the trailing edge.
Fig. 18 presents the vectors of the time-averaged wall pressure coefficient.For low-α conditions, such as α = 20 • (Fig. 18(b)), a strong suction force occurs on the suction side close to the leading edge.Since the normal direction of the wall surface in this area is close to the positive Y direction, this suction force efficiently contributes to the lift force.
Fig. 19 presents the wall shear stress distribution along the wall surface of the sectional profile for different α values.Figs.20 and 21 show the flow field around the sectional profile for different α values

Fig. 20 .
Fig. 20.V X distribution and velocity direction vectors around the foil from two-dimensional simulations for different α values.The inlet flow velocity is in the positive X direction.
Fig.20(g)).These differences cause the uneven distribution of C p on the suction side in Figs. 17 and 18 (the blue lines and shadows represent the two-dimensional results).However, three-dimensional simulations yield different results.Take

Fig. 21 .
Fig. 21.V X distribution and velocity direction vectors around the foil from three-dimensional simulations for different α values.Z = 0.5H.The inlet flow velocity is in the positive X direction.

Fig. 24 ,
Fig. 24, spanwise flow occurs within more than 10 m around the tip.For both sectional profiles, at the leading edge, two tip vortices develop independently at the pressure side and the suction side (Fig.25).The tip vortex that develops at the pressure side, i.e., the tip leakage vortex, is considerably stronger than the other vortex.Hence, at around the mid-chord, these two tip vortices begin to fuse together.The tip vortices generated by the wingsail with the crescent-shaped profile are stronger and dissipate more slowly than those induced by the NACA 0015 wingsail.For the wingsail with the crescent-shaped profile, the structure of the vortex becomes complex at the tip closer to the trailing edge, as illustrated by the sectional plane X = 0.41L c in Fig.25(a).In the wake region of the wingsail, at the sectional plane X = 0.63L c in Fig.25(a), a main large tip vortex is finally formed.However, for the wingsail with the NACA 0015 section, the main large vortex forms early around the sectional plane X = 0.91L c in Fig.25(b).A comparison of the Fig. 24, spanwise flow occurs within more than 10 m around the tip.For both sectional profiles, at the leading edge, two tip vortices develop independently at the pressure side and the suction side (Fig.25).The tip vortex that develops at the pressure side, i.e., the tip leakage vortex, is considerably stronger than the other vortex.Hence, at around the mid-chord, these two tip vortices begin to fuse together.The tip vortices generated by the wingsail with the crescent-shaped profile are stronger and dissipate more slowly than those induced by the NACA 0015 wingsail.For the wingsail with the crescent-shaped profile, the structure of the vortex becomes complex at the tip closer to the trailing edge, as illustrated by the sectional plane X = 0.41L c in Fig.25(a).In the wake region of the wingsail, at the sectional plane X = 0.63L c in Fig.25(a), a main large tip vortex is finally formed.However, for the wingsail with the NACA 0015 section, the main large vortex forms early around the sectional plane X = 0.91L c in Fig.25(b).A comparison of the (a)), and vertical vortex tubes can still be seen developing on the suction side, but vortex tubes extending in the horizontal direction can also appear to be induced by the trailing edge.When α = 75 • (Fig. 26(b)), the tip vortices are much weaker and dissipate more quickly.Fig. 26(c) illustrates the conditions for α = 90 • , i.e., when the chord line of the sectional profile is perpendicular to the apparent wind.Under this condition, tip vortices and vertical vortex tubes almost disappear, whereas horizontal vortex tubes induced by the edges become the main characteristic of the wake flow.

Fig. 23 .
Fig. 23.Time-averaged wall pressure coefficient distribution of the crescent-shaped profile from three-dimensional simulations with different boundary conditions.α = 20 • .

Fig. 25 .
Fig. 25.Q at different streamwise positions around the tip.The boundary condition is a freestream tip and symmetric bottom.α = 20 • for the crescent-shaped profile and α = 16 • for NACA 0015.The inlet flow velocity is in the positive X direction.

Fig. 26 .
Fig. 26.Iso-surface of Q = 5 s − 2 under deep-stall conditions, colored with ω X .The inlet flow velocity is in the positive X direction.

Fig. 29 .
Fig. 29.Polar plot of relative fuel consumption savings vs. θ TW at two V TW .

Fig. 31 .
Fig. 31.Scatter plot of the experiences with apparent wind angles and true wind speeds.

Table 2
Main dimensions of the case study ship.