Improving the Downwind Sail Design Process by Means of a Novel FSI Approach

The process of designing a sail can be a challenging task because of the difficulties in predicting the real aerodynamic performance. This is especially true in the case of downwind sails, where the evaluation of the real shapes and aerodynamic forces can be very complex because of turbulent and detached flows and the high-deformable behavior of structures. Of course, numerical methods are very useful and reliable tools to investigate sail performances, and their use, also as a result of the exponential growth of computational resources at a very low cost, is spreading more and more, even in not highly competitive fields. This paper presents a new methodology to support sail designers in evaluating and optimizing downwind sail performance and manufacturing. A new weakly coupled fluid–structure interaction (FSI) procedure has been developed to study downwind sails. The proposed method is parametric and automated and allows for investigating multiple kinds of sails under different sailing conditions. The study of a gennaker of a small sailing yacht is presented as a case study. Based on the numerical results obtained, an analytical formulation for calculating the sail corner loads has been also proposed. The novel proposed methodology could represent a promising approach to allow for the widespread and effective use of numerical methods in the design and manufacturing of yacht sails.


Introduction
In recent years, the exponential growth of computational resources at a very low cost has allowed for the widespread use of numerical methods in the nautical field as well as at an industrial level, not only for scientific research or highly competitive purposes. Nevertheless, to date, advanced numerical methods are commonly used in the nautical field only for competitive applications, but it is conceivable that in the future these methods could also be used in other application fields. Recent racing yachts, as demonstrated by the last edition of the America's Cup, have achieved very high performances thanks to the massive improvements in yacht and sail design, materials, and fabrication. The design of the hulls, sails, and rigging of competitive racing yachts needs more and more detailed research and development that can be obtained by combining advanced computational resources and experimental studies. In particular, with regard to the design of sails, based on traditional and empirical manufacturing processes, nowadays, the best sail designers aim to use more new research and development tools [1].
To improve sail design and to better understand sail aerodynamics, numerical simulations and experimental testing are the most used and reliable methods [2][3][4][5]. These methods have led to a better comprehension of the complex phenomena and the physics underlying sailing and, consequently, through them, a remarkable increase of the performances of yacht sails has been obtained in recent years [6,7]. Experimental full-scale testing is the most accurate and reliable method, but, because of cost constraints, its use is only 2 of 19 recommended for a few exceptional cases [7,8] or to validate numerical methods [2,9,10]. Instead, numerical simulations are cheaper than experimental methods, but in order to obtain reliable and accurate results, they need long setup times by highly qualified and experienced users and they require high-performing computational resources.
Among the numerical methods, computational fluid dynamics (CFD) simulations are the most widely used in the design process of high-performance sailing yachts [8]. CFD analyses can be effectively used to design high-performance hulls [11] and to simulate both upwind and downwind sail configurations [12][13][14][15][16][17]. Usually, the flow in upwind sailing is mostly attached and the turbulence effect is quite limited; the flow regime around a downwind sail can be highly turbulent [7] and, consequently, simulating real sailing conditions can be a very difficult task. For this reason, in the past, many studies have concentrated on upwind sails [6], and in recent years, many researchers have focused on simulating the complex downwind sailing conditions [6,18,19]. The first numerical studies on downwind sails frequently treated the topic in a simplified numerical way, considering only the flow effects around a rigid structure [17]. However, the physics of downwind sails is by far more complex than the upwind ones, not only because of the highly turbulent flow, but also because these kinds of sails have an inherent unsteadiness, even in quite stable conditions [7]. This happens mainly because downwind sails are made of very light and flexible cloth, and are attached to the yacht's structure at three points, namely, the head, the clew, and the tack [20]. The shape of a downwind sail is formed by self-generated aerodynamic forces that are strongly affected by the sail shape itself [21]. Moreover, the shape can change remarkably when sailing, depending on the trim settings and wind conditions [7]. All of these aspects must be taken into account when simulating the real sailing conditions of these kinds of sails, and the sail shape must be accurately measured in the flying condition [10,21]. For this purpose, nowadays, several authors focus on the issue of predicting the flying shape a sail develops under the impact of flow forces. In this context, it is well known that the fluid-structure interaction (FSI) is a key feature to study these kinds of problems [11,16,22,23]. Today, sail designers use specific software to define the constructed (molded) sail shape based on their experience, but they do not have any certainty about the real flying shape [6]. The possibility of predicting the real flying shape of a downwind sail could allow designers to improve the sailing performances of yachts [5,7,24,25], for instance, by maximizing the driving forces. However, knowing the real flying shape could also allow for optimizing the mechanical characteristics or the manufacturing process of the sails, in order to improve their stability, to reduce their weight, and to optimize their mechanical behavior [18,22,26,27]. Many papers have focused on evaluating the real shape of sails only to predict sailing performances [20,22,25,28]. However, to the best of our knowledge, no relevant papers have tried to measure the flying shape and the loads of a sail, taking into account both the point of view of sailors, whose primary interest is the sail performance, and of the sailmakers, who are also responsible for suitably dimensioning and manufacturing the sail.
In this work, a new approach based on the perspective of both sailors and sailmakers is proposed. In particular, the authors propose calculating the flying shape and the loads on downwind sails, not only to estimate propulsive forces [2,10], but also to help sailmakers in dimensioning and defining the best arrangement of the sail reinforcements. For this purpose, a weakly coupled FSI procedure [16] is used to find the real flying shape and to estimate the loads on the head, the clew, and the tack of the sail when the sailing conditions change. On the basis of the numerical results obtained, an analytical formulation for the calculation of the forces on the sail is also proposed.

Materials and Methods
In this work, to determine the flying shape and to evaluate the loads on downwind sails, a numerical FSI approach was used. Starting from a previous work of the authors [16], a new weakly coupled method was developed and tested. The new procedure is parametric and fully automated, and it only needs some input data to calculate the flying shape and the loads on the sail. As it is known, to perform a numerical FSI analysis, two embedded problems, aerodynamic and structural, need to be solved together [29]. In this study, to solve the aerodynamic and structural problems, the commercial CFD solver Ansys CFX and the Ansys Static Structural solver were used, respectively. JavaScript and a pythonbased programming language were used to develop the algorithm that managed the whole procedure, including the data exchange between the workbench and spreadsheet, which was developed ad hoc.

FSI Analysis: Setup of CFD Environment
The proposed procedure was developed to solve the FSI problem of a downwind sail but, with the aim to simulate the flow as accurately as possible, and so a mainsail was also considered in the CFD simulations. As mentioned, the procedure is parametric, in order to study different sailing configurations in a very simple and fast way. The main CFD parameters that must be first set to run the FSI analysis are reported in Table 1. The downwind sail trim angle, ∆α g , is a parameter that has been introduced to simulate different sail trims. In particular, the sheet angle of the downwind sail, that is the angle between the tack-clew line and the centerline of the boat (α g in Figure 1), is defined as α g = α g_design + ∆α g ; where α g_design is the design sheet angle of the downwind sail. In this way, once the specific sail type and its design sheet angle have been chosen, different sail trims can be simulated by simply changing the value of parameter ∆α g .

Materials and Methods
In this work, to determine the flying shape and to evaluate the loads on downwind sails, a numerical FSI approach was used. Starting from a previous work of the authors [16], a new weakly coupled method was developed and tested. The new procedure is parametric and fully automated, and it only needs some input data to calculate the flying shape and the loads on the sail. As it is known, to perform a numerical FSI analysis, two embedded problems, aerodynamic and structural, need to be solved together [29]. In this study, to solve the aerodynamic and structural problems, the commercial CFD solver Ansys CFX and the Ansys Static Structural solver were used, respectively. JavaScript and a python-based programming language were used to develop the algorithm that managed the whole procedure, including the data exchange between the workbench and spreadsheet, which was developed ad hoc.

FSI Analysis: Setup of CFD Environment
The proposed procedure was developed to solve the FSI problem of a downwind sail but, with the aim to simulate the flow as accurately as possible, and so a mainsail was also considered in the CFD simulations. As mentioned, the procedure is parametric, in order to study different sailing configurations in a very simple and fast way. The main CFD parameters that must be first set to run the FSI analysis are reported in Table 1. The downwind sail trim angle, Δαg, is a parameter that has been introduced to simulate different sail trims. In particular, the sheet angle of the downwind sail, that is the angle between the tack-clew line and the centerline of the boat (αg in Figure 1), is defined as αg = αg_design + Δαg; where αg_design is the design sheet angle of the downwind sail. In this way, once the specific sail type and its design sheet angle have been chosen, different sail trims can be simulated by simply changing the value of parameter Δαg. The domain used for the CFD simulation was made parametric by setting the positions of the main boundary surfaces (inlet, outlet, and top) as a function of the maximum sail chord, Cm, as defined in Table 2. The domain used for the CFD simulation was made parametric by setting the positions of the main boundary surfaces (inlet, outlet, and top) as a function of the maximum sail chord, C m , as defined in Table 2. Following a largely used approach, the boundary conditions at the inlet surfaces ( Figure 2) were set following the formulation proposed by Richards and Hoxey [30]. Concerning the other boundary surfaces, the outlet was set as "opening condition" [16], the bottom surface was set as a rough wall, and the top surface was set as the inlet. The flow velocity at the top surface was calculated using the same formulation of the velocity profile proposed in [30]. The k-ω SST formulation [31] was chosen as the turbulence model.  Following a largely used approach, the boundary conditions at the inlet surfaces (Figure 2) were set following the formulation proposed by Richards and Hoxey [30]. Concerning the other boundary surfaces, the outlet was set as "opening condition" [16], the bottom surface was set as a rough wall, and the top surface was set as the inlet. The flow velocity at the top surface was calculated using the same formulation of the velocity profile proposed in [30]. The k-ω SST formulation [31] was chosen as the turbulence model. Regarding the shapes of the sails, the geometry of the mainsail was fixed while the downwind sail was updated at each simulation according to the deformed shape calculated using the FEM analysis. To optimize the computational costs, the domain was divided into an inner region (grey in Figure 2), meshed using a mixed hexa-tetrahedral mesh with prismatic layers on the surfaces of the sails, and an outer region, where hexahedral elements were used. Mesh refinements were applied to the cells adjacent and close to the sail surfaces. The height of the first prismatic element attached to the surfaces of the sails was set to achieve the condition of y + ≈ 1.
A convergence study was performed to quantify the influence of the grid resolution on the CFD results. The ratio between the distances of the grid nodes is equal to The relative step ratio of the i-th grid is ℎ = √ 3 , with N being the number of cells, resulting in ℎ = 1.00 (base), 1.26, 1.59, and 2.00, respectively. The convergence trend is approximated with the following equation [32]: The values of coefficients S0, C, and p are computed with the least-squares method. Figure 3 and Table 3 show the drive force computed for the different analyzed grids and the results of the uncertainty estimation. A monotonic convergence can be observed in Figure 3. The value of the estimated uncertainty, Ug, for h = 1 (corresponding to 1.5 M cells grid) is 2.1%; this value can be considered acceptable. Regarding the shapes of the sails, the geometry of the mainsail was fixed while the downwind sail was updated at each simulation according to the deformed shape calculated using the FEM analysis. To optimize the computational costs, the domain was divided into an inner region (grey in Figure 2), meshed using a mixed hexa-tetrahedral mesh with prismatic layers on the surfaces of the sails, and an outer region, where hexahedral elements were used. Mesh refinements were applied to the cells adjacent and close to the sail surfaces. The height of the first prismatic element attached to the surfaces of the sails was set to achieve the condition of y + ≈ 1.
A convergence study was performed to quantify the influence of the grid resolution on the CFD results. The ratio between the distances of the grid nodes is equal to 3 √ 2. Four similar grids with about 1.5 M (base), 0.75 M, 0.38 M, and 0.18 M cells were investigated.
The relative step ratio of the i-th grid is h i = 3 N base N i , with N being the number of cells, resulting in h = 1.00 (base), 1.26, 1.59, and 2.00, respectively. The convergence trend is approximated with the following equation [32]: The values of coefficients S 0 , C, and p are computed with the least-squares method. Figure 3 and Table 3 show the drive force computed for the different analyzed grids and the results of the uncertainty estimation. A monotonic convergence can be observed in Figure 3. The value of the estimated uncertainty, U g , for h = 1 (corresponding to 1.5 M cells grid) is 2.1%; this value can be considered acceptable. ci. Eng. 2021, 9,    Basing on the grid verification results and the good performance hardware available, we used the 1.5 M cells grid. As a result of the unavailability of experimental data, it is not possible, to date, to present a complete validation of the numerical model based on a comparison of specific CFD results with corresponding experimental data. However, in Section 4.0, the results obtained with the proposed method are compared with others in the literature. A good level of agreement between the numerical and the experimental results has been noted; this demonstrates that the proposed method is able to give sufficiently reliable results. Of course, a complete validation will be carried out in future work.
The meshed domain and the details of mesh around the sail are shown in Figure 4.  Basing on the grid verification results and the good performance hardware available, we used the 1.5 M cells grid. As a result of the unavailability of experimental data, it is not possible, to date, to present a complete validation of the numerical model based on a comparison of specific CFD results with corresponding experimental data. However, in Section 4, the results obtained with the proposed method are compared with others in the literature. A good level of agreement between the numerical and the experimental results has been noted; this demonstrates that the proposed method is able to give sufficiently reliable results. Of course, a complete validation will be carried out in future work.
The meshed domain and the details of mesh around the sail are shown in Figure 4.   Basing on the grid verification results and the good performance hardware available, we used the 1.5 M cells grid. As a result of the unavailability of experimental data, it is not possible, to date, to present a complete validation of the numerical model based on a comparison of specific CFD results with corresponding experimental data. However, in Section 4.0, the results obtained with the proposed method are compared with others in the literature. A good level of agreement between the numerical and the experimental results has been noted; this demonstrates that the proposed method is able to give sufficiently reliable results. Of course, a complete validation will be carried out in future work.
The meshed domain and the details of mesh around the sail are shown in Figure 4.

FSI Analysis: Setup of FEM Environment
As a result of the assumption of a rigid mainsail, FEM analysis was performed only on the downwind sail. The FEM model is parametric; the main parameters are the sail geometry, the Young's modulus, the Poisson ratio, and the sail thickness. The geometry of the sail was automatically meshed using eight-node shell elements (Shell 281). This kind of shell element has six degrees of freedom at each node, and it is characterized by a mathematical formulation suitable for simulating thin structures (like sails). Suitable mesh refinements were applied on the head, the clew, and the tack.
A convergence study was performed with five structural meshes, with the following average sizes for the shell elements: 200, 100, 5, 30, and 20 mm. The different analyzed meshes had about 1540, 5700, 9200, 12,870, and 18,574 elements, respectively. The result of the convergence study is shown in Figure 11, where the value of the maximum displacement is plotted as a function of the average element size.

FSI Analysis: Setup of FEM Environment
As a result of the assumption of a rigid mainsail, FEM analysis was performed only on the downwind sail. The FEM model is parametric; the main parameters are the sail geometry, the Young's modulus, the Poisson ratio, and the sail thickness. The geometry of the sail was automatically meshed using eight-node shell elements (Shell 281). This kind of shell element has six degrees of freedom at each node, and it is characterized by a mathematical formulation suitable for simulating thin structures (like sails). Suitable mesh refinements were applied on the head, the clew, and the tack.
A convergence study was performed with five structural meshes, with the following average sizes for the shell elements: 200, 100, 5, 30, and 20 mm. The different analyzed meshes had about 1540, 5700, 9200, 12,870, and 18,574 elements, respectively. The result of the convergence study is shown in Figure 11, where the value of the maximum displacement is plotted as a function of the average element size.
ar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 9 of 19 Figure 11. Plot of the maximum displacement vs. average element size.
From Figure 11, it can be observed that convergence is obtained (relative differences between two consecutive meshes lower than 5%) for values of the element size lower than 100 mm. Given the result of this convergence study, the mesh with an average size of 50 mm was chosen. This represents a good compromise between the quality of the results and the computing time. From Figure 11, it can be observed that convergence is obtained (relative differences between two consecutive meshes lower than 5%) for values of the element size lower than 100 mm. Given the result of this convergence study, the mesh with an average size of 50 mm was chosen. This represents a good compromise between the quality of the results and the computing time.
The chosen mesh, composed of about 9200 elements, is shown in Figure 12. Figure 11. Plot of the maximum displacement vs. average element size.
From Figure 11, it can be observed that convergence is obtained (relative differences between two consecutive meshes lower than 5%) for values of the element size lower than 100 mm. Given the result of this convergence study, the mesh with an average size of 50 mm was chosen. This represents a good compromise between the quality of the results and the computing time.
The chosen mesh, composed of about 9200 elements, is shown in Figure 12. To calculate the loads on the rig, the fluid-dynamic pressure distribution calculated at each step by the CFD module was applied over both sail surfaces. The following boundary conditions were applied on the head and the tack (corners "B" and "C", respectively, in Figure 12): the displacements along the x, y, and z directions were fully constrained and only the rotations along the x, y, and z axes were allowed. The clew (corner "A" in Figure  12) was joined by a spherical joint to a link element, represented by a black dashed line in Figure 12, which simulates a rope. The end of the rope connected to the sail can move freely along the x, y, and z directions; the other end, instead, can only rotate, but no displacement is allowed. A multipoint constraint approach [16] was used to distribute the To calculate the loads on the rig, the fluid-dynamic pressure distribution calculated at each step by the CFD module was applied over both sail surfaces. The following boundary conditions were applied on the head and the tack (corners "B" and "C", respectively, in Figure 12): the displacements along the x, y, and z directions were fully constrained and only the rotations along the x, y, and z axes were allowed. The clew (corner "A" in Figure 12) was joined by a spherical joint to a link element, represented by a black dashed line in Figure 12, which simulates a rope. The end of the rope connected to the sail can move freely along the x, y, and z directions; the other end, instead, can only rotate, but no displacement is allowed. A multipoint constraint approach [16] was used to distribute the constraints over many nodes of the head, tack, and clew regions so to avoid peaks of stress due to the application of a load in a single node.

FSI Procedure
The flow chart of the developed procedure is shown in Figure 13. To launch the procedure, in addition to the CFD and FEM parameters, the following input parameters must be preliminarily defined: the sail design shape, S d ; -the number of steps, N s , that the FSI analysis is subdivided into; -the displacement threshold of the current deformed shape, T d .

FSI Procedure
The flow chart of the developed procedure is shown in Figure 13. To launch the procedure, in addition to the CFD and FEM parameters, the following input parameters must be preliminarily defined: -the sail design shape, Sd; -the number of steps, Ns, that the FSI analysis is subdivided into; -the displacement threshold of the current deformed shape, Td. Sd represents the shape of the sail as defined by the designer. The geometry (shape) file can be loaded in different file formats (iges, sat, dwg, etc.). As "classic" weak coupled procedures could be very unstable because of the large deformations of the sail under the wind load, in the new proposed procedure, the wind velocity was gradually increased in different subsequent steps, until the design value was reached. In this way, the sail shape gradually changed as the wind speed increased and, consequently, no convergence problems of the structural simulations arose. The parameter Ns is used to define the preliminary number of steps and influences the wind velocity, Vs, imposed at each step. The value of the velocity at the i-th step, Vs, is calculated as follows: Moreover, to better simulate the progressive deformation of the sail, thus avoiding any abrupt change in geometry between two subsequent steps, a check on the maximum value of the displacement was also introduced at the end of each step. For this purpose, a S d represents the shape of the sail as defined by the designer. The geometry (shape) file can be loaded in different file formats (iges, sat, dwg, etc.). As "classic" weak coupled procedures could be very unstable because of the large deformations of the sail under the wind load, in the new proposed procedure, the wind velocity was gradually increased in different subsequent steps, until the design value was reached. In this way, the sail shape gradually changed as the wind speed increased and, consequently, no convergence problems of the structural simulations arose. The parameter N s is used to define the preliminary number of steps and influences the wind velocity, V s , imposed at each step. The value of the velocity at the i-th step, V s , is calculated as follows: Moreover, to better simulate the progressive deformation of the sail, thus avoiding any abrupt change in geometry between two subsequent steps, a check on the maximum value of the displacement was also introduced at the end of each step. For this purpose, a threshold parameter, T d , has been defined. T d represents the admissible maximum value of the displacement of the deformed shape measured at each step.
A simple user interface developed in MS Excel ( Figure 14) allows for setting the numerical parameters and launching the FSI analysis.
As soon as it is started, the implemented procedure works iteratively in the following way: the sail geometry, which during the first iteration is the designed one (S d ), subjected to V s < AWS, is analyzed through the CFD module. The calculated fluid-dynamic pressure distribution is used as the boundary condition for the subsequent structural simulation performed by the FEM module. When the FEM simulation is completed, the maximum value of the total displacement of the sail, D max , is calculated and compared with the threshold value, T d . If D max is higher than T d , it means that the sail is deformed excessively, so the current velocity is decreased and the procedure restarts from the CFD simulation. If D max is lower than T d , the deformed shape of the sail is extracted from the FEM solver and is used as the updated input geometry for a new CFD analysis. Before starting the new CFD simulation, the current wind velocity, V s , is increased by an amount that depends on the value of N s . When V s = AWS, the procedure stops and the flying shape of the sail is obtained. threshold parameter, Td, has been defined. Td represents the admissible maximum value of the displacement of the deformed shape measured at each step. A simple user interface developed in MS Excel ( Figure 14) allows for setting the numerical parameters and launching the FSI analysis. As soon as it is started, the implemented procedure works iteratively in the following way: the sail geometry, which during the first iteration is the designed one (Sd), subjected to Vs < AWS, is analyzed through the CFD module. The calculated fluid-dynamic pressure distribution is used as the boundary condition for the subsequent structural simulation performed by the FEM module. When the FEM simulation is completed, the maximum value of the total displacement of the sail, Dmax, is calculated and compared with the threshold value, Td. If Dmax is higher than Td, it means that the sail is deformed excessively, so the current velocity is decreased and the procedure restarts from the CFD simulation. If Dmax is lower than Td, the deformed shape of the sail is extracted from the FEM solver and is used as the updated input geometry for a new CFD analysis. Before starting the new CFD simulation, the current wind velocity, Vs, is increased by an amount that depends on the value of Ns. When Vs = AWS, the procedure stops and the flying shape of the sail is obtained.

Case Study
A common all-purpose gennaker typically used to equip small size sailing boats, like dinghies, was studied. The sail area was 14.2 m 2 , the luff edge length was 5.95 m, and the maximum chord (Cm) length was 3.20 m. A nylon typically used for this kind of applications, characterized by a Young modulus of 5500 MPa, was chosen for this sail. As mentioned in Section 1.0, the aim of this paper is not only to present a new FSI approach to measure the flying shape of a downwind sail, but also to test the proposed method to evaluate the structural loads. In particular, it was investigated how the loads on the clew, the head, and the tack vary depending of different sailing conditions and, on the basis of the numerical results obtained, an attempt was made to find a possible analytical formulation that would allow for estimating loads as the wind speed varied. Based on the typical sailing conditions of these types of sails, we chose to study the gennaker at different values of AWA ranging between 90° and 130°, with 5° increments. Moreover, for each analyzed AWA value, we decided to investigate different configurations of the gennaker by varying its trim angle (Δαg) between 0° and 20° (with increments of 5°), in order to find the best trim depending on the wind angle. This resulted in about 60 different configurations to simulate. For this reason, to limit the computational time, a reduced wind velocity was chosen (AWS = 3 m/s). This choice allowed for investigating all of the sailing configurations in a reduced time using a Dual Intel Xeon E5-2670 processor (12 cores) workstation with 128 GB RAM. The mainsail boom angle (αm) was set to be constant and equal to 35°. The values of all of the parameters are presented in Table 4.

Case Study
A common all-purpose gennaker typically used to equip small size sailing boats, like dinghies, was studied. The sail area was 14.2 m 2 , the luff edge length was 5.95 m, and the maximum chord (C m ) length was 3.20 m. A nylon typically used for this kind of applications, characterized by a Young modulus of 5500 MPa, was chosen for this sail. As mentioned in Section 1, the aim of this paper is not only to present a new FSI approach to measure the flying shape of a downwind sail, but also to test the proposed method to evaluate the structural loads. In particular, it was investigated how the loads on the clew, the head, and the tack vary depending of different sailing conditions and, on the basis of the numerical results obtained, an attempt was made to find a possible analytical formulation that would allow for estimating loads as the wind speed varied. Based on the typical sailing conditions of these types of sails, we chose to study the gennaker at different values of AWA ranging between 90 • and 130 • , with 5 • increments. Moreover, for each analyzed AWA value, we decided to investigate different configurations of the gennaker by varying its trim angle (∆α g ) between 0 • and 20 • (with increments of 5 • ), in order to find the best trim depending on the wind angle. This resulted in about 60 different configurations to simulate. For this reason, to limit the computational time, a reduced wind velocity was chosen (AWS = 3 m/s). This choice allowed for investigating all of the sailing configurations in a reduced time using a Dual Intel Xeon E5-2670 processor (12 cores) workstation with 128 GB RAM. The mainsail boom angle (α m ) was set to be constant and equal to 35 • . The values of all of the parameters are presented in Table 4.  Figures 15 and 16 show the flying shapes of the gennaker obtained for two different sailing configurations. In particular, the presented results are the ones obtained for AWA = 90 • (∆α g = 0 • ) and AWA = 110 • (∆α g = 15 • ), respectively. The design shapes are colored in green and the flying ones are blue. The maps of the 3D deviations of the flying shapes from the design ones are also presented in Figures 15 and 16. Different point of views, two from the leech and the luff sides and one frontal, are presented in order to better detect the differences between the real and the flying shape.

Flying Shapes
Downwind sail trim angle (Δαg) 0° ÷ 20° Figures 15 and 16 show the flying shapes of the gennaker obtained for two different sailing configurations. In particular, the presented results are the ones obtained for AWA = 90° (Δαg = 0°) and AWA = 110° (Δαg = 15°), respectively. The design shapes are colored in green and the flying ones are blue. The maps of the 3D deviations of the flying shapes from the design ones are also presented in Figures 15 and 16. Different point of views, two from the leech and the luff sides and one frontal, are presented in order to better detect the differences between the real and the flying shape.

Mainsail boom angle (αm)
35° Downwind sail trim angle (Δαg) 0° ÷ 20° Figures 15 and 16 show the flying shapes of the gennaker obtained for two different sailing configurations. In particular, the presented results are the ones obtained for AWA = 90° (Δαg = 0°) and AWA = 110° (Δαg = 15°), respectively. The design shapes are colored in green and the flying ones are blue. The maps of the 3D deviations of the flying shapes from the design ones are also presented in Figures 15 and 16. Different point of views, two from the leech and the luff sides and one frontal, are presented in order to better detect the differences between the real and the flying shape.   It can be noticed that the largest positive displacements occur for AWA = 110 • and mainly occur along the leech edge. Moreover, it can be observed that for AWA = 90 • , the largest positive displacements mainly involve the upper part, while for AWA = 110 • , the lower part of the gennaker is the most deformed.

Loads on the Corners of the Sail: An Analytical Formulation
For all of the analyzed configurations, the driving force and the loads on the corners of the sail, which are the head, the tack, and the clew, were calculated during all of the steps of the FSI simulation. In this way, it was possible to know how these forces varied as the wind speed increased until the chosen value of AWS. To better investigate the structural loads on the head, the tack, and the clew of the sail, for each of the analyzed configurations, the polynomial trendlines of the corner loads and their analytical formulations were extracted by a specific tool of MS Excel. For example, Figure 17 shows the plots of the corner loads over the wind speed for AWA = 120 • and ∆α g = 20 • . The polynomial trendlines and their analytical formulations are also reported.

Loads on the Corners of the Sail: An Analytical Formulation
For all of the analyzed configurations, the driving force and the loads on the corners of the sail, which are the head, the tack, and the clew, were calculated during all of the steps of the FSI simulation. In this way, it was possible to know how these forces varied as the wind speed increased until the chosen value of AWS. To better investigate the structural loads on the head, the tack, and the clew of the sail, for each of the analyzed configurations, the polynomial trendlines of the corner loads and their analytical formulations were extracted by a specific tool of MS Excel. For example, Figure 17 shows the plots of the corner loads over the wind speed for AWA = 120° and Δαg = 20°. The polynomial trendlines and their analytical formulations are also reported. As soon as the analytical formulations of the trendlines were calculated for all of the configurations, a careful analysis of the values of the first and second order coefficients was performed. From this analysis, it emerged that, in all of the analyzed configurations, the first order coefficient was much smaller than the second order one. In particular, it was calculated that the first order coefficient was, on average, 3.7% of the second order coefficient. By taking this into account and in order to find a simplified analytical formulation for the calculation of the loads on the corners, we neglected the first order coefficient. For this reason, the following equation has been preliminarily proposed for the calculation of the load F on the corners: where CC [kg/m] is the load coefficient that has to be defined.
To determine the most suitable value of CC to make the Equation (3) usable in different sailing conditions, the second order coefficients were analyzed for all of the configurations. Figure 18 shows the values of the second order coefficients calculated for the clew, the tack, and the head in all of the analyzed configurations. As soon as the analytical formulations of the trendlines were calculated for all of the configurations, a careful analysis of the values of the first and second order coefficients was performed. From this analysis, it emerged that, in all of the analyzed configurations, the first order coefficient was much smaller than the second order one. In particular, it was calculated that the first order coefficient was, on average, 3.7% of the second order coefficient. By taking this into account and in order to find a simplified analytical formulation for the calculation of the loads on the corners, we neglected the first order coefficient. For this reason, the following equation has been preliminarily proposed for the calculation of the load F on the corners: where C C [kg/m] is the load coefficient that has to be defined.
To determine the most suitable value of C C to make the Equation (3) usable in different sailing conditions, the second order coefficients were analyzed for all of the configurations. Figure 18 shows the values of the second order coefficients calculated for the clew, the tack, and the head in all of the analyzed configurations. It can be seen that for each corner (clew, tack, and head), the coefficient values change as AWA and Δαg vary. For example, it can be noted that for the clew corner at AWA = 100° (red circle in Figure 18), the highest value of the coefficient was calculated when Δαg = 0°, while for the head corner at AWA = 120° (blue circle in Figure 18), the highest value was obtained when Δαg = 20°. To simplify the data analysis and to determine the most suitable values of CC to be used in Equation (2), we considered only the maximum values It can be seen that for each corner (clew, tack, and head), the coefficient values change as AWA and ∆α g vary. For example, it can be noted that for the clew corner at AWA = 100 • (red circle in Figure 18), the highest value of the coefficient was calculated when ∆α g = 0 • , while for the head corner at AWA = 120 • (blue circle in Figure 18), the highest value was obtained when ∆α g = 20 • . To simplify the data analysis and to determine the most suitable values of C C to be used in Equation (2), we considered only the maximum values of the second order coefficients at different values of AWA. Figure 19 shows the plots of the maximum values of the second order coefficients for the clew, the head, and the tack at different values of AWA. It can be seen that for each corner (clew, tack, and head), the coefficient values change as AWA and Δαg vary. For example, it can be noted that for the clew corner at AWA = 100° (red circle in Figure 18), the highest value of the coefficient was calculated when Δαg = 0°, while for the head corner at AWA = 120° (blue circle in Figure 18), the highest value was obtained when Δαg = 20°. To simplify the data analysis and to determine the most suitable values of CC to be used in Equation (2), we considered only the maximum values of the second order coefficients at different values of AWA. Figure 19 shows the plots of the maximum values of the second order coefficients for the clew, the head, and the tack at different values of AWA. Figure 19. Plots of maximum value of second order coefficients vs. AWA.
From Figure 19, it can be observed that the trends of the maximum values of the second order coefficients are substantially constant up to the value of AWA≈110°; beyond this value, a constant decrease can be. Based on this data, adopting a conservative approach, it was decided to use the average value of the second order coefficients calculated between AWA = 90° and AWA = 110° as the load coefficient CC for Equation (3). Table 5 reports the values of the second order coefficients and their average values calculated for the clew, the head, and the tack. From Figure 19, it can be observed that the trends of the maximum values of the second order coefficients are substantially constant up to the value of AWA≈110 • ; beyond this value, a constant decrease can be. Based on this data, adopting a conservative approach, it was decided to use the average value of the second order coefficients calculated between AWA = 90 • and AWA = 110 • as the load coefficient C C for Equation (3). Table 5 reports the values of the second order coefficients and their average values calculated for the clew, the head, and the tack. Following the proposed approach, the loads on the head, the tack, and the clew could be estimated using the following analytical formulations: F tack = C C_tack ·AWS 2 = 4.51·AWS 2 (5) F head = C C_head ·AWS 2 = 5.51·AWS 2 The values of the loads evaluated through CFD simulations and by Equations (4)-(6) are plotted in Figure 20.
The values of the loads evaluated through CFD simulations and by Equations (4)-(6) are plotted in Figure 20. It can be seen that for AWA < 110°, the values of the loads evaluated through the proposed analytical formulations are quite similar to the corresponding ones calculated using the CFD simulations. For values of AWA > 110°, as expected, the proposed equations returned slightly overestimated values for the loads, but, considering that a correct design approach should take into account the highest values of the loads on a structure, this could be acceptable.

Discussion and Conclusions
The study of the fluid-structure interaction is still one of the most interesting topics in order to predict the performances of sails. A sail, in fact, is usually made of thin plastic material or fabric-based composite, and it forms a three-dimensional camber shape with It can be seen that for AWA < 110 • , the values of the loads evaluated through the proposed analytical formulations are quite similar to the corresponding ones calculated using the CFD simulations. For values of AWA > 110 • , as expected, the proposed equations returned slightly overestimated values for the loads, but, considering that a correct design approach should take into account the highest values of the loads on a structure, this could be acceptable.

Discussion and Conclusions
The study of the fluid-structure interaction is still one of the most interesting topics in order to predict the performances of sails. A sail, in fact, is usually made of thin plastic material or fabric-based composite, and it forms a three-dimensional camber shape with wind pressure. The shape of the sail camber can remarkably change as the wind conditions vary, and this could result in a varied lift and drag performance by the sail [26]. These phenomena are relevant above all in the case of downwind sails, which are only attached to the yacht's structure at three points-the head, the clew, and the tack.
In this work, the issue of the FSI problem of downwind sails has been addressed by developing and testing a new weakly coupled procedure. Regarding the presented procedure, one of the main innovative features is related to the possibility of simulating the progressive deformation of the sail, through the gradual application of aerodynamic loads. This feature allows for reducing convergence problems and to avoid solution instabilities, in order to overcome the typical drawbacks of classic weak coupled procedures [19,33].
Another key feature of the new procedure is that it is parametric and fully automated. In fact, both CFD and FEM modules have been setup in such a way as to automatically update the models and the boundary conditions when simulation parameters change. Of course, the procedure is parametric within certain limits; if the type of sail to be analyzed changes remarkably, it is necessary to manually create the preliminary setup of the numerical models. Concerning the automation of the procedure, thanks to a specifically developed algorithm that fully manages the process and the data exchange between the CFD and the FEM modules, different sailing conditions can be simulated by simply introducing some parameters in a MS Excel user-interface.
The new procedure has been used to study a gennaker of a small sailing boat to evaluate its flying shape and to calculate the forces on the corners of the sail. The FSI study has been completed with no interruption or convergence problems within the considered sailing conditions. This demonstrates that, for the specific analyzed test case, the new procedure is effective and stable.
An extensive literature search was carried out to find experimental data to validate the numerical proposed methodology. Unfortunately, to the best of our knowledge, in the literature, there are no experimental studies on sails identical or very similar to the one studied here. However, some experimental studies performed on downwind sails show that the results obtained with the proposed methodology can be considered sufficiently reliable. As evidence of this, the values of the drive force coefficient, CFx, have been calculated at different AWAs using the numerical method proposed here ( Figure 21) and compared with the values measured experimentally by Motta et al. [34]. Some similarities can be evidenced.
be analyzed changes remarkably, it is necessary to manually create the preliminary setup of the numerical models. Concerning the automation of the procedure, thanks to a specifically developed algorithm that fully manages the process and the data exchange between the CFD and the FEM modules, different sailing conditions can be simulated by simply introducing some parameters in a MS Excel user-interface.
The new procedure has been used to study a gennaker of a small sailing boat to evaluate its flying shape and to calculate the forces on the corners of the sail. The FSI study has been completed with no interruption or convergence problems within the considered sailing conditions. This demonstrates that, for the specific analyzed test case, the new procedure is effective and stable.
An extensive literature search was carried out to find experimental data to validate the numerical proposed methodology. Unfortunately, to the best of our knowledge, in the literature, there are no experimental studies on sails identical or very similar to the one studied here. However, some experimental studies performed on downwind sails show that the results obtained with the proposed methodology can be considered sufficiently reliable. As evidence of this, the values of the drive force coefficient, CFx, have been calculated at different AWAs using the numerical method proposed here ( Figure 21) and compared with the values measured experimentally by Motta et al. [34]. Some similarities can be evidenced. The trends of the coefficient CF X as the AWA varies are very similar; in fact, they are almost constant with a slight increase around 100 • -110 • . Moreover, the order of magnitude of the coefficients is comparable, with 0.73 being the average value calculated with the numerical methodology proposed here and about 0.58 being the average value measured experimentally. Of course, as the compared sails are different, it is not possible to have identical or very similar values of CF X , but considering that the order of magnitude is identical, it could be deduced that the numerical results are sufficiently reliable. With reference to the loads calculated numerically with the proposed methodology, good levels of agreement were found with other case studies in the literature. In particular, it can be observed that the ratios of the loads on the sail corners and the trends of the forces, when AWA varies, are very similar to those found in other studies [2,20,34]. For example, the force ratios F clew /F head and F tack /F head calculated at AWA = 90 • using the new procedure are about 48% and 81%, respectively; the same ratios have been experimentally measured by Deparday et al. [2] for a similar case study and are equal to about 51% and 85%, respectively. Regarding the trends of the loads, it can be observed that, similarly to what was found with the numerical procedure proposed here, as well as in other experimental studies on similar kinds of sails [2,34], a decrease in the forces was detected when AWA > 110 • . This good level of agreement between the numerical data presented and other experimental data, even if it cannot be considered a direct validation of the numerical results, demonstrates the proposed method could give enough reliable results and can be used as a useful tool for sail design.
With regard to the analytical formulation for calculating the sail corner loads discussed in Section 3.2, the proposed Equation (3) of course cannot be used to calculate the corner loads for all types of sails. It applies only for the specific kind of gennaker; however, in our opinion, the proposed approach could be interesting because it could be effectively repeated to find other formulations for different sail types. Still, nowadays, in fact, most sailmakers design the reinforcements to apply on the head, the tack, and the clew, based on their experience. Analytical formulations, grouped by similar types of sails, could be an inexpensive and easy to use method to suitably dimension the reinforcements on the sail corners. Moreover, knowledge of the loads of the sail corners could allow for a better comprehension on wind/rig/sail interaction [2,10]. Such an approach could allow for remarkably improving the design and manufacturing process of sails, not just for competitive sailing yachts, but also for cruising ones.
For future work, the authors believe that further developments could be addressed to the setup and validation of a database of analytical formulations in order to evaluate the loads on the sail corners and to the improvements of the numerical models. A database of equations for an initial evaluation of the sail loads could be developed by analyzing numerous sails different in shape and size; in this way, sailmakers could use a simple tool to design the most suitable reinforcements for the sails, based not only on their experience, but also on numerical data. CFD simulations could also be improved by investigating unsteady sailing conditions, while FEM models could be enhanced by introducing additional information on different types of sail materials regarding structure non linearities, and by evaluating the rig deformation. In this way, very accurate and complete information could be obtained and the sail performance could be optimized. In fact, as the developed procedure is based on a weak coupled FSI method, the FEM module could be effectively interfaced with an optimization module to find the best structural parameters for the sails (e.g., optimal thickness of different regions of the sail and most suitable layout of composite materials).
In conclusion, in our opinion, the novel proposed methodology is surely promising, and further enhancements could allow for getting more information on the dynamics of downwind sails and could lead to the widespread use of this procedure in the design and manufacturing of yacht sails.  Acknowledgments: The authors gratefully acknowledge ANSYS Inc. for the support and the issue of academic licenses.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: