We consider the first order diffractionradiation problem of a ship moving on the horizontal plane. The two dimensional movement of the ship is identified by , where and are the linear and angular velocity of the ship. We assume that the flow is incompressible and irrotational, so that the velocity field , is derived from a potential function as: . The following assumptions are made on the order of magnitude of the velocity components and free surface elevation :

(1) 
Based on the previous assumption, the governing equations for the first order diffractionradiation wave problem are:

(2) 
where is the gradient in the horizontal plane, is the fluid domain, represents the wetted surface of the ship, is the water depth, is the local ship velocity of a point over the wetted surface, is the normal of the ship wetted surface pointing outwards the ship, is the water current, and is the free surface pressure. Once the velocity potential is obtained, the pressure field is obtained straightforward using the Bernoulli equation:

(3) 
In order to solve the governing equations provided in Eq.(2), a velocity potential decomposition is introduced. Let and be such that;

(4) 
Then, we can split the governing equations into the following sets of equations:
Set 1:

(5) 
Set 2:

(6) 
The first set of equation, Eq.(5), has an analytical solution (known as Airy waves):

(7) 
where are the wave amplitudes; are the wave frequencies; are the wave numbers; are the wave phases. From this point on, we will refer to as the incident potential and as the scattered potential. Also the dispersion relation for Airy waves holds:

(8) 
From Airy waves theory, we know that and . Then, the terms , , and . Neglecting these terms in Eq.(6), and using the analytical solutions of Airy waves, the governing equations for the scattered waves potential becomes:

(9) 
and

(10) 
Notice that the terms and accounts for the deviation of the incident Airy waves due to the fact that the incidents waves are transport by a nonuniform flow field. Also, the terms , , and are not subject to any kind of linearization, so that there will be no restriction in using nonsteady base flows.
From a numerical point of view, based on the numerical approaches used in this work, it is convenient to solve Eqs. (9)(10) in a fix frame of reference rather that on the global frame of reference. Therefore, the aforementioned equations will be solved in a frame of reference fix to the ship. Figure 1 shows the global and local frame of reference. This frame of reference is assumed to match the global frame at time zero. For a observer sitting in the ship, the flow field around the ship will be . Therefore, the governing equations in the local frame of reference become:

 style="width: 5px;textalign: right;whitespace: nowrap;"  (11) }
and

(12) 
where the incident wave potential and incident wave elevation must be transformed to the local frame of reference.
The flow is assumed to be incompressible and irrotational. Therefore the Laplace equation has to be fulfilled within the fluid domain. We adopt the finite element method for solving the Laplace equation in the fluid domain.
This formulation has been developed to be used in conjunction with unstructured meshes. The use of unstructured meshes enhances geometry flexibility and speed ups the initial modelling time. An automatic unstructured grid generator based on the advancing front method is used to generate triangular surface grids and tetrahedral volumetric meshes.
Let be the finite element space to interpolate functions, constructed in the usual manner. From this space, we can construct the subspace , that incorporates the Dirichlet conditions for the potential . The space of test functions, denoted by , is constructed as , but with functions vanishing on the Dirichlet boundary. The weak form of the problem can be formulated as:
Find , by solving the discrete variational problem:

(13) 
where , , and are the potential normal gradients corresponding to the Neumann boundary conditions on body, radiation boundary, free surface and bottom, respectively. At this point, it is useful to introduce the associated matrix form

(14) 
where is the standard laplacian matrix, and , and are the vector resulting of integrating the corresponding boundary condition term. The bottom boundary is imposed in a natural way when using FEM just doing . From this point forward we are omitting .
Solving the free surface boundary condition efficiently is the key point when dealing with water waves problem. The free surface conditions can be rewritten as:

(15) 
where is the base flow. Linearization respect to this based flow is quite common. Linearizations most commonly used are: the Kelvin linearization, which assumes that and hence ; and the double body, which assumes that , where and , being the scattered velocity potential obtained when solving the double body problem. As we mentioned before, in this work no linearization is assumed since it is considered that . Therefore, the convective velocity will be different in each time step, and the numerical scheme will be adapted accordingly. Three main issues must be kept in mind:
The numerical schemes adopted for solving the kinematicdynamic free surface boundary conditions are based on AdamsBashforthMoulton schemes, using an explicit scheme for the kinematic condition, and implicit one for the dynamic condition. Then can be imposed as a Dirichlet Boundary condition. The schemes read as follows:

(16) 
It was pointed out before that care must be taken when solving the free surface conditions, especially when evaluating the convective terms. In this work, the convective term is obtained by differentiating along streamlines:

(17) 
where denotes the derivative along the streamline. This streamline derivatives is estimated using a two points upstream and one point downstream differential operator inspired by the quickest scheme [1]. Figure 2 shows the tracing of the streamline at node C. The left (L) and forward left (FL) points are the upstream points, while the right (R) point corresponds to the downstream point. The values of the scattered velocity potential and scattered free surface elevation at L, FL and R points are obtained by linear interpolation between the nodes of the edges where they lie on. The stream line differential operator reads as:

(18) 
where the stencils are:

(19) 
While most potential codes based on boundary element methods rely on structured meshes, the above scheme allows SeaFEM to use unstructured meshes, which provides with high flexible to simulate complex geometries, as well as helps to reduce the number of elements by refining only in those areas of interest. Figure 3 provides a sample of an unstructured mesh used for SES simulations.
When modelling Aircushion vehicles (ACV), it is important to model properly the effect of the aircushion over the free surface and the vehicle. In this work, the aircushion is modelled as a pressure distribution over a specified free surface. Then evaluations of this pressure field on specific locations are introduced in Eq.(16) straightforward. Further details about the pressure field generated by the aircushion will be given later on.
Wave absorption is necessary in order to avoid wave reflection at the edges of the computational back to the ship. A simple an effective wave absorber will be used, based of transferring wave energy to the atmosphere. This is carry out by modifying the pressure over the free surface adding a term such that , where is an absorption factor. By setting , no wave absorption will be imposed. Then, increasing the pressure when the vertical component of the free surface velocity is positive, or vice versa, the free surface has to make a virtual work over the atmosphere, transferring energy from waves to the atmosphere.
The seals of air cushion vehicles are a very important part since keeping the pressure inside the aircushion depends on their interaction with the free surface. Because of this interaction with the free surface, modelling the behaviour of the seals is not a trivial issue.
In this work, the free surfaceseal interaction problem is formulated as finding a pressure field to be applied over the free surface such that the elevation of this one is limited by the location of the seal. That is to say, the seal act as an upper limit for the free surface elevation.
The free surface boundary conditions are applied in different ways depending on whether the free surface is in contact with the seal or not. If the free surface is not in contact, the boundary conditions are applied as if there is no seal, but if there is contact, the implementation will be different in order to ensure that the free surface does not penetrate de seal and the necessary pressure to fulfil this condition is calculated.
It will be said that the free surface node where the algorithm is to be applied is dry if the seal is not in contact with the free surface at that location, and wet if it is. Figure 5 illustrates the wet and dry concepts.
The main challenge for an algorithm like this is to be capable of capturing when a node goes from dry to wet and vice versa, as well as estimating the pressure field on the wet nodes. Figure 7 shows the wetting and drying algorithm. For a dry node, the implementation of both, the kinetic and dynamic boundary condition is the same as for any other node not interacting with the seal. However, for wet nodes, the free surface boundary condition is imposed via imposing that the free surface elevation matches the seal elevation, and ensuring that there is no flow across the seal. These two conditions are represented by the following equations:

(20) 

(21) 
The change from being a dry node to become a wet node is identified via the kinematic BC though the condition . On the other hand, the switch from being a wet node to become a dry one is carried out by comparing the dynamic pressure with the reference pressure plus a detachment condition. This detachment condition requires for the free surface to detach at a specific node that this node is not completely surrounded by attached nodes. Should this be the case, there would be not path for the air to move in. Therefore, if there is no connection with the air, pressure might drop below the reference pressure, which might even predict the happening of cavitations. For ACVs, the reference pressure will be the aircushion pressure. The dynamic pressure on wet nodes is obtained by the dynamic BC:

(22) 
In order to pass information from the hydrodynamics solver to the structural solver and viceversa, a strategy has been devised. First, an interpolation library is used to perform the coupling. This interpolation library is capable of interpolating data (scalar or vectorial fields) from one mesh to another mesh. The meshes need not be matching, but should be discretizations of the same surface. The interpolation library was developed by Compass and is a general purpose C++ code provided for this project.
Secondly, the interaction is performed using a TCL script invoked from the hydrodynamics solver. This script uses the interpolation library to pass pressure information from the hydrodynamics mesh to the structural mesh, and vertical displacements from the structural mesh to the hydrodynamics mesh.
This strategy assumes that the structure moves only in the vertical direction. This is not what the structure does, but is a good enough approximation to predict the forces of the water on the seals and the free surface elevation in the seals region.
In the Exhibit 1 the TCL script is shown along with the documenting comments. As an indication, we shall explain that 2 auxiliary meshes are used in this interpolation process. This is because the meshes of the structure and the coupling free surface are not coincident in geometry. That is, the mesh of the free surface lies at the flat free surface of the water. This is natural given the solving strategy for the fluid. Likewise, the mesh of the structure lies outside of the flat free surface of the water, as the seals are designed to seal the air cushion above the water.
Therefore, in order to perform the interpolation, we use an auxiliary mesh that matches the geometry of the seals to interpolate with the structure. And then a duplicate of that mesh projected to the plane Z=0. This duplicate mesh serves to interpolate with the hydrodynamics solver. The interpolation between the two auxiliary meshes is direct, as they are duplicates of each other. They keep the same nodes and numbering.
The boundary condition to be applied on the body subject to study is the usual no flux across the body surface. For this purpose, a normal flux is induced via a Newman boundary condition in order to cancel out the normal component of the incident velocity plus the velocity of the body at the specific location.
However, when considering transom sterns, flow detachment happens at the lower edge of the stern. While potential flow is incapable of predicting this sort of detachment, a transpiration model will be used to enforce it. To do so, a no flux boundary condition is not used. On the contrary, a flux is allowed in order to enforce that the detachment edge belong to the free surface stream surface. Figure 6 explains the transpiration model idea.
When designing an ACV, the aircushion plays an important role since it will carry some 90% of the weight, reducing significantly the total drag by reducing the viscous drag. On the other hand, the pressure distribution along with the free surface elevation inside the aircushion, will induce forces and moments in all directions. In fact, the resulting force in the forward direction is an extra drag that corresponds to the waves generated by the pressurized aircushion. These forces and moments induced over the ACV can be obtained by integrating the pressure over the free surface covered by the aircushion.

(23) 
These forces and moments have to be used to solve the dynamic of the ACV along with other loads. However, obtaining the pressure distribution within the aircushion is not trivial, and would require solving the dynamics of the aircushion, accounting for the blowers, airleaks underneath the seals and hulls, and solving the airflow inside the aircushion. Since solving the aircushion dynamics is out of the scope of this work, the aircushion will be modelled as a pressure distribution that, while simple, can be tuned into different pressure distribution aiming to be capable of reproducing a variety of models. The pressure distribution used in this work is as follows:

(24) 
where is the pressure at location and time t, is an average pressure value, represents temporal variations of the pressure, and represents a spatial distribution. In our implementation, is a constant and uniform value; is time dependent and it might depend on variables such that the volume inside the aircushion (for instant, to implement compressibility effects), airflow displaced by the free surface, etc; and is spatially variable and is meant to depend on the coordinates (x,y).
The software has been equipped with a new type of boundary condition in order to be capable of accounting for a pressurized free surface just like an aircushion. Figure 8 shows how the aforementioned aircushion model is easily implemented in the software interface.
The ship dynamics consist of using Newton’s laws to obtain the ship movements once external loads acting on it are known. In this work, it is assume that the ship trajectory is fixed in the horizontal plane, constraining the surge, sway and yaw movements. On the other hand, heave, roll, and pitch will be free.

(25) 
where is the mass matrix, is the restoring matrix, and represents all kind of external forces and moments. Eq. (25) can be written as:

(26) 
The equation of motions, Eq.(26), are solved using the well known Newmark´s scheme with parameters and (also known as the average acceleration method).

(27) 
When solving the body dynamics, care must be taken since not all loads are considered by default. Therefore, all those that are not considered by default, must be introduced externally via the available function editor for external forces and moments. By default, it is accounted for the weight and the buoyancy of the body, the hydrostatic restoring forces and moments, and the dynamic pressure integration over the wet surfaces. Any other loads, like those induced by the aircushion must be introduced manually via de function editor. Figure 9 shows an example of how external pitching moments are introduced within the interface. The aircushion pitching moment, for instance, is accessible via the My_Pfs variable. A wide range of internal variables are accessible for the function editor, which allows the implementation of a wide range of analytical models.
For the analysis of the flexible seals, two different alternatives were considered. First, to finish the development of a new rotationfree shell triangle [2] that holds the promise of avoiding numerical locking for very thin shells while maintaining accuracy. Second, to use an already developed membrane triangle [3].
The development of the shell element in the nonlinear regime progressed more slowly than expected, and did not reach a sufficient mature status before the end of the project. Therefore, the focus turned towards the membrane solver.
The preference was set rather in favour of the shell element because for a coupled problem of this nature, an implicit solver is preferred. And for this, it is known that membrane solvers usually suffer from singularities in the system matrix.
Nevertheless, the membrane element has finally been used (in an implicit scheme), as it is expected to perform well at least in the modelling of the stern seals that are mostly subjected to tension state stress. The membrane solver uses a total Lagrangian formulation and assumes a linearelastic behaviour of the material. The kinematic description is formulated in large displacements and large deformations, making use of the right CauchyGreen deformation tensor. For further reference, see the reference [3].
The loads on the membrane are of two types. First, let’s consider the loads imposed by the air cushion system. The air cushion system acts differently on the bow skirt and the stern seal. On the bow skirt, the air cushion applies a direct follower load on the inner side of the skirt. On the stern seal, a differential pressure between the inside of the lobes and the air cushion maintains the seals taut. All these loads are considered as follower loads. A follower load has the property of not maintaining the symmetry of the system matrix for the implicit solver. Therefore, we need to choose an adequate numerical solver for the system of equations.
Second, the loads applied by the water by means of the interaction algorithm. The structural solver has been modified in order to consider the loads applied by the liquid as follower loads. Failing to include this consideration makes the system as a whole to diverge.
The contacts between the bow skirt stiffeners have been modelled as elastic contacts. This allows for a more realistic behaviour than considering the stiffeners stitched together and at the same time keeps the computational cost controlled.
First of all, let’s define the local frame of reference to be used. Figure 10 explains the location of the frame of reference used. The origin of coordinates is located at the intersection of the middle plane, the plane of the transom stern, and the water plane.
The XR1B surface effect ship is used for this demonstration analysis. Figure 11 shows the geometry of the SES, and Table 1 provides information regarding the particular of the ship, as well as of the design conditions.
Length overall (m)  76.35 
Length waterline on cushion (m)  67.52 
Draft on cushion (m)  1.33 
Beam max. (m)  22.25 
Cushion width (m)  16.50 
Cushion length (m)  67.14 
Displacement (tonnes)  1560 
LCG from wet deck transom (m)  33.75 
TCG from starboard of centreline (m)  0 
VCG above water line (m)  2.3564 
Pitch gyradius (m)  21.6 
Roll gyradius (m)  8.01 
Aircushion average pressure (KPa)  11 
Water depth (m)  10 
In order to carry out simulations, an aircushion model is needed. In this work, we will adopt an analytical aircushion model of the type [4]:

(28) 
The above pressure distribution accounts for a smooth pressure drop in the areas of the seals as well as when getting closer to the hull. Figure 12 shows the effect of the aircushion pressure on the free surface at zero speed.
Let’s mention that the use of this specific aircushion is arbitrary, and results will highly depend on how well modelled the aircushion is. However, for the purpose of this demonstration, this arbitrary model is good enough since the purpose is to show how to use the tools developed, rather than modelling a specific real case. This model could be easily improved by using any experimental data of pressure distribution as well as pressure variations along with its governing equations. Should the aircushion behaviour be better known, then a better model could be implemented in the code. However, no realistic and verified model has been developed yet.
A number of forces are acting simultaneously when a surface effect ship is navigating. Figure 13 shows these forces. A brief description of each component is given next:
More complex model can be found in the literature to model some of the previously described loads. Also, there are other external loads acting over a surface effect ship, such that aerodynamic momentum drag, airleakage drag, and others [5]. However, for the sake of simplicity and for the demonstration purpose of this work, there is no need of including them in the analyses. In case of including them, they should be estimated in a reliable manner, which is out of the scope of this work. On the other hand, we will pay special attention to the estimation of the aircushion wave making loads and the loads resulting from the freesurface seal interaction, which are the main objectives of this work.
Moving pressurized free surface results in a wave train. Therefore, since some energy has been transferred to the fluid in the form of waves, a drag must be induced in the opposite direction to the movement of the pressurized free surface. This drag can be obtained via integration of the pressure distribution (excess of pressure over the atmospheric pressure) over the free surface. In case that the pressurize free surface is created by an ACV, the forces and moments acting over the ship are obtained via Eq.(23). Being capable of measuring these loads is a key point in the design of ACVs, and one of the main objectives of this work.
In [6] analytical solutions are found for the drag induced by several pressure distribution moving with constant forward speed in a canal. These solutions are linear and based on the Kelvin linearization. Figure 14 compares the wave drag coefficient induced by a rectangular pressure distribution (LxB) obtained by the proposed methodology versus the analytical solution provided in [6]. It can be observed that numerical results agree quite well overall.
The lower the Froude number (or the higher the dimensionless number ), the higher are the difference between the numerical and analytical solutions. This discrepancy is not surprising since the waves generated by a pressure patch becomes more nonlinear at lower Froude numbers [4], leading to more complex flows to be solved numerically. In fact, the lower the Froude number, the more complicated is obtaining a steady solution. In particular, at Froude number equal to 0.2, an unsteady wave is generated at the stern edge of the pressure patch. This wave moves back and forward in time, following an unsteady, but cyclical, process where at some times the wave becomes too steep to be solved under the numerical approach of this work. This steep wave generates some numerical dispersion that dissipates afterwards while the wave starts smoothing out till the cycle starts over. Despite of the limitations of the numerical approach used in here, it is capable of predicting the happening of these physical unstabilities, which might be critical in the design of ACVs.
In this section, a first analysis on the XR1B navigation properties will be carried out. In a first step, numerical simulations of the XR1B with forward speed and all its degrees of freedom constraint are carried out. The main objective of these simulations is to obtain the loads induced by the dynamic pressure of the water over the hull as well as by the aircushion. All the other loads are either known or estimated via analytical models. Table 2 shows the drag, vertical forces and pitching moments
Dynamic pressure loads  Aircushion loads  
0  0  0  0  0  10388780  20660000 
0.1  50200  220000  600000  70000  10388780  20550000 
0.2  60000  250000  2800000  280000  10388780  19900000 
0.3  66000  1010000  1270000  220000  10388780  20135800 
0.4  112400  1320000  33550000  503600  10388780  19407100 
0.6  95810  533400  26000000  282390  10388780  19834520 
0.8  117594  1059660  12136880  276460  10388780  19816880 
1.2  178965  1447142  3634880  183518  10388780  20142400 
1.6  260156  1545982  14873440  116532  10388780  20354600 
2  358270  1480790  22941400  76991  10388780  20468000 
The viscous drag will be estimated following ITTC directions using the ITTC friction line.

(29) 
where is the dimensionless viscous drag coefficient, is the fluid density, is the wet surface, and is the form factor which, for this demonstration, will be assumed to be equal to 0.025. will be assumed to be horizontal and located at a height of . Then, the viscous load will induce a pitching moment .
The seal loads are mainly frictional loads. In this first approximation, we will estimate using the ITTC viscous drag coefficient. The wet surface is estimated as , where is the aircushion width, and is the wet length of the seal. The wet length for the bow and stern seals will be assumed to be and respectively. Both frictional forces will be assumed to be horizontal and acting on the z=0 plane. Therefore, they will induce a pitching moment .
0  0  0  0  0  0 
0.1  1.75 10^{8}  1.972 10^{3}  3900  0  11783 
0.2  3.5 10^{8}  1.795 10^{3}  14197  0  42894 
0.3  5.25 10^{8}  1.702 10^{3}  30291  0  91521 
0.4  7 10^{8}  1.641 10^{3}  51903  0  156819 
0.6  1.05 10^{9}  1.559 10^{3}  110997  0  335365 
0.8  1.4 10^{9}  1.505 10^{3}  190488  0  575540 
1.2  2.1 10^{9}  1.434 10^{3}  408231  0  1233428 
1.6  2.8 10^{9}  1.386 10^{3}  701596  0  2119804 
2  3.5 10^{9}  1.351 10^{3}  1068261  0  3227644 
0  0  0  0  0  0 
0.1  2.98 10^{6}  3.75 10^{3}  2.36E+02  0.00E+00  5.55E+02 
0.2  5.95 10^{6}  3.29 10^{3}  8.28E+02  0.00E+00  1.95E+03 
0.3  8.93 10^{6}  3.06 10^{3}  1.73E+03  0.00E+00  4.08E+03 
0.4  1.19 10^{7}  2.91 10^{3}  2.93E+03  0.00E+00  6.90E+03 
0.6  1.79 10^{7}  2.72 10^{3}  6.16E+03  0.00E+00  1.45E+04 
0.8  2.38 10^{7}  2.59 10^{3}  1.04E+04  0.00E+00  2.46E+04 
1.2  3.57 10^{7}  2.43 10^{3}  2.20E+04  0.00E+00  5.19E+04 
1.6  4.76 10^{7}  2.33 10^{3}  3.75E+04  0.00E+00  8.83E+04 
2  5.95 10^{7}  2.25 10^{3}  5.66E+04  0.00E+00  1.33E+05 
0  0  0  0  0  0 
0.1  6.75 10^{6}  3.22 10^{3}  4.59E+02  0.00E+00  1.08E+03 
0.2  1.35 10^{7}  2.85 10^{3}  1.63E+03  0.00E+00  3.83E+03 
0.3  2.03 10^{7}  2.66 10^{3}  3.42E+03  0.00E+00  8.06E+03 
0.4  2.7 10^{7}  2.54 10^{3}  5.81E+03  0.00E+00  1.37E+04 
0.6  4.05 10^{7}  2.39 10^{3}  1.23E+04  0.00E+00  2.89E+04 
0.8  5.4 10^{7}  2.28 10^{3}  2.08E+04  0.00E+00  4.91E+04 
1.2  8.1 10^{7}  2.15 10^{3}  4.41E+04  0.00E+00  1.04E+05 
1.6  1.08 10^{8}  2.06 10^{3}  7.53E+04  0.00E+00  1.77E+05 
2  1.35 10^{8}  2 10^{3}  1.14E+05  0.00E+00  2.68E+05 
The aerodynamic profile drag will be estimated as . We will assume a projected area , and a drag coefficient of . will be assumed to be horizontal and located at . Therefore, the induced pitching moment is .
0  0  0  0 
0.1  313  0  1139 
0.2  1250  0  4555 
0.3  2813  0  10248 
0.4  5000  0  18218 
0.6  11250  0  40991 
0.8  20000  0  72872 
1.2  45000  0  163962 
1.6  80000  0  291488 
2  125000  0  455450 
The propulsion system has to provide the thrust required to compensate the drag. Therefore, the horizontal thrust to be provided can be estimated by summing up all the drags given above. We will assume a propulsion system to be equivalent to an horizontal force (thrust) acting on the intersection of the base line and the stern, that is to say, acting at (XE, YE, ZE)=(0,0,1.33). Then, the propulsion system will induce a pitching moment of
0  0  0  0 
0.1  125107  0  461194 
0.2  357901  0  1319366 
0.3  324256  0  1195338 
0.4  681637  0  2512787 
0.6  518857  0  1912715 
0.8  635830  0  2343923 
1.2  881890  0  3251001 
1.6  1271010  0  4685453 
2  1799020  0  6631907 
Figure 15 and Figure 16 show the estimated drag components and cumulative drag respectively. It can be observed that the aircushion drag is especially important at lower Froude numbers. However, as the Froude number increases, frictional drag of the hulls becomes dominant, followed by wave making and seal drags, while the aircushion drag becomes less important.
In order to obtain the equilibrium position (sinkage and trimming angle), of the XR1B for each velocity condition, Eq. (25) must be solved in the absence of the inertial term:

(30) 
In order to obtain the equilibrium sinkage and trimming when the XR1B is moving forward in the absence of waves, it is necessary to know the vertical forces and pitching moments acting on the ship, as well as the hydrostatic properties of the ship. Figure 17 and Figure 18 shows the vertical forces and pitching moments at different Froude numbers, and Table 8 provides the hydrostatic properties.
Notice that based on the local reference system of the XR1B, negative pitch moments induce sinkage of the bow and negative values of the pitch angle. However, we will say that the trim angle is negative when it induces the sinkage of the bow. That is to say, the trim angle is the pitch angle but in opposite direction.
Solving Eq. (30) for the loads estimated above, we obtain the sinkage and trimming for each Froude number. Figure 19 and Figure 20 Shows the evolution of sinkage and trimming across a range of Froude numbers.
A set of numerical tests for a range of wave periods and wave directions have been carried out for Fr=1.6 and forward speed. The aim of these tests is to estimate response amplitude operators (RAOs) in heave, roll and pitch, so that complicated navigating conditions could be anticipated. The wave direction is measured respect to the OX axis and positive clockwise. That is to say, a direction of 180º means head seas, while 0º means following seas.
Figure 30, Figure 31, and Figure 32 provide the results from the numerical test. It can be observed that stern quartering waves with wave periods around 9s10s are the most dangerous ones out of the range analyzed. Strong resonance effect might appear leading to an unstable behaviour of the ship.
In order to assess the manoeuvring capabilities of the XR1B, a set of test runs performing a turning circle have been carried out at Fr=1.6. The test runs have been performed in different conditions, based on the depth, turning diameter and drift angle (see Figure 33).
The aim of this test is threefold: First, to evaluate the seakeeping behaviour while manoeuvring at different turn diameters; second, to evaluate the influence of the water depth; and third, to evaluate the influence of the drift angle of the ship.Figure 34, Figure 35, and Figure 36 show the wave pattern for a turn with D=700m for each test.
Table 9 contains the seakeeping results obtained. While no big differences are observed between conditions 1 and 2, big differences are observed from these ones with respect to condition 3. The effect of the drift basically consists of introducing a side slip velocity. This side slip velocity, which is in the order of 9% of the forward speed for a 5º drift, highly modifies the flow field around the ship, inducing a much different dynamic pressure distribution which induces larger movements.
D
(m) 
Condition 1  Condition 2  Condition 3  
H=10m  Drift=0º  H=100m  Drift=0º  H=10m  Drift=5º  
Heave
(m) 
Roll
(º) 
Pitch
(º) 
Heave
(m) 
Roll
(º) 
Pitch
(º) 
Heave
(m) 
Roll
(º) 
Pitch
(º)  
300  0.948  2.509  0.320  1.031  2.419  0.406  5.059  4.815  12.694  
350  0.671  1.995  0.323  0.739  1.912  0.414  2.483  2.426  0.683  
400  0.495  1.626  0.311  0.556  1.552  0.409  1.183  1.469  7.848  
500  0.289  1.178  0.254  0.346  1.123  0.360  0.109  0.854  14.227  
700  0.123  0.800  0.206  0.177  0.762  0.309  0.416  0.767  17.768  
1000  0.033  0.519  0.158  0.087  0.499  0.281  0.491  0.943  18.593  
inf  0.033  0  0.118  0.020  0  0.299  0.166  1.471  15.642 
D
(m) 
Condition 1  Condition 2  Condition 3  
H=10m  Drift=0º  H=100m  Drift=0º  H=10m  Drift=5º  
Fx
(T) 
Fy
(T) 
Mz
(Tm) 
Fx
(T) 
Fy
(T) 
Mz
(Tm) 
Fx
(T) 
Fy
(T) 
Mz
(Tm)  
300  21.51  10.20  7158.41  21.16  13.36  7107.42  16.74  773.62  34686.36  
350  23.32  30.59  6118.30  23.10  32.75  6061.19  20.86  399.02  22224.31  
400  24.47  41.50  5384.10  24.34  42.76  5322.92  23.17  207.17  14727.86  
500  26.00  42.93  4252.22  25.88  43.34  4181.86  25.45  53.33  7389.40  
700  27.23  32.79  2907.50  27.12  32.94  2856.23  26.92  21.30  2419.15  
1000  27.94  24.17  1998.64  27.87  24.33  1967.03  27.46  48.52  101.42  
inf  28.55  0  0  28.45  0  0  27.07  119.67  6029.33 
D
(m) 
Condition 1  Condition 2  Condition 3  
H=10m  Drift=0º  H=100m  Drift=0º  H=10m  Drift=5º  
Fx
(T) 
Fy
(T) 
Mz
(Tm) 
Fx
(T) 
Fy
(T) 
Mz
(Tm) 
Fx
(T) 
Fy
(T) 
Mz
(Tm)  
300  23.43  16.01  2477.91  26.48  15.19  2487.09  41.29  131.95  9969.88  
350  20.36  4.39  2173.66  23.56  3.61  2364.72  30.98  77.17  6812.31  
400  18.15  2.74  1951.13  21.44  3.47  1956.84  24.47  48.28  4871.71  
500  15.50  7.61  1599.48  18.94  8.15  1600.95  18.05  19.44  2815.15  
700  13.66  7.16  1121.69  17.19  7.49  1122.71  13.95  2.56  1130.23  
1000  12.64  5.96  780.08  16.23  6.16  781.10  13.02  15.74  87.01  
inf  11.80  1.01  8.71  15.46  0.77  6.44  14.31  36.66  2070.12 
Looking at Figure 30, Figure 31, and Figure 32, unstable behaviour can be expected when navigating at high speeds, leading to the incapability of navigating with waves. On the other hand, still water result might not be that useful since they don’t take into account the effect that there always are small waves that might work as small perturbation that might excite large movements in a SES.
A seakeeping test run has been carried out to simulate the performance of the XR1B when performing a turn of tactical diameter D=700m in TD=55s (Fr=1.6). The manoeuvre is performed in the presence of small waves corresponding to a Pierson Moskowitz spectrum with significant wave height of 0.1 meter and mean period of 9 seconds. Figure 37, Figure 38, and Figure 39 show the heave, roll and pitch motion obtained during the simulation. At the beginning of the turn, waves propagate as following seas and covering an angular sector of 30º.
This section aims at showing the effect of navigating in shallow waters. The XR1B have been simulated navigating at 10m/s (Fr=0.4 approx.) in three different water depths of 2.216 m, 10 m, and 100m. Figure 40, Figure 41, and Figure 42 show the wave pattern for each water depth.
Table 12 compares drag, lifting force and pitching moment induced by the hydrodynamic pressure field over the hulls and by the aircushion.
Depth (m)  2.216  10  100  
Hydrodynamic pressure  Fx (T)  8.12  11.46  8.39 
Fz (T)  83.21  134.60  28.55  
My (Tm)  4625.43  3421.15  1693.75  
Air cushion  Fx (T)  27.02  51.35  26.83 
Fz (T)  1059.36  1059.36  1059.36  
My (Tm)  37851.87  1978.97  37841.67 
As expected, the drag and sinkage are larger at shallow water depths than at infinite or very shallow waters. Regarding the pitching moment, it is dominated by the aircushion, resulting in a big reduction of its pitch moment at shallow waters, leading to the immersion of the bow.
In this section, a set of test runs are performed for the XR1B with rigid seals. Rigid seals behave as planning surfaces when in contact with the free surface. Therefore, a free surface rising in the front face of the seal is expected, as well as a pressure profile with a peak pressure where the free surface impacts the seal. Hydrodynamic pressure integration over the wet surface of the seals will provide the forces and moments induced on the XR1B. It will be showed that these forces and moments induced by rigid seals working as planning surfaces are much bigger than the actual forces and moments obtained in reality, and therefore, rigid seals shouldn’t be assumed to be a good approximation for flexible seals. Moreover, the shape of the rigid seal must be guest beforehand to approximate the final shape of the flexible seals in steady conditions, which is quite unpredictable and depends on many working conditions.
Despite of considering that rigid seals are not a realistic approximation for a real SES, the main reason to carry out these set of tests is to prove the validity of the free surfaceseal interaction algorithm explained in section 2.2.4. This algorithm is a key point in order to solve the free surface and flexible seals interaction, since the pressure field obtained via this algorithm will feed the structural solver that will analyze how the seals will deform, and will feedback the new shape of the seal to the algorithm.
In this section, we will analyze two types of bow seals, called plate and fingers. The former is assumed to be a flat plate when undeflected, while the latter is assumed to have the typical finger shape used in bow seals. In both cases, the seal is assumed to be undeflected (linear) in the dry (upper) side, and deflected in the wet (lower) side. The deflection in the wet side is assumed to be cubic and to fulfil tangency condition with the undeflected side, as well as horizontal tangency at the lower end. This adhoc deflection of the bow seal has been adopted based on [7].
However, this model required of a guess of when the transition from the wet to the dry side will happen. As shown in Figure 43: the water level has been assumed to rise up to 1 meter in height. It is shown next that this estimation is quite close for high Froude numbers, when the seal is behaving as a planning surface.
Figure 45 shows the dimensionless pressure distribution obtained in a mid section of the flat plate bow seal. For the lower Froude number, Fr=0.2, the water impact the seal around x=68, which corresponds to a water raise of about 0.8 meters (based on Figure 45). For Fr=0.4, the water impact moves forward up to 68.7 m, corresponding to a water raise over 1.5 meters. On the other hand, as the Froude number increases, the water impact stays very stable around x=68.25 m, corresponding to a water rise of 1 meter, as assumed beforehand. Also, for the higher Froude numbers, the dimensionless pressure distribution corresponds to a typical pressure distribution of a planning surface, reaching up to one right after the impact point, and decreasing to the reference free surface pressure towards the trailing edge.
Towing tests in a range of Froude numbers have been simulated considering including the above mentioned rigid seals. The free surfaceseal interaction algorithm has been used to simulate the planning surface behaviour of the seals. As a result, the free surface elevation and pressure field underneath the seals are obtained. Hydrodynamic pressure integration over the seals surfaces provide the forces and moments induced on the XR1B. Table 13 shows the results obtained in the simulation.
Figure 46 shows the geometry of the XR1B along with the rigid seal geometries used in the simulations. Figure 47 shows the pressure distribution resulting over the bow seal when the ship is moving with constant forward speed and restrained in movement.
Fr  Fx (N)  Fz (N)  My (Nm) 
0  0  0  0 
0.1  9.15E+03  2.23E+04  1.13E+05 
0.2  8.14E+04  1.64E+05  3.42E+06 
0.3  1.97E+05  1.99E+05  1.38E+07 
0.4  5.84E+05  5.27E+05  3.71E+07 
0.6  1.52E+06  1.44E+06  1.00E+08 
0.8  2.56E+06  2.52E+06  1.75E+08 
1.2  5.68E+06  5.64E+06  3.92E+08 
1.6  9.84E+06  9.92E+06  6.87E+08 
2  1.55E+07  1.54E+07  1.08E+09 
Comparing the values in Table 13 with Table 4 and Table 5, it is observed that loads induced by rigid seals working as planning surfaces are much higher that those induced by viscous friction. Moreover, these loads are much higher than any other component as speed increases. Therefore, the use of rigid seals seems not to be appropriate for modelling flexible seals. Figure 47 shows the pressure field induced on the finger of the bow seal located in the port side. Negative pressures indicate areas where pressure drops below atmospheric pressure. The free surface at those nodes cannot detach unless it is surrounded by, at least, one detached node, allowing the entrance of air.
Fr=0.1  Fr=0.2 
Fr=0.3  Fr=0.4 
Fr=0.6  Fr=0.8 
Fr=1.2  Fr=1.6 
Fr=2 
It is important to understand that the high pressure induced in the seal should reshape the seal such that the pressure should drop to a value where dynamic equilibrium exists. Should this equilibrium be reached, the final shape of the seal would be found, and this shape should be parallel to the free surface. Then, based on the amount of wet surface of the seal and its final position, viscous friction could be estimated.
In this section, the dynamics of the XR1B SES is analyzed when in the presence of a regular wave and using rigid seals. While viscous drag, aerodynamic profile drag, and propulsion forces and moments are estimated using the previously used analytical models, forces and moments induced by the hydrodynamic pressure on the hull, aircushion pressure, and hydrodynamic pressure on the seals are calculated in the hydrodynamic solver and accounted for in the dynamic solver in each time step.
Since rigid seals behave as planning surfaces, they have a very large influence on the dynamics of the SES, and the higher the speed of the ship, the larger the influence. Therefore, the dynamics of a SES with rigid seals is not an appropriate model for a real SES with flexible skirts. However, it is a previous step aiming to analyze the free surfaceseal interaction scheme proposed in this work.
Figure 48, Figure 49, and Figure 50 show the XR1B motions for one hundred seconds when moving with forward speed in the presence of a regular wave. It is noticeable that as the speed increases, the movements of the ship become larger. These large movements are induced by the large pressure values obtained on the wet side of the rigid seals. These large pressure values are the result of the high speed of the craft in combination with the seals working as planning surfaces, which stagnates the flow around the impact area and increases the pressure up to the stagnation value.
Figure 51, Figure 52, and Figure 53 show the temporal variation of the forces and moments resulting from pressure integration over the rigid seals. Figure 54, Figure 55, and Figure 56 show snapshots of free surface elevation and pressure over seals.
Figure 48: Heave and pitch motions for XR1B SES with rigid seas at Fr=0.4 and in the presence of a regular wave with H=0.5m and T=6s.
Figure 49: Heave and pitch motions for XR1B SES with rigid seas at Fr=0.8 and in the presence of a regular wave with H=0.5m and T=6s.
Figure 50: Heave and pitch motions for XR1B SES with rigid seas at Fr=1.6 and in the presence of a regular wave with H=0.5m and T=6s.
Figure 54: Snapshots of XR1B with rigid seals at Fr=0.4. From left to right: time=5s, time=10s, and time=15s. From top to bottom: free surface elevation, pressure over bow seal, pressure over stern seal.
Figure 55 : Snapshots of XR1B with rigid seals at Fr=0.8. From left to right: time=5s, time=10s, and time=15s. From top to bottom: free surface elevation, pressure over bow seal, pressure over stern seal.
Figure 56 : Snapshots of XR1B with rigid seals at Fr=1.6. From left to right: time=5s, time=10s, and time=15s. From top to bottom: free surface elevation, pressure over bow seal, pressure over stern seal.
In this section, the behaviour of the seals of the XR1B is analysed. Different types of simulations have been performed. On one hand, the effect of the flexible bow seals has been analysed. Afterwards bow and stern seals have been included in the analyses. The cushion model using in these tests is that presented in previous sections.
The fluidstructure interaction model used for these analyses has been described in sections 2.2.4 and 2.2.5. The ingredients of this model are the following:
The seals structural models used in these analysis consisted on 12976 triangular membrane elements (6864 nodes) for the bow seal and 18608 elements (9428 nodes) for the stern seals.
The following pictures show the boundary conditions using in the structural model, for the bow seals. These include fixed constraints in upper part and nonlinear elastic constraints in the sides of the seals. The later simulate the mechanic contacts of neighbouring fingers.
Figure 57. Fix constraints (left) and nonlinear elastic constraints in bow seals (right). 
The boundary conditions for the stern seal include fixed constraints in upper part and rigid boundary constraints. The later simulate the mechanic contact between the upper and lower lobe. The boundary conditions for the bow seal are shown in the following picture.
Figure 58. Fix constraints in stern seals (left) and Rigid boundary constraints in stern seals (right). 
In this section, the towing test of the craft is analysed, including flexible fingers (bow seals). Three velocities are run, Froude numbers 0.4, 0.8 and 1.2.
In the following sequences, the deformed bow seals, wave pattern and cushion free surface deformation for the different analyses are shown. In the different results, it can be seen how the seals interaction creates an unstable free surface deformation pattern in the cushion, for the higher Froude numbers, while the solution is quite steady for the smaller figure. In this latter case, the deformation of the seals is very small, since the dynamic pressure acting on the seals is quite similar to the cushion pressure.
It is important to remark that the seals deformation obtained for the higher speeds qualitatively agrees with the experimental information available. Unfortunately, no quantitative data exists to further validate these results.
Figure 59. Bow seals deformation sequence: t=2, 4, 6, 8, 10, 12 sec. (from left to right and top to bottom) for Fr=0.4. 
Figure 60. Wave pattern and cushion free surface deformation sequence: t=2, 4, 6, 8, 10, 12 sec. (from left to right and top to bottom) for Fr=0.4. 
Figure 61. Bow seals deformation sequence: t=1, 2, 3, 4, 5, 6 sec. (from left to right and top to bottom) for Fr=0.8. 
Figure 62. Wave pattern and cushion free surface deformation sequence: t=1, 2, 3, 4, 5, 6 sec. (from left to right and top to bottom) for Fr=0.8. 
Figure 63. Bow seals deformation sequence: t=1, 2, 3, 4, 5 sec. (from left to right and top to bottom) for Fr=1.2. 
Figure 64. Wave pattern and cushion free surface deformation sequence: t=1, 2, 3, 4, 5 sec. (from left to right and top to bottom) for Fr=1.2. 
In this section, the towing test of the craft is analysed, including flexible bow and stern seals for a Froude number of 1.2.
In the following sequences, the deformed bow and stern seals, wave pattern and cushion free surface deformation are shown. The results show how the bow seals create an unstable free surface deformation pattern in the cushion area, while the stern seal is deformed until it is adapted to a configuration in which it almost slides on the water.
Again, the bow seals deformation obtained qualitatively agrees with the experimental information available. With respect to this, figure 67 shows a comparison between the computed seal deformation and a picture obtained in a test, taken from reference [9].
Figure 65. Bow seals deformation sequence: t=1, 2, 3, 4 sec. (from left to right and top to bottom) for Fr=1.2. 
Figure 66. Stern seals deformation sequence: t=1, 2, 3, 4 sec. (from left to right and top to bottom) for Fr=1.2. 
Figure 67. Stern seals deformation sequence at 2.5 sec and test data obtained of Reference [9]. 
Figure 68. Wave pattern and cushion free surface deformation sequence: t=1, 2, 3, 4 sec. (from left to right and top to bottom) for Fr=1.2. 
Some of the main concerns when carrying out computational fluid dynamics (CFD) simulations are the computational resources and time required. In our case study, the analysis of the seakeeping and manoeuvring capabilities of a surface effect ship is not an exception. On the contrary, it requires to simulate several complex phenomena coupled between them, such as wave making, aircushion dynamics, free surfaceseal interaction, and ship dynamics.
Different approaches with different types of accuracy can be taken to cope with such an analysis. For instance, in [8], the semicoupled two phase URANS CFD ShipIowa V4.5 is sued for evaluating the performance of the XR1B SES. While the outcome of the analyses is outstanding, the CPUtime reported is quite unaffordable for being used systematically during design stages.
In this work, a potential flow approach along with first order Stokes waves and small movement approximation is used. While this approximation is much simpler than using RANS computations, significant outcomes can be obtained as well, allowing to significantly reducing computational time by 2 or 3 orders of magnitude even when computing on a regular desktop or laptop.
The use of potential flow requires the use of semianalytical formulation in order to account for viscous effects, such as viscous drag and skinfriction drag in seals. On the other hand, the number of degrees of freedom is reduced from four, three 3 velocities components and pressure, to just one (the velocity potential). Moreover, having no need of modelling turbulence allows for having a much coarser discretization.
On the other hand, the use of the Stokes approximation for treating the free surface, along with small movements, leads to having fix boundary along the calculations, which reverberates in having high computational savings.
Despite of starting off with lesser computational requirements than more complex approaches like in [8], special attention has been paid to reducing CPUtime. Looking at the computational time of different part of the code, it was realized that most computational effort was spent in the linear solver. Therefore, the main objective towards reducing computational time has been oriented to accelerating the iterative solver.
Some iterative solvers such as preconditioned conjugate gradient and preconditioned biconjugate gradient failed to solve the linear system resulting from the discretization of the governing equations. However, the preconditioned stabilized biconjugate gradient succeeded and therefore, it was the solver used along with an ILU0 preconditioner.
In order to accelerate the iterative solver, a build in house GPU solver library has been used. Despite of the limitation to used just the diagonal preconditioner (no availability of ILU0 like in CPU computing mode), using the GPU results in speeding up calculations in a range of 2 to 5 times, depending on the case study and on the CPU/GPU architectures used.
The GPU solver library has been built on the Nvidia´s platform called CUDA, requiring for Nvidia GPU units to be used. The cublas library has been used for carrying out linear algebra operations such as dot product, while the cusparse library is used for sparse matrix vector multiplication.
All results reported in this section where obtained using GPU computation. The specific hardware used is the Nvidia GTX 280, which is a three years old graphic processing unit that costs around 200$.
Table 14 provides de time step and computational time required for calculating the towing test cases. In these cases, only half model is considered due to symmetry respect to the plane OXZ. Moreover, the XR1B is fixed, and hence there is no need of solving the body dynamics, which implies that just once iteration per time step is needed.
Fr  Time step (s)  Simulation time (s)  Computational time (s)  Real time ratio (s/s) 
0.1  0.012747  300  3535  11.8 
0.2  0.012747  300  3658  12.2 
0.3  0.012747  300  3598  12.0 
0.4  0.012747  300  3602  12.0 
0.6  0.008061  300  4040  13.5 
0.8  0.005700  30  714.2  23.8 
1.2  0.004031  6  272  45.3 
1.6  0.002549  5  331  66.2 
2  0.002549  4  274.4  68.6 
On the other hand, when simulating the turning circle manoeuvre, the full geometry must be considered, which doubles the computational domain requirements. Furthermore, since we are interested in the seakeeping capabilities, the XR1B is free to move in response to all the forces and moments acting on the vessel. Therefore, an iterative procedure is carry out within each time step to ensure the convergence of the ship dynamics. The computational time for simulating the performance of the XR1B when performing a turn during 110 seconds (see section 5.14) was 39160 seconds, which implies a “real time ratio” of 356. Despite of appearing to be a very high “real time ratio”, it is not compare to RANS solver which usually require real time ratios in the order of 10^{4} 10^{5}.
[1] Leonard, B. P. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Computer methods in applied mechanics and engineering 19 (1979) 5998.
[2] P.A. Ubach and E. Oñate, New rotationfree finite element shell triangle accurately using geometrical data. Computer methods in applied mechanics and engineering. 199 (2010) 383391.
[3] R.L. Taylor, E.Oñate and P.A. Ubach, Finite element analysis of membrane structures. Textile composites and inflatable structures. Springer (2005) 4768.
[4] Maki, K. J., Broglia, R., Doctors, L. J., and Mascio, A. Nonlinear wave resistance of a twodimensional pressure patch moving on a free surface. Ocean Engineering. 39 (2012) 62–71
[5] Yun, L., and Bliault, A. Theory and designo of aircushion crafts.Elsevier ButterworthHeinemann, Oxford, UK. 2000.
[6] Newman, J. N., and F.A.P. Poole, F. A. P. Wave resistance of a moving pressure distribution in a canal. Department of the Navy, David Taylor Model Basin.Report 1619. March 1962.
[7] Doctors, L. J., and Chris B. McKesson, C. B..The resistance components of a surface effect ship. Proceedings of The TwentySixth Symposium on Naval Hydrodynamics. Rome, Italy, September 1722, 2006.
[8] Mousaviraad, S. M., Bhushan, S., and Stern, F. CFD Prediction of FreeRunnin SES/ACV Deep and Shallow Water Maneuvering in Calm Water and Waves. MARSIM 2012. Singapore, April 2327, 2012.
[9] Wiggins, A. D., Zalek , S. F., Perlin, M. and Ceccio, S. L. Development of large scale surface effect ship bow seal testing plataform. 11^{th} International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawai. September 2001.
Published on 14/03/18
Accepted on 14/03/18
Submitted on 14/03/18
Licence: CC BYNCSA license