Rapid aerodynamic method for predicting the performance of interacting wing sails

Rapid performance prediction tools are required for the evaluation, optimization, and comparison of different wind propulsion systems (WPSs). These tools should capture viscous aerodynamic flow effects in 3D, particularly the maximum propulsion force, stall angles, and interaction effects between the lift-generating units. This paper presents a rapid aerodynamic calculation method for wing sails that combines a semi-empirical lifting line model with a potential flow-based interaction model to account for 3D interaction effects. The method was applied to a WPS that consisted of several wing sails with considerable interaction effects. The results were compared to CFD RANS simulations in 2D and in 3D. For the evaluated validation cases, the interaction model improved the prediction considerably compared to when the interaction was not accounted for. The method provided acceptable driving force, moments, and stall predictions, with negligible computational cost compared to 3D CFD simulations.


Introduction
Global trade has depended on shipping for centuries, to the extent that shipping has become a synonym for freight transport in general.Currently, the global fleet of ships continues to rely almost exclusively on fossil fuels; the pathway toward its decarbonization in this field remains undetermined.In fact, long-haul shipping is still largely considered a hard-to-abate sector and is at risk of contributing an increased share of global emissions, as other industry sectors transition to sustainable energy sources.Electrification, fuel cells, nuclear power, or new types of fuels could offer carbon-neutral energy on certain trading routes, but each of these technologies faces several challenges that are yet to be resolved (DNV, 2021).
The reintroduction of wind propulsion for commercial ships, which directly utilizes the wind to drive the vessel forward, offers a large potential to reduce emissions from global shipping (Chou et al., 2020).If used in combination with other propulsion systems, it can extend the sailing range and lower the energy consumption both of new-builds and of retrofitted existing ships.In a world with volatile energy prices and limited global availability of several energy carriers, the possibility of directly harvesting wind energy on board decreases operational costs (Gerhardt et al., 2022), thereby reducing financial risks as well as emissions.Modeling aerodynamics using the classical inviscid potential flow theory is an option, as it has high computational efficiency.Methods based on Prandtl's lifting line theory (LLT; see Katz and Plotkin (2001) or Anderson (2017)) have been widely used in the aeronautic industry for the last hundred years.The strength of LLT is its simplicity in relation to potential flow-based panel methods and viscous simulation tools, such as computational fluid dynamics (CFD).Several adaptions of the theory have been made -for example, Phillips and Snyder (2000) and Cheng and Wang (2018) extended the classical LLT to several lifting surfaces in their study of tandem airplane wings; they concluded that the LLT-based methods could predict the inviscid loads generated by high-aspect-ratio lifting surfaces at low angles of attack with an accuracy as good as that obtained from modern panel codes or from CFD. Spall et al. (2013) and Bordogna et al. (2018) successfully used LLT to model interacting thin sails for pre-stall angles.The LLT-based methods have, along with other potential flow methods, the downside of not including viscous effects, such as viscous drag and stall, and cannot be used to predict maximum lift.
CFD and wind tunnel testing are the two main tools for simulating interacting sails.For example, Nakashima et al. (2011) used large eddy simulation CFD and Ouchi et al. (2011) used Reynolds-Averaged Navier-Stokes (RANS) CFD to model interacting WPSs in 3D.Further, 2D RANS CFD has been combined with different optimization algorithms to identify the optimal sheeting angles (see Chapin et al. (2006) and Lee et al. (2016)) and to optimize parametric design of suction wing sails with active circulation control (Cairns et al. (2021)).This 2D CFD optimization approach, which relies on surrogate models and response surfaces, does not consider the 3D effects on the sail-sail interaction and sheeting optimization.Extending the 2D CFD optimization to full 3D CFD would come at a (too) high computational cost.The list of variables is long: a large number of 3D simulations is required to evaluate the effects of different WPSs, deck placements, aerodynamic profiles, heights and planforms for a large number of wind conditions and apparent wind angles.
Several WPS concepts have been tested in wind tunnels by for example Ingham and Tersløv (1985), Fujiwara et al. (2005), Li et al. (2015) and Bordogna et al. (2018).Marimon Giovannetti et al. (2022) performed wind tunnel tests for the interacting wing sails of the Oceanbird.Wind tunnel testing remains an important tool for validating simulations and final concepts, but it generally comes at too high a cost to be used to evaluate multiple concepts at early design stages or for optimization.Additionally, certain aerodynamic effects may not be accurately captured at model scale.
Due to the shortcomings of the aerodynamic performance prediction tools discussed above, there is a need for new methods.Graf et al. (2014) presented an efficient method for modeling a single wing sail by combining a 3D non-linear lifting line model with viscous 2D data.Based on their study, Persson et al. (2019) developed another lifting line algorithm called SILL, which assumes an elliptical lift distribution along the span.2D CFD data were used as input, and the methods were compared with 3D RANS CFD simulations.It was concluded that the SILL method had good potential to predict both generated forces and the point of stall for a single wing sail.
However, Persson et al. (2019) did not consider the interaction effects between several sails.Such effects may be considerable, as indicated, for example, by the wind tunnel tests with interacting rigid wing sails performed by Bordogna et al. (2018).Changes in sail loading due to such interaction were seen when the sail distances ranged from one to three times the chord lengths.The optimal sheeting angles also changed; in their study of the performance of multiple sails, Li et al. (2015) found that the thrust of the interacting sails increased by 42.3% when the sails were sheeted independently.
Thus, a computationally efficient aerodynamic method that considers both the viscosity and the sail-sail interaction effects is needed.Such a method, where a 3D potential flow-based interaction model is combined with the SILL code, was first outlined by Malmek et al. (2020) in a preliminary report and subsequently detailed in her Licentiate Thesis (Malmek, 2023).Moreover, Tillig and Ringsberg (2020) also suggested a potential flow-based method to model interaction effects for rotor sails in the 2D horizontal plane.However, this method cannot be applied to wing sails.
This paper describes the 3D sail-sail interaction model proposed by Malmek (2023).The 3D interaction model is combined with the single-sail code SILL to model multiple interacting wing sails.The new combined method is called Interacting Sectionally Integrated Lifting Lines (ISILL).The ISILL method is described and validated against RANS CFD simulations.
The setup of the validation study is presented in Section 2. This is followed by a description of the proposed aerodynamic method in Section 3. The validation results are reported and discussed in Section 4. The conclusions and recommendations for further research are presented in Section 5.

Methodology
This section describes the setup of the study.Its main focus is to validate the ability of the ISILL method to improve the sail-sail interaction modeling compared to SILL.The test case, described in Section 2.1, is delimited to wing sails standing in a row.Section 2.2 presents the selected validation cases, covering upwind to beam reach apparent wind angles.The loads predicted by ISILL are compared to RANS CFD and the setup of the CFD simulations are described in Section 2.3.A further discussion of the study limitations is provided in Section 4.5.

Test case geometry and definitions
The WPS test case in this paper is based on the Oceanbird research project (2018)(2019)(2020)(2021)(2022), where a concept for a wind-powered car carrier was developed by Wallenius Marine, SSPA (now RISE), and the KTH Royal Institute of Technology.A conceptual design of Oceanbird is presented in Fig. 1, with several wing sails standing in a row.The length of the ship is 200 m and the beam is 40 m.
The wing sail used as a test case has a slightly simplified platform compared to the conceptual design (see the left-hand side of Fig. 2).The section shape is a NACA0015 airfoil with a modified thinner trailing edge.The wing sail geometry is equal to that of the validation case used by Persson et al. (2019) when developing the SILL method.Each sail has a geometric planform area () of 1,844 m 2 , a mean chord length () of 23 m, and an aspect ratio () of 3.47.The sails stand in a row at the ship's centerline with an equal leading edge to leading edge distance () of 43.2 m.The center of rotation is fixed at the quarterchord (i.e., 25% of the sail chord from the leading edge).The sails are assumed to stand on an infinite flat plate, which is treated as a symmetry plane.Therefore, the the effective span of the wing sail () is doubled compared to the sails' height.
Figs. 2 and 3 depict the angles, forces, and different coordinate systems.The sails are indexed from the foremost sail aftwards.The angle of attack   , (the subscript refers to the sail number), is defined as the difference between the apparent wind angle   and the sheeting angle, .The ships' coordinate system is defined with the  axis in the forward direction and the  axis positive to port.The origin is at the midpoint between the rotation centers of the first and last sails.The  ′  ′  ′ −coordinate system is related to each sail and has its origin at the quarter-chord.Moreover, the  ′ -axis is parallel to the direction of the incoming wind.Fig. 3 presents the total driving force coefficient,   , and total side force coefficient,   , generated by the whole WPS.The longitudinal center of effort (LCE) is located where the forces of all sails combined apply and is measured forward of the origin of the ship's coordinate system.
All the force coefficients are defined by dividing the corresponding force by 1 ∕2 2 ∞ , where  is the density.In 2D, corresponding to the −plane, the sail area is replaced by the chord length.The moment coefficients are defined by dividing the moment by 1 ∕2 2 ∞   , where   is the number of sails.The moment coefficients are calculated in relation to the coordinate system of the ship and the sail forces are assumed to apply at the quarter-chord.The force coefficients are denoted by capital letters in 3D and by lowercase letters in 2D.

Selected validation cases
Table 1 provides an overview of the validation matrix.The ISILL model was evaluated both in 2D (with three sail profiles in a row) and in 3D (with four wing sails in a row).The 2D evaluation was performed in the horizontal plane since the interaction effect was strongest in this plane and it isolated the interaction of the bound vortices.By limiting the simulations to 2D, it is possible to perform multiple CFD simulations at a low computational cost.This enables the performance of several sheeting sweeps with an increasing angle of attack, from preto post-stall.Three different apparent wind conditions were studied: (1) sailing upwind close-hauled at   = 30 • ; (2) sailing upwind on a close reach at   = 60 • ; and (3) beam reach at   = 90 • .The apparent wind speed ( ∞ ) was set at 10 m/s, thereby yielding a mean Reynolds number of ≈ 15 ⋅ 10 6 .

CFD simulation setup
In this study, CFD simulations are used in two ways: (1) as 2D input   ∕  data to SILL/ISILL and (2) as validation for ISILL.An overview of the different CFD simulations is presented in Table 2.In 2D, single sail computations were first performed for CFD verification and validation.This step was followed by generating a   ∕  -table for model input.Simulations with three sails were used to validate ISILL in 2D.In 3D, the CFD settings were first verified and then used to generate the validation data with four sails.
All simulations were run with Simcenter STAR-CCM+ (Siemens, 2022), using the steady RANS-equations along with Menter's - shear stress transport (SST) model.The grid verification studies were conducted using the numerical uncertainty analysis tool provided by MARIN based on studies by Eça and Hoekstra (2014) and Eça et al. (2019).
RANS CFD simulations were selected due to their relative computational efficiency.In general, the accuracy of steady RANS simulations is low beyond stall.However, from a sail performance prediction perspective, this area of the lift curve is of less importance since it is located above the optimal sheeting angle.The 2D CFD validation below shows that pre-stall forces and the point of stall are well predicted.

2D CFD simulations
In the 2D CFD calculations, the domain had four boundaries.The incoming wind was parallel to the top and bottom boundaries, which were treated as symmetry planes.The airfoil profiles were initially placed horizontally and then rotated to the correct angle of attack.A low y+ wall treatment and the gamma transition model were used.
Fig. 4 presents the convergence of the lift coefficient with the grid density.The step size on the horizontal axis is relative to the finest grid (Grid 1) and is defined as [(No. of cells Grid 1)/(No. of cells Grid )] 1∕2 .The number of cells in Grids 1 to 5 were 1,119,542; 591,772; 309,864; 171,743; and 96,096.The computed uncertainty bars for each grid are also depicted in the figure.The finest grid had an uncertainty of 1.8%; however, Grid 3 was deemed sufficiently accurate, with an uncertainty of 4%.For the viscous drag,   , the resulting numerical uncertainty for Grid 3 was 17.5%.This is a large uncertainty and was caused by the fact that the verification was performed at  = 15 • , which is close to stall.However, the viscous   is at least one order of magnitude smaller than   and has little effect on the driving force   .
The CFD computations were validated against 2D wind tunnel experiments (Jacobs et al., 1933) for NACA0015 at Reynolds number of 3.2 ⋅ 10 6 .The result is presented in Fig. 5.The error in the maximum lift coefficient is 6% and the angle of maximum lift is overpredicted by 1 • .For lower angles of attack, the RANS simulations and the experimental data agree very well.Beyond stall, the error increases drastically, which is expected for RANS simulations, as the flow becomes highly unsteady.
The same set up as that for the the single sail was used to generate data for the 2D ISILL validation cases (see Table 1).The selected mesh had a maximum surface cell size of 0.0025 at the profile and 36 Fig. 5.The 2D RANS CFD simulations compared to wind tunnel test (WTT) data (Jacobs et al., 1933) for the NACA0015 airfoil at Re = 3.2 ⋅ 10 6 .prism layers.Additionally, the area between the sails was refined to a maximum cell size of 0.025.The three profile simulations had a mesh with ca 390,000 cells.

3D CFD simulations
In the 3D simulations, the domain has six boundaries.The top and bottom (floor) of the domain were set as symmetry planes, and the four sides of the domain as velocity inlets and pressure outlets.A high y+ wall treatment and no transition model were used.
The grid verification, presented in Fig. 6, resulted in a numerical uncertainty of 1.9% for   and 5.4% for   for the selected Grid 3. The step size on the horizontal axis is relative to the finest grid (Grid 1) and is defined as [(No. of cells Grid 1)/(No. of cells Grid )] 1∕3 .The number of cells in the grids from Grid 1 to Grid 5 were 25,424,790; 17,450,590; 12,194,230; 8,565,199; and 6,090,272.The same setup was used for generating the 3D validation data with four sails.The mesh had approximately 43,000,000 cells and 11 prism layers.There were several refinement areas, as displayed in Fig. 7.A single simulation required ca 5 h to converge using 40 nodes on the high-performance computing cluster Tetralith (provided by the National Academic Infrastructure for Supercomputing in Sweden).Each node has 2 CPUs (2x Intel Xeon Gold 6130) with a total of 32 cores, 384 GB RAM, and a 960 GB local disk.

Aerodynamic interaction method
This section presents the proposed aerodynamic method, ISILL, for modeling interacting wing sails.The single-sail code SILL (Persson et al., 2019) is extended with a potential flow vortex model for the interaction effects between the sails.Moreover, an optional viscous correction of the model is suggested to manage boundary layer effects.
The SILL method for a single sail is introduced in Section 3.1.The sail-sail interaction model is described in Section 3.2.A correction to the effects of a changed boundary layer development caused by the interaction is proposed in Section 3.3.Finally, the steps in the full ISILL method, which combines SILL and the interaction model, are summarized in Section 3.4.

The single sail model (SILL)
In SILL, the sail is divided spanwise into several strips.Based on the wind conditions at each strip, corresponding   and   are obtained from the pre-tabulated 2D data.This enables the calculation of the related forces that act on each strip.Since the local force components obtained from the coefficients are oriented relative to the local angle of attack   , new components are obtained that represent the lift and the drag force relative to the average angle of attack of the sail.These 2D forces are then integrated to obtain the total sail force, denoted by  2 and  2  , to indicate that they do not consider the 3D effects.Then, the 3D   is estimated from  2 by assuming an elliptical lift distribution over the sail (Anderson, 2017): Based on the lift and aspect ratio, the induced drag    may be calculated as The induced drag is added to the integrated viscous drag to obtain the total drag: (3) The method includes approximations, particularly the assumption of an elliptical lift but has been shown to yield acceptable results for a single wing sail (Persson et al., 2019).Fig. 8 presents the prediction of the lift and drag forces for a single wing sail.The SILL results are compared with 3D RANS computations.
An alternative to using SILL for estimating the 3D effects on the single sail is to include the vortex system of the sail itself in the interaction model presented below.This approach could improve the prediction of the stall angle; however, it increases the computational cost and the risk of numerical instability.

The sail-sail interaction model
The iterative sail-sail interaction model is explained below.Each sail is discretized in strips with a corresponding horseshoe vortex, as illustrated in Fig. 9.The vertical vortex filament of each horseshoe vortex, the bound vortex, is placed at the quarter-chord.The horizontal vortex parts of each horseshoe vortex-the free vortex-are, in theory, infinitely long, but have a finite length in the model.The free vortices follow the free stream direction.
To determine the strength of each vortex, an elliptical lift distribution1 is assumed.Based on the total lift coefficient, the maximum strength of the circulation at mid sail span,  0 , can be calculated in the following manner (Anderson, 2017): where |  | is the average local wind strength over all strips.The total circulation strength, summing the contribution from all bound vortices, at a strip  is calculated as where ℎ is the distance to the center of the strip from midspan.The free vortex that originates at the intersection between two strips  − 1 and  has the strength  () =  ( − 1) −  ().
Modeling the sails in the above manner, the total induced velocity from the surrounding sails on each strip  on a sail  can be calculated.This implies that a local wind condition   (, ) can be determined [i.e., the global wind conditions experienced by a sail in undisturbed flow,  ∞ (), plus the induced flow due to interaction,   (, )].The model for the four sails is exemplified in Fig. 10.According to the Biot-Savart law, the induced velocity   at point  at distance  from the center of a segment  of a straight vortex is given by where  is the vortex strength.By using the definitions in Fig. 11 and by integration, Eq. ( 6) can be modified (Katz and Plotkin, 2001, p. 38-41) to the following form: where  0 =  1 −  2 .This form enables the calculation of the the induced velocity of a straight vortex segment between two points and is used in the ISILL code.
In theory, the free vortices behind the sails are of semi-infinite length, but to use Eq. ( 7), they are assumed to have a finite length of 400.The total interaction induced velocity   (, ) in a point is calculated by summing the contributions from all bound and free vortices that represent the surrounding sails.Note that the SILL approximation is used to consider the finite (3D) sail effects on sail .Therefore, the change in local wind angle and magnitude induced by the wake sheet behind sail  itself is not included in the calculation of   (, ).
As mentioned, the bound vortices are located at the quarter-chord.This is a correct choice for airfoil sections in an undisturbed flow since that is where the lift acts for thin symmetric sections.However, in a flow disturbed by other sails, this is not necessarily the best choice.To identify the best position, several 2D cases with three sails were run.The effect of interaction on the pressure distribution and the position of the force center was studied.Moreover, the velocity field generated by the sail using potential flow was compared with CFD for different bound vortex positions in ISILL.The conclusion was that the quarter-chord remains the best choice for the vortex position.
To use the pre-calculated   () and   () data tables, an angle of attack is required.For an airfoil in a disturbed flow, where the local velocity changes its direction along the chord, the angle of attack is not obvious.Where should the angle be determined?Having fixed that point, it appears natural to also use the velocity magnitude at that point to evaluate the forces.To optimize the position of the point, referred to as  , in Fig. 10, systematic computations with varying positions were conducted.It turned out that a suitable position is at the threequarter chord (75% of the chord from the leading edge).Incidentally, this is also the point at which the sum of the induced velocity from the bound vortex on a single airfoil in an undisturbed flow and the free stream velocity is parallel to the chord line (Weissinger, 1947).

Boundary layer effects and postponed stall (ISILL+BL)
When two lifting airfoils are close together, such that the trailing edge (TE) of the forward foil is in the low-pressure zone on the suction side of the aft foil, the adverse pressure gradient on the suction side of the forward foil is reduced.The boundary layer (BL) then grows more slowly and is less prone to separation.Consequently, larger angles of attack may be attained before the sail stalls.Therefore, this section suggests a correction for these effects.The ISILL method with this correction is called ISILL+BL.
In its basic implementation, ISILL does not consider the interaction effect on the boundary layer, as it uses the local wind conditions at the control point to look up pre-calculated lift and drag coefficients.These coefficients have been established in an undisturbed environment.When the sails are oriented in a row, this causes an underprediction of the maximum lift force and stall angle in upwind conditions.In a few of the 2D validation cases, the maximum lift force on the forward sail was found to be underpredicted by approximately 10% by ISILL.
Here, we suggest a boundary layer correction to ISILL by adjusting the 2D lift curve   () when the TE is near the low-pressure region of the following sail.This is done by blending a theoretical inviscid lift coefficient    and the lift coefficient obtained from the pre-calculated table,    .The inviscid lift coefficient is based on the thin airfoil theory: According to the inviscid potential flow theory, the sail does not stall.Therefore, the stall angle and maximum lift force can be increased by interpolating between the pre-tabulated viscous coefficient and the inviscid coefficient: where the relation between the fractions is   = 1 −   .However, a control mechanism for the interpolation is required.
It is assumed that the fraction of the inviscid lift,   , can be determined based on the change in the adverse pressure gradient on a sail profile caused by the interaction.To obtain a metric for this change, the pressures at the TE and at the three-quarter-chord chord are utilized.The difference between these pressures is a measure of the rear part's pressure gradient.By computing the change in this pressure difference due to the surrounding sails, it is possible to obtain a metric for the interaction effect.This metric is denoted as   and calculated in the following manner: The induced pressure coefficient   is obtained from where   =   +  ∞ .Only the induced velocities created by the surrounding sails are considered when calculating   .
Although there are good physical reasons for a coupling between   and   , the exact relation is unknown.To determine this, a reference upwind case where the postponed stall effect is strong was selected (v.c.2.1).  was then adjusted to yield a good correspondence between CFD and ISILL+BL; see Fig. 12. Setting    to 0.45 yields a balanced correction, which implies that 45% of the lift coefficient comes from the inviscid    and 55% from the viscous    .The corresponding metric for the reference case is    .The assumption is that the   for any case can be established by where  is an unknown function.
The simplest relation is a linear dependence.This yields a considerable improvement in the prediction of maximum   and the corresponding stall angle.However, the connection appears to be stronger; thus a quadratic dependence was tested.The results were then further improved, and the following relation was adopted: To avoid unrealistic corrections,   is limited to the range 0.0-0.45 and the change in stall angle due to the correction is limited to a maximum of 5 • .In this study, the correction reference values,    and    , have been calibrated in 2D.The same reference values are applied in 3D.

Full iteration loop
Because the SILL method requires the local wind conditions at each strip to estimate   and the local wind conditions are updated by the vortex interaction model, the final   needs to be calculated iteratively.To initiate the iterations, the induced velocities of the surrounding sails are assumed to be zero.Steps 1 and 2 below are then repeated until the induced velocities have stabilized.
1.The following steps are repeated for each sail  to establish the circulation, lift force, and drag force: (a) For a set apparent wind angle   () and sheeting angle of sail ,   , the angle of attack  , is calculated for each strip .(b) The angle of attack is corrected by the induced velocities of the surrounding sails, calculated in Step 2. In the first iteration step, the induced velocities are 0. The new local angle of attack,   (, ), is then used to retrieve the local force coefficients in the 2D data table:   (, ) ⇒   (, ),   (, ). 2 (c) Following the steps in the SILL method,   () and   (𝑖) of the 3D finite sail are estimated by integrating the 2D forces and applying the assumption of an elliptical lift distribution (Eqs.( 1)-( 3)).(d)  0 () is calculated with Eq. (4) using the sail average local wind speed |  |.(e) The circulation strength at each strip,  (, ), is estimated by assuming of an elliptical lift distribution (Eq.( 5)).
2. To update the induced velocity at point  , at the three-quarterchord of each strip  on each sail , the following steps are taken: (a) For each sail, the strengths of the horseshoe vortices representing the sail are determined based on the circulation strength established in Step 1e.
2 If the boundary layer correction introduced in Section 3.3 is applied, the inviscid lift coefficient    is calculated using the local angle   (𝑖, 𝑗).The interpolation between    and the pre-tabulated coefficient,    , is determined by Eqs. ( 9)-( 13).
(b) The total induced velocity,   (, ), is calculated using Eq. ( 7) by summing the contributions from all surrounding sails.When applied to the study test case, the ISILL iteration procedure above converges within a few seconds on a standard computer.

Results and discussion
In this section, the comparison of the ISILL method to CFD simulations is presented and discussed.First, the effect of the boundary layer correction is highlighted in Section 4.1.Section 4.2 presents the 2D stall validation study with several sheeting sweeps at different   .The 3D results are divided into two parts: Section 4.3 presents the results for fixed sheeting at pre-stall angles and Section 4.4 presents the 3D stall study at   = 30 • .

Boundary layer correction
Figs. 12 and 13 illustrate the effect of the suggested boundary layer correction.In the figures, the predictions without (ISILL) and with (ISILL+BL) boundary layer correction are compared with the CFD results and the results using SILL.
Fig. 12 depicts the variation in lift coefficient for Sail 1 for three different apparent wind angles.Sail 1 is swept whilst the other two sails are stationary [see Table 1 (validation cases 2.1, 2.4, and 2.7)].As expected, the strongest boundary layer and postponed stall effects are seen at   = 30 • .At   = 90 • , the interaction effects are weaker and there is no need for correction, as the trailing edges of the forward wings are away from the low pressures of the following sail.Compared to SILL, there is a large improvement in the lift predictions even without the boundary layer correction; but when the correction is included, a further improvement is noted for the upstream sails around the maximum lift angle when sailing upwind.This is evident in Fig. 13, which depicts the lift coefficient for each of the three sails at   = 30 • when sweeping with Sail 3 (v.c. 2.3).
The boundary layer correction has a physical foundation: the influence of the adverse pressure gradient on separation.However, it requires CFD or other data for calibration for a given geometry.On the other hand, only a few 2D CFD cases are required to obtain a correction applicable for all sheeting angles and all possible positions of the 3D wings relative to one another.The quadratic relationship in Eq. ( 13) is not well founded and a linear relation also yields an improvement.In future research, further investigations into this relationship should be made.

Sheeting sweeps up to the point of stall in 2D
Figs. 14-16 present the total 2D forward driving force coefficient and yaw moment coefficient for the 2D validation cases defined in Table 1.In each figure, one of the three sails is swept from a low to a high angle of attack past the optimal sheeting angle at which the maximum driving force is generated.The other two sails are fixed, set at sheeting angles close to the maximum lift force so that these two stationary sails generate a strong interaction effect, thereby altering the velocity field around the sweeping sail.
In the figures, the blue lines show the results for the ship sailing at an apparent wind angle of 90 • ; the black lines for 60 • ; and the red lines for 30 • .All the CFD results are presented with full lines, the ISILL results with dash-dotted lines, the ISILL+BL results with dashed lines, and the SILL results with dotted lines.
Studying the difference between ISILL and SILL-that is, with and without considering interaction-the difference increases when the ship points higher toward the wind.Sailing upwind also generates the largest differences in moments.The highest driving force is generated for beam reach (  = 90 • ).
For all test cases, the interaction model improves the prediction of the driving force coefficient.The model improves the predicted absolute value and angle of attack for which the maximum value of   is reached.The largest effects of the interaction are seen in the two upwind cases (  = 30 • and   = 60 • ), which ISILL+BL predicted.In general, ISILL without compensating for the boundary layer effects also improves the prediction of the driving force around stall compared to SILL, but not as much as ISILL+BL.The force prediction at beam reach is good, but this condition has the weakest interaction effects.There are a few cases where the predicted driving force deviates from the CFD results, see the sweep with Sail 2 in Fig. 15(a).For Sail 2, some of the interaction effects from the sail in front and the sail behind cancel each other out, but the flow still differs considerably compared to a sail in undisturbed flow.In Fig. 14(a) and Fig. 16(a), presenting the sweeps with the forward sail and the aft sail, it is evident that the sail-sail interaction has a strong effect on the location of the maximum   for these cases, shifting it by several degrees.
The yaw moment in the upwind cases is strongly affected by the interaction effect in all three sweeps, with the overall load center shifting forward because the lift force generated by the foremost sail increases and the force on the aftmost sail decreases.As evident from Fig. 15, even when the change in the total driving force due to the interaction is relatively small, the effect on the yaw moment in upwind conditions can be large.This effect is captured both by ISILL and ISILL+BL.However, without the boundary layer correction, the loads (a) Driving force coefficient.
(b) Yaw moment coefficient.on the foremost sail are underestimated, thereby resulting in a toolow yaw moment in the upwind cases.Studying the moments for   = 30 • in Fig. 14(b), it is evident that the difference between the results without considering interaction effects (SILL) and the results from CFD are large.This shift in the center of longitudinal effort has important implications for the ship's maneuverability.In addition, from Fig. 14(b), it is apparent that sheeting in the foremost sail, increasing  1 and the load on Sail 1, increases the absolute yaw moment until the sail stalls.The opposite effect on the yaw moment is seen in Fig. 16(b) when sheeting in the aftmost sail.The ISILL+BL method provides very good prediction of this behavior.
The results of the 2D sheeting sweeps are summarized in Tables 3-5.The relative errors in maximum driving force (  ) and maximum moment (  ) have been calculated for each validation case.Moreover, the deviation in the angle of attack (  ) for the sweeping sail at which the maximum driving force is predicted to occur is also presented.Consistent with the discussion in Section 4.1, the largest differences between ISILL and ISILL+BL are seen in the upwind cases, particularly in validation case 2.3, which is detailed in Fig. 13.
Study of pre-stall interaction effects in 3D Figs. 17-19 present the driving and side force coefficients for validation cases 3.1, 3.2, and 3.3, defined in Table 1.In each of these validation cases, all four sails are sheeted at  = 15 • to highlight the effect of the sail-sail interaction on each sail.Since the sails are sheeted equally, all four sails generate an equal force when the sailsail interaction is not considered.When the interaction is considered, there is generally a large variation in forces between the sails.
Fig. 17 (validation case 3.1) presents the results when the ship is sailing upwind at   = 30 • .Consistent with the 2D results, this upwind case shows strong sail-sail interaction effects.Fig. 17(a) depicts how the side force changes drastically from fore (Sail 1) to aft (Sail 4) due to the interaction.The ability to predict a correct moment has great implications when it comes to sailing in practice.Both versions of ISILL predict this change.There is little difference between the two versions because the foremost sails are not sufficiently close to their stalling points for the correction to have any effect.The negative effect on the downwind sails (Sails 2-4) is somewhat underestimated, thereby resulting in a small overprediction of the forces.
Fig. 18 (validation case 3.2) presents the results for an apparent wind angle of   = 60 • , corresponding to close reach.The forces are well predicted on all sails for this validation case.Fig. 19 (validation case 3.3) presents the force coefficients when the ship is sailing at a beam reach.Here, the forces on the upwind sails (Sail 1-3) are slightly overestimated, particularly the side force seen in Fig. 19(a).In this case, the SILL model without interaction actually corresponds better with the CFD results, but the difference is small.
As expected, it is evident that the total driving force decreases and the side force increases as the ship points closer to the wind.Moreover, in each validation case, the total driving force predicted by ISILL is lowered due to the interaction.One could draw the conclusion that the interaction effect is always unfavorable in terms of generated (a) Driving force coefficient.
(b) Side force coefficient.driving force.However, proper optimization of the sheeting angles could increase the driving force.Table 6 summarizes the 3D results by presenting the errors in total driving force (  ) and longitudinal center of effort (  ).To calculate the percentage error for the location of the LCE, the difference in meters has been normalized by dividing it by the distance between the rotation centers of Sails 1 and 4. As validation cases 3.1-3.3all have pre-stall effective angles of attack, there is little difference between ISILL with and without the boundary layer correction.

A sheeting sweep up to the point of stall in 3D
This section presents a sheeting sweep with Sail 1, corresponding to validation case 3.4 in Table 1.Here, the sheeting angles of Sail 2-4 have been adjusted to more realistic settings for an upwind case (  = 30 • ).In Fig. 20, the lift and drag force coefficients of each sail are plotted as a function of the angle of attack of Sail 1 ( 1 ).The angle is increased from 15 • to 20 • in order to find the point of stall of Sail 1.
In the CFD simulations, the maximum lift force is generated at  1 = 19 • .To illustrate the point of stall, the first point post-stall is also included, even though it has a high uncertainty since the RANS simulation no longer converges properly.This occurs as the forces fluctuate, and a steady-state simulation can no longer approximate the flow.The average of the forces predicted using RANS CFD are plotted in gray to indicate that the sail has stalled.The fluctuations are illustrated using an error bar based on the maximum and minimum values.In Fig. 20(a), it is evident that the difference between the lift force predicted by ISILL and ISLL+BL increases with the angle of attack.The method with the boundary layer correction predicts the same point of stall as the CFD simulations.Without the correction, ISILL predicts the point of stall at 17 • .Both versions improve the predicted maximum lift force.As the flow separates from Sail 1 the load on Sail 2 increases, evident from Figure 20(b).This effect is stronger in the CFD predictions than in the potential flow models, as the lift force decreases more drastically on Sail 1 in the CFD results.A small effect is also noted for Sails 2 (Fig. 20(c)) and 3 (Fig. 20(d)).When the flow separates on the forward sail, the angle of attack increases on the downstream sails.
Fig. 21 presents the total driving force and moment coefficient of each sheeting combination.The total driving force, presented in Fig. 21(a), is well predicted by ISILL with or without correction.In fact, the results without correction are slightly better.The largest improvements compared to SILL are presented in Fig. 21(b), which displays the total moment.Here, the advantage of the boundary layer correction is also very clear.
In the model, the settings of the boundary layer correction (   = 45% and    = 0.0025) were based only on 2D CFD simulations.

Study limitations and method generalizability
The main focus of this study was to evaluate the new sail-sail interaction model.The evaluation of SILL was of secondary interest since it had already been studied by Persson et al. (2019).It is possible to combine the suggested interaction model with alternative representations of the single sail forces and lift distribution, but that is not included in this study.
The validation was limited to wing sails standing in a row on a flat plate with a uniform incoming flow.The sails stand close together at a fixed distance, and the interaction effects are strong.The effect of anonuniform inflow, due to the presence of an atmospheric boundary layer or the hull and superstructure of the ship, is not considered.In cases where the sails are placed close to the ship's side, or the gap between the deck and sail is large, the assumption of a sail-force mirror plane in the deck should be challenged.
ISILL is only validated for wing sails, but since the sail-sail interaction model is based on circulation strength, the underlying theory should also apply to other lifting devices.However, due to the lifting line representation of the sail, the applicability of the model needs to be investigated in cases in which the distance between the sails is decreased.
Further, the tested apparent wind angles were limited to a range from upwind sailing to beam reach.These angles were selected since headwinds are predominant on the assumed Northern Atlantic route and when the ship is motor sailing, combining wind and propeller propulsion.Downwind sailing and apparent wind angles encountered when maneuvering are not considered.In fact, in its current form, the sail-sail interaction model does not consider the effect of the wake directly striking a downwind sail.If so, the viscous wake and the viscous damping at the center of the free vortices should be included   in the model.These effects must also be considered in a beam reach case if the sails are positioned side by side.For example, Bordogna has suggested a correction of the velocities in the wake region, which could be a suitable addition to the method presented here.In this study, ISILL was validated against RANS CFD simulations.As revealed in Section 2.3.1, pre-stall angles and the point of stall were well predicted when comparing the CFD simulations to wind tunnel tests in 2D.This is the relevant region for normal sailing and sheeting optimization.The uncertainty levels in the CFD predictions increase around and above stall.However, it can be noted that CFD simulations with similar settings are used to generate the validation data and the 2D   () and   () input to ISILL.This implies that the validation study, regardless of the accuracy of the absolute values, compares the ability of the ISILL method to model the 3D sail-sail interaction to that of a viscous high-fidelity method-that is, RANS CFD.

Conclusions
This paper presented a rapid method for predicting the performance of interacting rigid wing sails in 3D.The lifting line-based one-sail SILL method (Persson et al., 2019) was extended to multiple sails by introducing a potential flow-based interaction model.Two versions of the new method, ISILL and ISILL+BL, were evaluated, where ISILL+BL considers additional viscous interaction effects on the boundary layer and separation.The ability of the method to predict aerodynamic force and moment coefficients has been validated against RANS CFD simulations.
In the 2D validation study, sheeting sweeps at varying apparent wind angles with three wing sail profiles were carried out.They confirmed that the maximum driving force, the yaw moment, and the point of stall are strongly affected by the sail-sail interaction.The interaction model considerably improved the force predictions compared to SILL (not considering any interaction effects).In a few cases, ISILL underestimated the peak force, as the onset of stall was predicted at an angle of attack that was too low.Adding the boundary layer correction, ISILL+BL predicted the stall angle in all validation cases within an accuracy of 2 • , whereas SILL was off by up to 7 • .The largest benefit of introducing the boundary layer correction was seen in the yaw moment predictions.However, in certain upwind cases, an error in the yaw moment remained in the order of 10%.This indicates that even though ISILL is suitable for predicting the ship's speed and sail sheeting angles, complementary simulations with higher fidelity tools are required to control the final yaw balance of the ship.
The validation of the method in 3D consisted of two parts: (1) a prestall study of the interaction effects when sailing at different apparent wind angles and (2) a single sheeting sweep study with the foremost sail in an upwind condition.The pre-stall study revealed that both versions of ISILL predicted the effect of the sail-sail interaction on the generated sail forces well.The interaction affected the sail-generated forces for all tested apparent wind angles, but the effect was strongest for the upwind validation case.For the pre-stall sheeting angles, there was no difference between ISILL and ISILL+BL.In the 3D sheeting sweep study, the foremost sail reached its point of stall.Consistent with the 2D sweep studies, the boundary layer correction improved the predicted loads close to the stall angle.
The presented results suggests that combining viscous 2D profile data with a 3D potential flow-based interaction model is a suitable approach for rapid prediction of the driving force and the corresponding approximate optimal sheeting angles for a WPS with wing sails.The ISILL code, in its current version that was written without focus on efficiency, runs in seconds on a standard computer.The total computational cost, including the cost of generating 2D CFD input data and calibrating the boundary layer correction, is negligible compared to the cost of multiple 3D CFD simulations.The method fulfills the criteria of being sufficiently rapid to perform large numbers of simulations whilst maintaining acceptable force and stall angle predictions.
In developing the interaction model, care was taken to ensure its applicability to general WPSs.Yet, as discussed in Section 4.5, ISILL

Fig. 2 .Fig. 3 .
Fig. 2. Definitions of the coordinate systems, force vectors, apparent wind angle   , angle of attack , and sheeting angle .The apparent wind  ∞ is indicated by an arrow.

Fig. 4 .
Fig. 4. Grid convergence for   at  = 15 • for the 2D simulations.Grid 3, with a relative step size of 1.90, was used in the calculations.The uncertainty indicated in the figure is for the finest grid.

Fig. 6 .
Fig. 6.Grid convergence for   at  = 13 • for the 3D simulations.Grid 3, with a relative step size of 1.28, was used in the calculations.The uncertainty indicated in the figure is for the finest grid.

Fig. 7 .
Fig. 7.An example (validation case 3.4) of the mesh for four sails in a plane at  = 75 m.

Fig. 9 .
Fig. 9. Representation of a sail in the interaction model.Each sail is modeled by several superimposed horseshoe vortices (red), mirrored in the deck plane at midspan.The elliptical distribution of the circulation strength  () (green) is illustrated on the left side of the sail.

Fig. 10 .Fig. 11 .
Fig. 10.Illustration of the assumed vortex system in ISILL.The affected strip  on sail  is illuminated.

Fig. 13 .
Fig.13.Lift coefficient for Sail 1, Sail 2, and Sail 3 at an apparent wind angle of   = 30 • when sweeping with Sail 3 (v.c.2.3).Sails 1 and 2 are sheeted close to their maximum lift angles, thereby highlighting the effect of the BL-correction.

Fig. 20 .
Fig. 20.Predictions of the lift and drag coefficients when increasing the angle of attack on Sail 1 from  1 = 15 • to  1 = 20 • at   = 30 • (v.c.3.4).The other sails are sheeted to  2 = 12 • ,  3 = 10 • , and  4 = 8 • .To illustrate the point of stall, the results of the simulations above stall have been included even though they have a high uncertainty.

Fig. 21 .
Fig. 21.Aggregated results combining the loads on all four sails for a sheeting sweep with Sail 1 at   = 30 • .

Table 1
Test matrix displaying the different validation cases (v.c.).Every case has a set apparent wind angle (  ) of the ship and defined angles of attack (  ) for every sail .The italicized numbers indicate that   is swept over a span of angles.

Table 2
Overview of the RANS CFD computations.The chord length  (23 m) is used as the reference length.
(),   ()) The sail's own vortex system is excluded.The 3D effects on sail  are handled in Step 1c.(c) The new   (, ) gives an updated local angle of attack,   (, ), and local wind speed, |  |, which are used as inputs in Step 1.

Table 3
2D validation results for the sheeting sweeps with Sail 1.

Table 4
2D validation results for the sheeting sweeps with Sail 2.

Table 5
2D validation results for the sheeting sweeps with Sail 3.