# A fast and robust approach to compute nonlinear Fluid-Structure Interactions on yacht sails

Matthieu Sacher, Jean-Baptiste Leroux, Alain Nême, Christian Jochum

► To cite this version:

Matthieu Sacher, Jean-Baptiste Leroux, Alain Nême, Christian Jochum. A fast and robust approach to compute nonlinear Fluid-Structure Interactions on yacht sails: Application to a semi-rigid composite mainsail. 24e Congrès Français de Mécanique, Aug 2019, Brest, France. hal-03271369

HAL Id: hal-03271369

https://hal.science/hal-03271369

Submitted on 25 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A fast and robust approach to compute nonlinear Fluid-Structure Interactions on yacht sails – Application to a semi-rigid composite mainsail

M. SACHERª, J. B. LEROUXᵇ, A. NÊMEᶜ, C. JOCHUMᵈ

a.–d. ENSTA Bretagne, CNRS UMR 6027, IRDL, 2 rue Francois Verny, 29806 Brest Cedex 9, France

a. matthieu.sacher@ensta-bretagne.fr

b. jean-baptiste.leroux@ensta-bretagne.fr

c. alain.neme@ensta-bretagne.fr

d. christian.jochum@ensta-bretagne.fr

## Abstract

A numerical method to compute the nonlinear Fluid-Structure Interactions (FSI) of yacht sails is presented in this paper. The inviscid lifting-line flow model is used by including a quadrature method to efficiently compute sail pressure loads. The structural calculation is performed with a quasi-static resolution by using a dynamic backward Euler scheme, in order to improve the computation convergence. A specific thickness sail approach is also proposed to make the FSI solving easier. The assembly of these flow and structural methods leads to a fast and robust strategy to compute nonlinear FSI on yacht sails, and the proposed approach is applied on a complex semi-rigid composite mainsail.

**Keywords: Fluid, structure, numerical, fast, robust, sail, membrane.**

## 1 Introduction

In most cases, yacht sails are made of thin fabrics that can be subjected to large structural displacements and deformations due to the complex wind aerodynamic load pressure [30]. The flying shape of the sail depends on the wind, the mast and sail geometries and structural properties, in addition to the rig trimming. Similarly, the wind pressure field depends on the sail shape, and it is therefore necessary to consider a nonlinear coupled problem to correctly model the sails Fluid-Structure Interactions (FSI) physics.

FSI modeling of yacht sails has now been studying for a few years [5, 8, 28]. The numerical state of the art shows that FSI solvers have been compared and validated with experiments [2, 3], and they are now also applied to several optimization problems [22, 25]. However, due to the very complex aero-elastic phenomena of yacht sails, such nonlinear FSI computations remain challenging. The modeling issues can concern the fluid and structural solvers, in addition to the coupling strategy. Structural solvers are usually based on finite element methods, while for flow solvers, inviscid models [10] are widely used, because of their fast computing time in comparison to the viscous flow ones [12]. In many cases these
computations are still very expensive, in addition to convergence problems sometimes. As a result, they are used little or not at all in industrial projects, and there is a real need today, for some companies, to have fast and robust FSI yacht sails models, for precise design processes, in a moderate engineering time.

Following the FSI model concern, the present paper focuses on an original fast and robust approach to compute FSI effects on a yacht mainsail. Specifically, the proposed strategy uses the lifting-line theory [21] combined with a semi-analytical Chebyshev-Gauss quadrature method [9, 19, 20], that allows efficient aerodynamic sail pressures computations. Regarding the structural part, the present method uses the Abaqus™ 2017 finite element software. A quasi-static resolution with a dynamic backward Euler scheme [13] is carried out. It facilitates the convergence of nonlinear static equilibrium owing to dynamic relaxation and the natural damping of the scheme. A standard Newton-Raphson nonlinear static resolution is then performed to ensure the exact static convergence. The sail membrane is modeled by shell elements and a specific thickness approach is employed in order to make the nonlinear sail FSI resolution solving [18] easier.

The present work is applied within the "Solid Sail 2.0" project, and enables sail FSI computing on a real industrial case of study. Indeed, this project is led by the French "Chantiers de l'Atlantique" company, who's purpose is to develop new designs of rigid sails. The main benefit of such sails concerns the life-cycle which is significantly increased and then makes it possible to consider this kind of sail design for large merchant ships. However, the industrial design of rigid sails is a complex problem that requires accurate numerical models to avoid large weights in the top rig parts.

A 300 m² mainsail prototype was built by the "Chantiers de l'Atlantique", and it has the specificity to be composed of 21 articulated composite panels (see Figure 1). Each panel is set with 2 battens and

<figure><img src="image_2.png" /><figcaption>Figure 1: Semi-rigid prototype mainsail hoisted on the "Ponant" yacht.</figcaption></figure>

it is linked to its neighbours through flexible loops at both luff and leech. Synchronization straps are also used in order to reduce possible relative transverse displacements. The sail is to be considered as "semi-rigid" since membranes flexibility's are close to a conventional soft sail, as membrane thickness at the central part of the panel is about 0.5 mm.

The following sections present the computation approach for the FSI problem of the prototype mainsail, and are organized as follows. Section 2 presents the numerical modeling methods by including the proposed FSI strategy. Numerical results are then reported in Section 3, and finally conclusions of this work and future development paths are provided in Section 4.
# 2 Numerical modeling

Numerical methods used to compute FSI on the prototype mainsail, as shown in Figure 1, are presented in this Section 2. In the current work, the study is limited to an upwind sailing case, mast deformations are neglected and the jib is not modeled. We also assume that the yacht mainsail FSI problem has a steady solution, which allows us to couple flow and structural solvers with a quasi-static approach. The flow and structural solvers are described below.

## 2.1 Inviscid flow solver

The flow solver relies on the well known lifting-line theory [21] and is firstly summarized in the following Section 2.1.1. Then, Section 2.1.2 describes the leading idea to use a semi-analytical method to efficiently compute the aerodynamic sail pressure loads.

### 2.1.1 Lifting-line theory

The lifting-line theory (LLT) was established by Prandtl [21], based on previous Lanchester [17] research work. This model allows to compute the spanwise distribution of section lift and induced drag of a three-dimensional finite wing. For an incompressible and inviscid flow, the wing (or sail) is modelled as a single vortex line, located at the quarter chord position, and an associated shed vortex sheet. Figure 2 illustrates this assumption and introduces related parameters: the distance y across the sail span measured from the boom, the sail half-span b, the circulation Γ(y), the reference flow velocity U₀(y), the induced flow velocity ωᵢ(y), the effective flow velocity Uₑ(y), the angle of attack α(y), the induced angle of attack αᵢ(y) and finally the effective angle of attack αₑ(y). As shown in Figure 2(b), the effective

<figure><img src="image_3.png" /><figcaption>Figure 2: Lifting-line model.</figcaption></figure>

flow velocity and angle of attack at each sail section is affected by the three dimensional shed vortices, also known as downwash effects. The Biot–Savart’s Law gives rise to the induced velocity definition with

$$ \omega_i(y) = -\frac{1}{4\pi} \int_{-b}^{b} \frac{d\Gamma(y_0)/dy}{y - y_0} dy_0, \quad -b < y < b, \qquad (1) $$
and by assuming low $\alpha(y)$ values, the linearization $\alpha_i \approx -\frac{\omega_i}{U_0}$, allows the effective angle of attack computation with

$$ \alpha_e(y) = \alpha(y) - \alpha_i(y), \qquad (2) $$

$$ = \alpha(y) + \frac{1}{4\pi U_0(y)} \int_{-b}^{b} \frac{d\Gamma(y_0)/dy}{y - y_0} dy_0. \qquad (3) $$

The lift spanwise distribution per unit length is computed by

$$ L'(y) = -\rho_0 U_e(y) \Gamma(y), \qquad (4) $$

$$ = \frac{1}{2} \rho_0 c(y) U_e^2(y) C_L(y), \qquad (5) $$

where $\rho_0$ is the reference fluid density and $C_L(y) = 2\pi \sin(\alpha_e(y) - \alpha_0(y))$ is the section lift coefficient. The fundamental Prandtl's integro-differential equation is then obtained from the linearization of the Kutta relation (4) and (5). It is given by

$$ \Gamma(y) = -K(y)c(y)U_0(y) \left( \alpha(y) - \alpha_0(y) + \frac{1}{4\pi U_0(y)} \int_{-b}^{b} \frac{d\Gamma(y_0)/dy}{y - y_0} dy_0 \right), \qquad (6) $$

where $K(y) \approx \pi$, and $\alpha_0(y)$ is the two-dimensional sail section zero-lift angle.

By neglecting the sweep and dihedral sail effects, an analytical solution of (6) can be obtained in terms of a Fourier sine series [11, 16]. For a given sail geometry with a known spanwise distribution of $U_0(y)$, $c(y)$, $\alpha(y)$, $\alpha_0(y)$, the circulation distribution is written as a truncated Fourier series where the Fourier coefficients are determined by forcing the lifting-line equation to be satisfied at N specific sail sections along the span. In these conditions and from the change of variables, $y = b \cos(\theta)$, the circulation distribution is then given by

$$ \Gamma(\theta) = -4bU_0(\theta) \sum_{n=1}^{N} A_n \sin(n\theta), \quad 0 < \theta < \pi, \qquad (7) $$

while the induced velocity is

$$ \omega_i(\theta) = -U_0(\theta) \sum_{n=1}^{N} n A_n \frac{\sin(n\theta)}{\sin(\theta)}. \qquad (8) $$

Finally, after substituting the circulation distribution and induced angle of attack in (6), the Fourier coefficient $A_n$ are determined from the relation

$$ \sum_{n=1}^{N} \sin(n\theta) \left( \frac{4b}{K(\theta)c(\theta)} + \frac{n}{\sin(\theta)} \right) A_n = \alpha(\theta) - \alpha_0(\theta). \qquad (9) $$

The zero-lift angle of attack, at each $\theta$ in (9), can be computed from the thin sections of the sail geometry. Indeed, Abbott and Von Doenhoff [1] suggest the following expression to compute the zero-lift angle of attack of a thin profile section,

$$ \alpha_0 = \int_0^1 \frac{\eta_c(x)}{c} \mu\left(\frac{x}{c}\right) \frac{dx}{c}, \quad 0 < \frac{x}{c} < 1, \qquad (10) $$
where $\eta_c(x)$ is the profile camber law and $\mu(\frac{x}{c})$ is given by

$$ \mu(\frac{x}{c}) = \frac{1}{\pi(1 - \frac{x}{c})\sqrt{\frac{x}{c}(1 - \frac{x}{c})}} $$

In the present work, the spanwise distribution of $\alpha_0$ is computed with (10), by using $\eta_c(x)$ interpolations, based on Chebyshev polynomials [19].

From (9), the lift and induced drag coefficients are obtained with

$$ C_L = \pi \lambda A_1, \qquad (11) $$

$$ C_{D_i} = \pi \lambda \sum_{n=1}^{N} n A_n^2, \qquad (12) $$

where $\lambda = \frac{b^2}{S}$ is the sail aspect ratio.

In this form, the lifting-line method does not allow to evaluate the aerodynamic sail pressure loads. Indeed, global loads, such as lift and induced drag, are computed by the LLT, but it is necessary to add a calculation step to get the aerodynamic pressures, that are locally applied to the sail. This step is presented in the next Section 2.1.2.

## 2.1.2 Semi-analytical method for the aerodynamic pressure field

As shown previously, the LLT method alone is not able to compute the local aerodynamic sail pressure loads. To do this, the present work uses small-disturbance flow over the two-dimensional sail sections, based on singularity element distributions. Specifically, a vortex distribution based model is implemented, as suggested for thin lifting airfoil by Katz and Plotkin [16].

A two-dimensional sail section is presently considered, with a camber law $\eta_c(x)$, a chord length $c$, and an effective angle of attack $\alpha_e$. The flow velocity is $U_e$ and vortex elements are distributed along the section chord axis. The unknown vortex distribution $\gamma(x)$ has to satisfy the zero normal flow boundary condition on the sail section. In other words, the flow remains tangent to the camberline, and leads to following equation

$$ \frac{1}{2\pi} \int_{0}^{c} \gamma(x_0) \frac{dx_0}{x - x_0} = U_e \left( \frac{d\eta_c(x)}{dx} - \alpha_e \right), \quad 0 < x < c, \qquad (13) $$

which is the integral equation for $\gamma(x)$ to be determined. The sail section pressure differences $\Delta_P$ can then be computed by the steady-state Bernoulli equation [16] with

$$ \Delta_P(x) = \rho_0 U_e \gamma(x). \qquad (14) $$

The computation of $\gamma(x)$ in (13) is now presented. Equation (13) is firstly rewritten in the form

$$ \int_{0}^{c} \frac{\Phi(x_0)}{x_0 - x} dx_0 = F(x), \qquad (15) $$
where,

$$ 
\begin{aligned}
\Phi(x_0) &= \frac{\gamma(x_0)}{2\pi}, \quad 0 < x_0 < c, \\
F(x) &= -U_e \left( \frac{d\eta_c(x)}{dx} - \alpha_e \right), \quad 0 < x < c.
\end{aligned}
 $$

From the following variables change,

$$ 
\begin{aligned}
t &= \frac{2x_0}{c} - 1, \quad -1 < t < 1, \\
\chi &= \frac{2x}{c} - 1, \quad -1 < \chi < 1,
\end{aligned}
 $$

and,

$$ \phi(t) = \Phi(x_0), \quad f(\chi) = F(x), $$

Equation (15) becomes

$$ \int_{-1}^{1} \frac{\phi(t)}{t-\chi} dt = f(\chi), \quad -1 < \chi < 1. \qquad (16) $$

This last Equation (16) is of a specific form, that is called a singular integral equation of the first kind [20]. If $\phi$ is unbounded in $\chi = -1$ and bounded in $\chi = 1$, to agree with the Kutta condition at the trailing edge, a complete analytical solution [20] of (16) is given by

$$ \phi(\chi) = -\frac{1}{\pi^2} \sqrt{\frac{1-\chi}{1+\chi}} \int_{-1}^{1} \sqrt{\frac{1+t}{1-t}} \frac{f(t)}{t-\chi} dt, \quad -1 < \chi < 1. \qquad (17) $$

However, the analytical calculation cannot be performed further, since the camber shape is not always given as a known function. It is then necessary to rely on quadrature methods, such as the Chebyshev-Gauss ones proposed by Eshkuvatov et al. [9]. Specifically, Chebyshev polynomials of the fourth kind are used to build an approximation $\phi_n$ of the unknown function $\phi$ in (16). This polynomial function is defined by

$$ \phi_n(\chi) = w(\chi) \sum_{i=0}^{n} \beta_i W_i(\chi), \qquad (18) $$

where $\beta_i$ are unknown coefficients. The weight function $w$ and the Chebyshev polynomials $W_i$ in (18) are defined by

$$ w(\chi) = \sqrt{\frac{1-\chi}{1+\chi}}, \quad W_i(\chi) = \frac{\sin\left((i - \frac{1}{2}) \cos^{-1}(\chi)\right)}{\sin\left(\frac{1}{2} \cos^{-1}(\chi)\right)}. $$

Substituting the approximate solution (18) in (16) gives

$$ \sum_{i=0}^{n} \beta_i \int_{-1}^{1} \sqrt{\frac{1-t}{1+t}} \frac{W_i(t)}{t-\chi} dt = f(\chi), \quad -1 < \chi < 1. \qquad (19) $$
From the following equation [19],

$$ \int_{-1}^{1} \sqrt{\frac{1-t}{1+t}} \frac{W_i(t)}{t-\chi} dt = -\pi V_i(\chi), \qquad (20) $$

where $V_i$ are the Chebyshev polynomials of the third kind given by

$$ V_i(\chi) = \frac{\cos\left((i - \frac{1}{2}) \cos^{-1}(\chi)\right)}{\cos\left(\frac{1}{2} \cos^{-1}(\chi)\right)}, $$

and by using the zeros of $V_{n+1}$ polynomials, given by

$$ \chi_k = \cos\left(\frac{(k - \frac{1}{2})\pi}{n + \frac{3}{2}}\right), \quad k = 1, \dots, n + 1, \qquad (21) $$

the $\beta_i$ coefficients are solved with

$$ \sum_{i=0}^{n} \beta_i V_i(\chi_k) = -\frac{f(\chi_k)}{\pi}, \quad k = 1, \dots, n + 1. \qquad (22) $$

The approximate solution of (16) is finally given by

$$ \phi(\chi) \approx w(\chi) \sum_{i=0}^{n} \beta_i W_i(\chi), \qquad (23) $$

that is used to compute sail pressure loads in (14) with

$$ \gamma(x) = 2\pi\phi\left(\frac{2x}{c} - 1\right), \quad 0 < x < c. \qquad (24) $$

## 2.2 Structural solver

The structural model is built within the Abaqus™ 2017 finite element analysis (FEA) software. This powerful commercial code is able to compute nonlinear FEA over a wide range of element types, using various integration schemes and including possible complex contacts under transient loads.

### 2.2.1 Numerical model

The numerical model is based on the mainsail previously presented in Figure 1. It is composed of 21 composite panels connected to each other thanks to spring elements (see Figure 3(a)). The sail battens are modeled with equivalent beam stringer elements and “STR13” thin shell Kirchhoff elements [6] are used for the modeling of composite panel membranes. Draping layers composed of carbon and glass fibers are applied on each panel built from the membrane cutting plan that can be observed in Figure 3(a). A rigid link is applied between the luff of the panels and associated references points that can rotate and move along the rigid mast rail. The mainsail head is hooked and the first bottom panel is linked to a rigid boom with 2 axial springs, connected to 2 references points, that have only rotations as degrees of freedom. Regarding the mainsail loads, the gravity is taken into account and the aerodynamic loads are applied by using uniform pressures on each shell element face.
<figure><img src="image_4.png" /><figcaption>Figure 3: Structural mainsail model.</figcaption></figure>

Structural computations of thin sail membranes are usually complex due to nonlinear behaviours. Indeed, sail deformations generate possible buckling related wrinkles [29], and make challenging problem-solving because of numerical instabilities. Wrinkles models have been developed [7, 14, 31] when using membrane elements, to account for local deformations with a reasonable number of elements. In this work, shell elements that are able to numerically reproduce the wrinkling development are used. However, because of the challenging problem-solving in presence of wrinkles, a specific shell thickness approach is implemented to enhance the computation stability. This strategy allows to correctly capture the sail deformations by neglecting the low structural effects due to small wrinkle scales.

From the plate theory and under the Kirchhoff hypothesis [23, 27], the plates are supposed to be thin and straight lines perpendicular to the mid-surface $(\vec{x}, \vec{y})$ before deformation, and to remain straight after deformation. By considering that the relations between stresses $\sigma$ and strains $\varepsilon$ are Hooke's law, the plane stress-reduced constitutive equations for such thin isotropic plates can be expressed as

$$ \begin{pmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{pmatrix} = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & 1 - \nu \end{bmatrix} \begin{pmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{xy} \end{pmatrix}, \qquad (25) $$

where $E$ and $\nu$ are the Young's modulus and the Poisson's ratio. The forces and moments per unit length, for a thin plate thickness $h$, are given by

$$ \begin{pmatrix} N_{xx} \\ N_{yy} \\ N_{xy} \end{pmatrix} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \begin{pmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{pmatrix} dz, \qquad \begin{pmatrix} M_{xx} \\ M_{yy} \\ M_{xy} \end{pmatrix} = - \int_{-\frac{h}{2}}^{\frac{h}{2}} \begin{pmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{pmatrix} z dz, \qquad (26) $$

and by using (25), the relations in (26) become

$$ \begin{pmatrix} N_{xx} \\ N_{yy} \\ N_{xy} \end{pmatrix} = C \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & 1 - \nu \end{bmatrix} \begin{pmatrix} \varepsilon_{p_{xx}} \\ \varepsilon_{p_{yy}} \\ \varepsilon_{p_{xy}} \end{pmatrix}, \qquad \begin{pmatrix} M_{xx} \\ M_{yy} \\ M_{xy} \end{pmatrix} = D \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & 1 - \nu \end{bmatrix} \begin{pmatrix} \chi_{xx} \\ \chi_{yy} \\ \chi_{xy} \end{pmatrix}, \qquad (27) $$
where $\epsilon_p$ is the plane strain-reduced of the mid-surface, $\chi$ is the curvature matrix of the mid-surface, and the extensional stiffness $C$ and bending stiffness $D$ are written as

$$C = \frac{Eh}{1 - \nu^2}, \qquad (28)$$

$$D = \frac{Eh^3}{12(1 - \nu^2)}. \qquad (29)$$

The dimensionless ratio $\frac{D}{CS_p}$, where $S_p$ is a reference surface of a thin plate, allows to quantify the ratio of bending stiffness to extension. For a fixed $S_p$ and by modifying the thickness $h$, the $\frac{D}{CS_p}$ ratio can vary. If it is very low, the shell mechanical comportment is getting similar to the membrane one and numerical instabilities may appear during the problem-solving. If $\frac{D}{CS_p}$ is large, numerical instabilities drop but the shell is no more close to membrane behaviour. In the present work, the sail panels are of elongated shapes between the luff and the leech. The reference surface is then replaced by $S_p = L_p^2$, where $L_p$ is the averaged luff to leech panel length. After several shell to membrane comparative tests (not shown in this paper), a good compromise of $\frac{D}{CL_p^2} \approx 10^{-7}$ was found and used to define the fictive shell thickness panels $h_f$ of the mainsail. Indeed, the fictive shell thickness, for a panel reference length $L_p$, is given by

$$h_f = L_p \sqrt{12 \cdot 10^{-7}}, \qquad (30)$$

and the effective shell mechanical properties $E$, $G$ and $\rho$ are replaced by fictive magnitudes as following, in order to keep constant the membrane properties,

$$E_f = \frac{Eh}{h_f}, \quad G_f = \frac{Gh}{h_f}, \quad \rho_f = \frac{\rho h}{h_f}, \qquad (31)$$

where $E_f$, $G_f$ and $\rho_f$ are respectively the fictive Young's modulus, the fictive shear modulus and the fictive density.

### 2.2.2 Problem-solving

The nonlinear problem-solving is performed with a quasi-static resolution based on a dynamic backward Euler scheme [13]. The natural damping of this scheme helps the convergence towards a nonlinear static equilibrium. An initial modal analysis computation step is also carried out to extract the first natural frequency $F_1$. This frequency is then used to set up the total computation time $T_{end} = \frac{3}{F_1}$, the initial time step $T_0 = \frac{0.05}{F_1}$ and the minimal time step value $T_{min} = T_0 \cdot 10^{-4}$, that are required for the automatic time step increments [24] of the dynamic scheme. This specific calculation time set up, which is based on the expertise of the authors, is an important point to succeed of computing high nonlinear quasi-static problems. Finally, a standard Newton-Raphson nonlinear static resolution is performed to ensure the exact static convergence.

As shown previously in Figure 3(b), the initial mainsail geometry is fully flat. This non-deformed sail is a very complex starting configuration to perform FSI computations, additionally to the fact that mainsail panels do not have any membrane pretension. To overcome this, the structural analysis requires a total of 7 steps detailed hereafter and performed prior to FSI iterations.
**Load case #1.** The mainsail head is free to move along the mast rail, all nodes are set in the same plan $(\vec{x}, \vec{z})$ and a halyard pretension is imposed;

1. Linear modal analysis;

2. Linear static analysis;

3. Nonlinear implicit dynamic quasi-static analysis with calculation time set up given by 1.;

4. Nonlinear static analysis.

**Load case #2.** The mainsail head is hooked, planar node set and the halyard pretension are removed and a low uniform pressure field is imposed on shell element faces;

5. Nonlinear modal analysis;

6. Nonlinear implicit dynamic quasi-static analysis with calculation time set up given by 5.;

7. Nonlinear static analysis.

Indeed, a direct calculation with real boundary (Load case #2.) conditions did not work because of too many numerical instabilities. To avoid this kind of issue, Load case #1. was performed first to enable a smoother start of the calculation. Finally, FSI iterations of the coupled problem (see next Section 2.3) can restart from the nonlinear converged static solution obtained in step 7.

## 2.3 Coupled problem

The FSI problem is solved by coupling together the Abaqus™ 2017 finite element software and the lifting-line flow model. The numerical procedure, achieved through implemented Python scripts, is presented below.

In the present study, an upwind point of sail is considered including a 1/6 wind power law profile [15]. The boat and true wind speeds, at a reference height of 50 m, are 10 kts and 28.9 kts respectively, the true wind angle is 60 deg and a boom angle of 27.3 deg is also taken into account. The air density is set to $\rho_0 = 1.225 \text{ kg/m}^3$.

A sail Delaunay triangulation is firstly performed on the unstructured grid mesh provided by the structural solver (see Figure 4). Horizontal $(\vec{x}, \vec{y})$ plane intersections with the triangulated sail are then used to extract a sail section distribution along $\vec{z}$. As it can be seen in Figure 4, these sections (red points) are distributed along $\vec{z}$ according to a cosine law, in order to be directly transferred to the LLT model. The $\eta_c(x)$ camber law of each section is then interpolated with Chebyshev polynomials [19], and allows the computation of $\alpha_0(\theta)$ with (10). After solving (9) and computing $\omega_i(\theta)$ with (8), the aerodynamic pressure field of each sail section was determined with the semi-analytical method proposed in Section 2.1.2. The aerodynamic pressure field is finally computed for each structural shell element centroids by using linear interpolations. An example of such $\Delta_P$ reconstruction is given in Figure 5.

The nonlinear structural solving is then performed with a restart from step 7. (see end of Section 2.2.2). This is achieved by a first nonlinear implicit dynamic quasi-static analysis and a second nonlinear static analysis. From the deformed sail shape, a new pressure field distribution can be computed, and the iterative process is repeated until a convergence is reached. In practice, FSI iterations are stopped when
<figure><img src="image_5.png" /><figcaption>Figure 4: Mainsail triangulation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Example of mainsail aerodynamic pressure field (Pa) on deformed structure.</figcaption></figure>
relative differences between consecutive iterations of different relevant macroscopic mechanical parameters decrease as a power law of the iteration number.

For this mainsail FSI problem, it is interesting to note that there is no need to use any relaxation methods on one hand, and that the typical number of FSI iterations before convergence is about 5 to 10 on the other hand.

## 3 FSI results

Numerical results of the mainsail FSI problem are now presented. The considered point of sail has been defined in Section 2.3. The purpose of these results is to illustrate the quickness and robustness of the FSI approach presented in this study. A mesh convergence study has been carried out to determine the required number of shell elements and is presented in Section 3.1. The impact of the halyard pretension has also been investigated, in terms of maximum sail displacement and reaction forces, in addition to local loop and membrane shell forces. Results are presented in the Section 3.2.

### 3.1 Structural mesh convergences

The convergence study is performed with 7 structural meshes, having 6 249, 6 755, 7 924, 9 895, 25 574, 75 418 and 126 966 shell elements respectively. For each one, the mainsail FSI problem is solved for both 5 T and 10 T (metric ton) of halyard pretension, and under a number of iterations set at 10. Results are analyzed in terms of relative differences to the finest mesh and of relative differences between 2 consecutive FSI iterations.

Figure 6 shows relative differences to the finest mesh, for both 5 T and 10 T of halyard pretension, in terms of maximum displacement "U_MAX", total strain energy "ALLSE", clew reaction force "R_CLEW", tack reaction force "R_TACK", head reaction force "R_HEAD" and side force "R_Y". A convergence trend is seen to appear at the fifth mesh (25 574 shell elements) for both halyard pretension settings.

<figure><img src="image_7.png" /><figcaption>Figure 6: Structural mesh convergences.</figcaption></figure>

Computing times needed for solving FSI problems with 10 FSI iterations, on a classical 4 CPU desktop computer, for the first to the seventh meshes are about 0.51 h, 0.56 h, 0.61 h, 0.71 h, 1.46 h, 1.94 h and
2.42 h respectively. These times are relatively low given the complex FSI problems to be solved, and they can even be reduced more by reducing the number of FSI iterations.

Indeed, as shown in Figure 7, where relative differences between 2 consecutive FSI iterations for the finest mesh are plotted, the solving of the FSI problems appears to be already converged at the fifth or sixth iteration. The robustness of the proposed approach, without a need of relaxation methods, is also

<figure><img src="image_8.png" /><figcaption>Figure 7: Convergences of FSI iterations for the finest mesh.</figcaption></figure>

highlighted in Figure 7, where straight convergence lines in log-scale are visible.

Based on convergence results and solving times, the sixth structural mesh was chosen to investigate the impact of the mainsail halyard pretension in Section 3.2.

## 3.2 Impacts of the mainsail halyard pretension

The mainsail halyard pretension is an important parameter and it can furthermore vary during sailing. The sail flying shape will depend on this pretension, as will the forces on the rig, and it must therefore be taken into account during the design process. In the present study, clew and tack mainsail points do not have any pretension.

Figure 8 shows global impacts of the halyard pretension when it varies from 1.25 T to 10 T. It should be noted that a pretension of about 1 T is very low regarding the mainsail model weight, which is presently of about 0.9 T. The maximum sail displacement is linearly decreasing of about 20% when the pretension increases (see Figure 8(a)). This is related to the general decrease of the mainsail twist, that also affects aerodynamic forces such as the side force "$R_Y$", plotted in Figure 8(b), which shows a linear increase of about 20%. The mainsail head and tack points have a different trend, as reactions forces have a linear increase from 1.25 T to approximately 5 T of pretension, before a visible change in slope. This change in slope is due to the low tension values in luff connections of bottom panels, when the pretension is below of about 5 T. Above 5 T luff connections are all fully loaded and a part of the mainsail halyard pretension can pass from the head to the tack points by the luff. The clew reaction force is linearly increasing with the halyard pretension. When the sail is under aerodynamic pressure, an internal load transfer is performed trough panel membranes from the leech to the luff. This load transfer increases
<figure><img src="image_9.png" /><figcaption>Figure 8: Global impacts of the halyard pretension.</figcaption></figure>

with the halyard pretension, that produces a faster increase in reaction forces at head and tack mainsail points in comparison to the clew point.

Figure 9 illustrates forces in panel loop connections at both luff and leech as functions of the halyard pretension and loop vertical locations. Force functions are plotted with colored surface contours thanks to linear interpolations of loop forces, that are represented by dashed black curves. This representation

<figure><img src="image_10.png" /><figcaption>Figure 9: Loop forces for different halyard pretensions.</figcaption></figure>

allows to visualize the increase of low luff force values in bottom panels locations, when the halyard pretension is increasing (see Figure 9(a)). For low halyard pretensions, loop forces in the last third of the mainsail height appear to be almost constant for both luff and leech. Then, as long as the halyard pretension increases, loop forces in this mainsail region are no longer constant, and the luff presents a local maximum of loop forces at a height of about 25 m (≈ 63% of the overall height), while it is a local minimum of loop forces for the leech. This result shows the importance of considering an advanced
nonlinear mainsail model to correctly capture nonlinearities that occur in structural loads. The loop design must take into account both halyard pretension and loop mainsail location.

Internal loads of mainsail panel membranes are presented in the next following Figures, where two halyard pretensions of 1.25 T and 10 T are considered. Figure 10 shows $N_{11}$ membrane stress forces, direction 1 is the projection of the global X axis onto the deformed membrane surface. The first half height panels appear to be more $N_{11}$ loaded and the load transfer increases from the leech to the luff are clearly visible when the halyard pretension change from 1.25 T to 10 T. A diagonal direction loading, starting from the leech, is also shown for the first bottom panels. It is a typical result of panel wrinkles phenomenon.

<figure><img src="image_11.png" /><figcaption>Figure 10: Membrane stress forces $N_{11}$ (N/m) on deformed structure.</figcaption></figure>

Figure 11 shows $N_{22}$ membrane stress forces, where direction 2 is deduced from 3 × 1, where direction 3 denotes the normal of the deformed membrane surface. The loop forces distributions previously shown in Figure 9 are also qualitatively described. In addition to $N_{22}$ loads arising from loop connections, it is shown that a part of $N_{22}$ goes through synchronization straps, in the first half mainsail height.

## 4 Conclusion

A fast and robust strategy to solve nonlinear FSI on yacht sails has been presented in this paper. This study is part of the "Solid Sail 2.0" industrial project, who's purpose is to develop new designs of rigid sails. A 300 m² semi-rigid composite mainsail prototype, composed of 21 articulated panels, was numerically modeled and the proposed solving strategy was applied.

For the flow solver, an implemented lifting-line model including a semi-analytical method for the aero-dynamic pressure field computation is proposed. This model relies on an efficient Chebyshev-Gauss quadrature to approximate the two-dimensional vortex distribution of sail sections. Pressure distributions at each section are thus quickly determined and the overall mainsail pressure field is computed using
<figure><img src="image_12.png" /><figcaption>Figure 11: Membrane stress forces $N_{22}$ (N/m) on deformed structure.</figcaption></figure>

linear interpolations. Mainsail panels are modeled by thin shell elements including a specific thickness approach to enhance the numerical stability that is possibly affected by membrane winkles. The non-linear solving is performed by a quasi-static resolution with a dynamic backward Euler scheme. The assembly of these flow and structural methods leads to a fast and robust approach to compute nonlinear FSI on yacht sails, as no relaxation methods are needed.

A mesh convergence study has been carried out to determine the required number of shell elements for computing the mainsail prototype FSI problem. Results have shown the computing efficiency of the proposed method in terms of both time and robustness. Impacts of the mainsail halyard pretension have also been investigated. The pretension varied from 1.25 T to 10 T, and overall dependent results have been shown, in addition to the nonlinear response of structural connection loops. These results demonstrated the robustness of the nonlinear structural method, in particular for computations with low halyard pretension values.

An experimental campaign definition with a 50 m² composite mainsail is on progress. Methods described in this paper will be applied to compute the structural response of this sail and comparisons to experiments of numerical results will be carried out.

Future developments will concern the aerodynamic jib effects and mast deformations [4] considerations. The two-dimensional RANS/lifting-line method proposed by Graf et al. [12] seems also promising and could be well adapted to the present work too. The overall mainsail design may also be improved, for example with several weight reductions, through advanced optimization methods [26].

## Acknowledgments

This study takes part in the "Solid Sail 2.0" project funded by the "Région Bretagne", and approved by the "Pôle EMC2" and the "Pôle Mer Bretagne Atlantique". This project involves a consortium of 3
industrial partners, "Chantiers de l'Atlantique", "GSea Design" and "Multiplast" in collaboration with "ENSTA-Bretagne" a graduate and post-graduate school of engineering as an academic partner. The "Chantiers de l'Atlantique", "GSea Design" and "Multiplast" companies are gratefully acknowledged by the authors for their technical support and successful collaboration.

## References

[1] I. H. Abbott and A. E. Von Doenhoff. *Theory of wing sections, including a summary of airfoil data*. Courier Corporation, 1959.

[2] N. Aubin, B. Augier, P. Bot, F. Hauville, and R. Floch. Inviscid approach for upwind sails aerodynamics. how far can we go? *Journal of Wind Engineering and Industrial Aerodynamics*, 155:208–215, 2016. ISSN 0167-6105.

[3] B. Augier, P. Bot, F. Hauville, and M. Durand. Experimental validation of unsteady models for fluid structure interaction: Application to yacht sails and rigs. *Journal of Wind Engineering and Industrial Aerodynamics*, 101:53–66, 2012.

[4] S. Bak and J. Yoo. Fsi analysis on the sail performance of a yacht with rig deformation. *International Journal of Naval Architecture and Ocean Engineering*, 11(2):648–661, 2019. ISSN 2092-6782.

[5] S. Bak, J. Yoo, and C. Y. Song. Fluid-structure interaction analysis of deformation of sail of 30-foot yacht. *International Journal of Naval Architecture and Ocean Engineering*, 5(2):263–276, 2013.

[6] G. Bezine. On a method of comparison for plate elements in finite element engineering software programs. *Mechanics Research Communications*, 29(1):35–43, 2002. ISSN 0093-6413.

[7] M. Durand. *Interaction fluide-structure souple et légère, application aux voiliers*. PhD thesis, Ecole Centrale de Nantes (ECN), 2012.

[8] M. Durand, A. Leroyer, C. Lothodé, F. Hauville, M. Visonneau, R. Floch, and L. Guillaume. FSI investigation on stability of downwind sails with an automatic dynamic trimming. *Ocean Engineering*, 90:129-139, 2014.

[9] Z. Eshkuvatov, N. N. Long, and M. Abdulkawi. Approximate solution of singular integral equations of the first kind with cauchy kernel. *Applied Mathematics Letters*, 22(5):651–657, 2009. ISSN 0893-9659.

[10] M. Fluck, F. C. Gerhardt, J. Pilate, and R. G. J. Flay. Comparison of potential flow-based and measured pressure distributions over upwind sails. *Journal of Aircraft*, 47(6):2174-2177, 2010.

[11] H. Glauert. *The elements of aerofoil and airscrew theory*. Cambridge University Press, 1983.

[12] K. Graf, A. Hoeve, and S. Watin. Comparison of full 3d-rans simulations with 2d-rans/lifting line method calculations for the flow analysis of rigid wings for high performance multihulls. *Ocean Engineering*, 90:49–61, 2014. ISSN 0029-8018. Innovation in High Performance Sailing Yachts - INNOVSAIL.
[13] E. Hairer, S. P. Nørsett, and G. Wanner. *Solving ordinary differential equations I: nonstiff problems*, volume 8. Springer Science & Business Media, 2008.

[14] P. Heppel. Accuracy in sail simulation: Wrinkling and growing fast sails. In *High Performance Yacht Design Conference*, pages 4–6, 2002.

[15] S. Hsu, E. A. Meindl, and D. B. Gilhousen. Determining the power-law wind-profile exponent under near-neutral stability conditions at sea. *Journal of Applied Meteorology*, 33(6):757–765, 1994.

[16] J. Katz and A. Plotkin. *Low-Speed Aerodynamics*. Cambridge Aerospace Series. Cambridge University Press, 2001.

[17] F. W. Lanchester. *Aerodynamics*, 1907.

[18] A. Maison, A. Nême, and J. B. Leroux. De la problématique de dimensionnement de grands kites. In *115ème session de l'Association Technique Maritime et Aéronautique (ATMA)*, Paris, 2017.

[19] J. Mason and D. Handscomb. *Chebyshev Polynomials*. CRC Press, 2002.

[20] A. D. Polyanin and A. V. Manzhirov. *Handbook of integral equations*. CRC press, 1998.

[21] L. Prandtl. Theory of lifting surfaces. Technical report, NASA, 1920.

[22] R. Ranzenbach, D. Armitage, and A. Carrau. Mainsail planform optimization forIRC 52 using fluid structure interaction. In *The 21st Chesapeake Sailing Yacht Symposium, SNAME*, 2013.

[23] J. N. Reddy. *An Introduction to Nonlinear Finite Element Analysis: with applications to heat transfer, fluid mechanics, and solid mechanics*. OUP Oxford, 2014.

[24] D. F. Rossi, W. G. Ferreira, W. J. Mansur, and A. F. G. Calenzani. A review of automatic time-stepping strategies on numerical time integration for structural dynamics analysis. *Engineering Structures*, 80:118–136, 2014. ISSN 0141-0296.

[25] M. Sacher, F. Hauville, R. Duvigneau, O. Le Maître, N. Aubin, and M. Durand. Efficient optimization procedure in non-linear fluid-structure interaction problem: Application to mainsail trimming in upwind conditions. *Journal of Fluids and Structures*, 69:209–231, 2017.

[26] M. Sacher, M. Durand, Élisa Berrini, F. Hauville, R. Duvigneau, O. L. Maître, and J.-A. Astolfi. Flexible hydrofoil optimization for the 35th america's cup with constrained ego method. *Ocean Engineering*, 157:62–72, 2018. ISSN 0029-8018.

[27] S. P. Timoshenko and S. Woinowsky-Krieger. *Theory of plates and shells*. McGraw-hill, 1959.

[28] D. Trimarchi, S. Turnock, D. Chapelle, and D. Taunton. Fluid-structure interactions of anisotropic thin composite materials for application to sail aerodynamics of a yacht in waves. In *The 12th Numerical Towing Tank Symposium, Italy*, pages 4–6, 10 2009.

[29] D. Trimarchi, M. Vidrascu, D. Taunton, S. Turnock, and D. Chapelle. Wrinkle development analysis in thin sail-like structures using mitc shell finite elements. *Finite Elements in Analysis and Design*, 64:48–64, 2013. ISSN 0168-874X.
[30] I. M. Viola and R. G. Flay. Sail aerodynamics: Understanding pressure distributions on upwind sails. *Experimental Thermal and Fluid Science*, 35(8):1497–1504, 2011. ISSN 0894-1777.

[31] W. Wong and S. Pellegrino. Wrinkled membranes iii: numerical simulations. *Journal of Mechanics of Materials and Structures*, 1(1):63–95, 2006.