Time-Domain Simulation of Lifting Bodies Acting at or Near the Free Surface with Vortex Particle Wakes

Boundary Element Method (BEM) potential-flow solvers are regularly used in industrial applications due to their quick setup and computational time. In aerodynamics, Vortex Particle Methods (VPM) are widely used with BEM potential-flow solvers for modeling lift. However, they are seldom applied to the ocean environment. This paper discusses the implementation of a VPM into Aegir, an existing time-domain, seakeeping, medium-fidelity, BEM potential-flow solver. The wake in the VPM is modeled using both a small dipole buffer wake sheet as well as vortex particles. It has been observed that this method captures both the details of complex wake patterns behind lift-producing surfaces and the expected lift force, thus improving the accuracy of the solution. Two new contributions presented in this paper include the extension of the VPM from previous source-based methods to a potential formulation and full interaction with free surface waves.


INTRODUCTION
Boundary Element Method (BEM) potential-flow solvers require less setup and computational time compared to higher-fidelity solvers while still capturing the most significant physical effects of the simulation. However, common methods for modeling lift in BEM solvers often require the inclusion of a dipole wake sheet which increases the setup difficulty for users. Due to geometric restrictions, such as the need to prevent boundary elements from overlapping or twisting, it is difficult to evolve dipole wake sheets with the flow. Typically, these wakes are treated in a linear approach with the position of the wake prescribed by the user. More critically, the interaction of a wake sheet with a downstream body requires very careful treatment and does not lead to a robust solution. This places restrictions and difficulties on the application of the method. Vortex lattice methods remove some of the geometric complexity, but still have large limitations of their own.
Many of these shortcomings can be alleviated by applying a Vortex Particle Method (VPM). This paper describes the development and implementation of an unsteady VPM to model the lift and induced drag on lifting bodies operating at or near the free surface in Aegir, a time-domain, medium-fidelity, potential-flow seakeeping program that uses an advanced, Non-Rational Uniform B-Spline (NURBS) based, highorder BEM. The VPM implemented in Aegir couples the influence between the vortex particles and the potential-flow boundary elements. This method only requires that the user provide hull and lifting surface geometries, environmental conditions, and run scenario to establish the boundary and initial conditions for the simulation.
The objectives of this paper are to: (1) present the formulation of a VPM, (2) discuss the steps of coupling the VPM with Aegir's high-order boundary elements, (3) present the proper particle advection with an observed wake shape, (4) verify the resulting lift forces against existing theory, and (5) discuss the future work needed to further the capabilities of the VPM, including continued testing of interactions with downstream bodies and the inclusion of the free surface.

METHODS
This new VPM uses a combination of dipole panels (or NURBS-based boundary elements) and vortex particles to model the flow in the wake of a lifting surface, which is comparable to the representation of the wake introduced by Willis [1]. However, there are some key differences between the two approaches, the VPM in this paper: (1) has been strongly coupled to free surface waves, (2) uses a potential-based source and dipole distribution as well as high-order discretization as opposed to the source-only distribution on the body and constant dipole panels in the wake used by Willis, and (3) has body-conforming particles.

Formulation
Implementing the VPM in Aegir required the addition of particle-to-particle, particle-to-body, particle-to-free surface, body-to-particle, and free surface-to-particle influences in Aegir's Boundary Integral Equation (BIE), as depicted in Equation 1, since the body-to-body, body-to-free surface, free surface-to-free surface and free surface-to-body influences already existed [2,3].
The exact methods for including each of these influences will be discussed in the Particle Evolution section of this paper.

Buffer Wake Region
The buffer wake region is a truncated dipole wake sheet attached to the trailing edge of each lifting surface component. The inclusion of a buffer wake region provides the ability to utilize current Kutta solvers in Aegir, a geometric guide for the generation of vortex particles, a simplified method for computing initial particle strength, and a more accurate computation of pressures at the trailing edge of the lifting surface. Furthermore, the user is not required to provide any additional inputs for the creation of the buffer wake region as this buffer dipole wake sheet is automatically generated behind each lifting surface and split at geometric discontinuities in the trailing edge once the simulation is started. Additionally, each buffer wake sheet is automatically discretized to match the number of elements in the span-wise direction of the lifting surface component to which it is attached. The method of particle creation used in this VPM is inspired by literature derived from source-only distribution on the body and constant dipole panels in the buffer wake region, but takes advantage of the potential-based source and dipole distribution as well as high-order discretization in Aegir.
Each buffer wake sheet extends two panels downstream. The total wake length, , for each buffer wake sheet, described in Equation 2, is a function of the free stream velocity, , and the simulation time step, .
Both of the stream-wise panel strips, depicted in Figure 1, have a specific purpose. The first strip of panels, shown in green, is used to compute the potential jump between the top and bottom patches at the trailing edge of the lifting body by satisfying the Kutta condition. Though any Kutta condition may be used with the implemented VPM, Aegir gives the user the option of either solving a linear (Morino) or nonlinear (Iterative Pressure) Kutta condition [7]. The potential distributed on the first panel strip propagates directly to the second panel strip, depicted in blue, in the next time step, and the potential of the first strip is recomputed. The potential on the second panel strip is used to compute the initial strength of the particles created at the following time step. Both of the wake strips have lengths in the stream-wise direction equal to * .

Wake Particles
Once the buffer wake sheets have been generated, the particles are created using the approach described below. After all of the particles are created, they are advected downstream in a twostep process, discussed in the Particle Evolution section: first, the influences, or induced velocities, are calculated, then the particles evolve according to the calculated influences.
Particle strengths are computed and stored in preparation for their creation during the following time step. The initial particle strength is calculated by finding the product of the tangential derivative of the potential at the center of each of the second strip buffer wake panels and the area of the panel, as in Equation 3 below.
The gradient of the potential at any point on the dipole wake can be evaluated exactly due to Aegir's highorder panel method.
As the potential from the previous time step is propagated downstream from the first panel to the second, the potential on the second panel strip in the wake is propagated by creating vortex particles one panel length, U∆ , downstream of the second wake strip. The initial particle strength was calculated from the potential of the second wake strip at the previous step. After the particle is created, it will travel freely for the rest of the simulation by evolving in both strength and position.
Particles are created downstream of the buffer wake at each time step in the simulation by converting the dipole panels on the second panel strip into equivalent vortex particles and then advecting them with the free stream velocity one time step downstream. This process is illustrated in Figure 2.

Canceling Particles
The downstream edge of any dipole wake sheet will result in an effective strong end vortex. For regular dipole wake sheet methods, best practices of wake length convergence have been established so this end vortex propagates sufficiently far downstream as to not cause instabilities. When a dipole wake sheet is in the proximity of the free surface, these best practices suggest that the wake sheet extend beyond the extent of the free surface as to ignore its artificial effects.
For the current approach, the end vortex would have a large influence on the vortex particles as well as the free surface and any nearby bodies. Therefore, the artificial vortex at the end of the buffer wake is canceled once particles are created. This is done by creating stationary particles with strength updated at every time step to be equal and opposite of the end vortex at that time to counteract its effects.
These particles are placed at the midpoint locations on the downstream edge of each panel on the second wake strip and are created by treating the downstream edge of each of the second wake strip panels as separate vortex lines and then collapsing the vortex lines into particles. The strength magnitude of these particles is calculated by multiplying the length, L, of the vortex line the particle represents by the potential evaluated at the particle position.

Velocity Influences
As mentioned in the Formulation section of this paper, there are five particle-related influences that need to be computed. Two of these, the particle-toparticle and -body influences, involve the vorticityinduced velocity, Equation 5, which calculates the influence each particle subjects on all other particles or panels in the simulation [1].
In this equation, ����⃗ is the vorticity-induced velocity acting on point , ����⃗ is the position of point , and ⃗ and ⃗� ⃗ � are the position and strength, respectively, of the particle influencing point .
However, it can be seen in the denominator of Equation 3 that as particles approach other particles or panels, the velocity calculations become increasingly unstable. To combat this issue, a regularized 3D high-order algebraic smoothing function, Equation 6, which is dependent upon the distance between points of evaluation, , is added [4].
In order to add the smoothing function, a core radius, , needs to be assigned to each vortex particle, thus yielding Equation 7. Winckelmans [4] recommends a core radius of 1.3 times the average distance between particles which is taken as 1.3 * * in the presented VPM.
The particle-to-particle influences are found by calculating the sum of the vorticity-induced velocity, Equation 7, in which each particle subjects an influence on all other particles in the simulation [1].
Similarly, the particle-to-body influences are accounted for explicitly by updating the velocity used in the no flux body boundary condition with the sum of the vorticity-induced velocity, �⃗ , each particle imposes upon each of the body panels, as depicted in Equation 8.
The linearized kinematic free surface boundary condition, Equation 9, which is handled explicitly, must be updated to account for the vertical flux the particles induce on the free surface by adding an additional term, , to account for the change in vertical memory potential.
In this equation, ��⃗ is the mean body velocity, is the basis flow potential, is the local, or instantaneous perturbation potential, ψ is the wave memory perturbation potential, and ς is the linearized wave elevation, according to Kring [2,3].
The linearized dynamic free surface boundary condition, Equation 10, which is handled implicitly, is updated by splitting the memory component into two parts: the first for the memory effects due to waves, and the second for those due to the lift generated by the vortex particles.
Lastly, the body-to-particle and free surface-toparticle influences are determined by evaluating a second set of boundary integrals which calculates the source-and dipole-induced velocities that the body and free surface panels subject onto the particles. It should be noted that the buffer wake sheets are treated as part of the body to which they connect.

Particle Evolution
Once all the influences have been calculated, each particle's new velocity, ⃗ � ⃗( ), � , is found by summing the free stream velocity, the vorticityinduced velocity, and the source-and dipole-induced velocities acting on the particle. The position of each particle is then evolved as a function, Equation 11, of their new velocity [1].
Additionally, the strength of each particle must be updated using Equation 11 due to the effects of stretching [1].

RESULTS AND DISCUSSION
Testing can be classified into three categories. The first verifies the calculated lift force, the second reviews the resulting wake patterns, and the third examines the ability of the VPM to handle more complex problems involving a free surface or downstream bodies.

Results
The lift force verification involved performing a preliminary steady lift test using the time domain solver as well as a set of simulations of a NACA0012 cross-section airfoils forced to heave. A forced instantaneous pitch run was used for wake pattern testing in addition to the mentioned heaving simulations. The VPM was used to simulate a steady lifting hydrofoil upstream of a non-lifting hydrofoil as well as a hydrofoil forced to pitch at various depth to chord ratios beneath a free surface to examine its application to more complex problems.

3.1.1
Lift Force Verification This section will review two test cases in which the VPM lift force outputs were verified against theory.

Steady Forcing Lift Force
Aegir's option to solve steady forcing in the body boundary condition was enabled for initial timedomain lift force verification tests. This allowed for a direct comparison to the steady-state solver and aided in identifying any sensitivities in the formulation. The time history of the lift force in these simulations was expected to converge near the steady solution when simulated without incident waves. The validity of the resulting lift force from the VPM for these simulations was analyzed by performing a combined spatial and temporal convergence test to verify that it converged to a solution.
This convergence test was conducted to check that the coefficient of lift approached that predicted by Prandtl's Lifting Line theory [6] as the number of panels used to discretize the foil and wake increased. Aegir simulations involving the pre-existing dipole wake method are easier to control in convergence testing as the panel size on the wake approaches that on the body. With the VPM, performing combined spatial and temporal convergence tests is necessary to set the panel size on the body equal to that on the buffer wake, ∆ . The runs were performed using the new VPM on an aspect ratio 8 hydrofoil with a NACA4412 cross-section (chord of 0.5 ft; span of 4 ft) fixed at a 10° angle of attack and traveling at a forward speed of 1 ft/s using the steady forcing capability in the time-domain solver. The convergence test was performed using time steps 0.125, 0.100, and 0.075 seconds corresponding to panel lengths in the buffer wake and on the foil trailing edge of 0.125, 0.100, and 0.075 ft respectively. The coefficient of lift is observed to converge towards a value near that predicted by Prandtl, which is identified by the black, dashed horizontal line in Figure 4. It is not expected that the solution would converge directly onto the two-dimensional coefficient of lift predicted by Prandtl.

Forced Heaving Motions
To confirm the effect that the particles have on the pressure distribution of their lifting body, a set of forced heaving simulations were performed on airfoils of increasing aspect ratio with a NACA0012 crosssection. The lift from these were compared to Theodorsen's theoretical 2D lift [8]. To do this, the quasi-2D lift was calculated by integrating the pressure along only the root of the foil. The chord length was held fixed for this test matrix at 1.0 m while the span of the foil was increased resulting in lengths of 1, 2, 4, 8, and 12 m (aspect ratios 1, 2, 4, 8, and 12 respectively). The foils were forced in a sinusoidal heaving trajectory of amplitude 0.0159 m and oscillation period of 2 s while being held fix in all five other degrees of freedom.
It can be seen in Figure 5 that the quasi-2D lift converges towards Theodorsen as the aspect ratio of the foil increases (or as the 3D run better mimics a 2D run at that root location).

Instantaneous Pitch Motion
To confirm the wake pattern resulting from the particle advection, the wake behind an aspect ratio 8 hydrofoil with a NACA0012 cross-section (chord of 1.0 ft; span of 8 ft) was forced with an instantaneous pitching amplitude of 5° and traveling with a forward speed of 1 ft/s. The pattern of wake vortices, Figure  8, is captured in the particle wake of the hydrofoil and is compared to the test performed by Willis, Figure 9, [1]. Figure 6. Wake pattern of an instantaneous pitch run from presented VPM using the High-Order Approach Figure 7. Wake pattern of an instantaneous pitch run from Willis [1] The resulting particle wake represents the same volume of fluid as that shown in the results by Willis along with a very similar shape.

Wake Roll-Up
Additional simulations were performed similar to that described above in section 3.1.2.1 of instantaneous pitch. To observe the effect of wake roll up on the body solution, these runs were repeated with various pitch angles 1˚, 3˚, and 5˚ and therefore increasing wake roll-up as shown in Figure 8. The pressure distributions at the root of the foils were computed using Aegir's dipole wake sheet method and compared to those found using the VPM.
As shown in Figure 9, as the magnitude of the angle of attack increases, the pressures computed by the VPM, deviate from the ones calculated using the dipole wake sheet method. Since higher angles of attack result in increased wake roll up, this confirms that wake roll up effects are being captured in the potential distribution solution on the body. Figure 8. Pressure distributions at the root of the hydrofoil of instantaneous 1˚, 3˚, and 5˚ pitch runs from presented VPM compared to the pre-existing dipole wake method

Forced Heave
A clear pattern of wake vortices is captured in the particle wake of the hydrofoil discussed in section 3.1.1.2. It is observed to oscillate in the rotational direction as anticipated, as shown in Figure 8. In order to confirm that the particles are traveling at the expected period, the distance between the particle trajectory troughs is measured and compared to the forced heave wavelength, 2ft, which is found by taking the product of the forced heave period and the forward speed.

Complex Problems
The addition of the free surface and bodies downstream of the lifting body were simulated in order to begin testing the influences and interactions between the particles and these boundaries.

Free Surface Effect
To observe the influence of the particles on the free surface and vice versa, a set of simulations with a NACA0015 cross-section aspect ratio 8 hydrofoil were performed. These hydrofoils underwent a forced sinusoidal pitch motion and were placed at three unique depths corresponding to depth to chord ratios of 1.0, 0.75, and 0.5.
In Figure 9, the results are compared from decreasing depth to chord ratios with the deepest location (a) to the shallowest location (c). The color pattern on the free surface depicts the wave elevations resulting solely from the presence of the foil undergoing 1.0 m/s forward speed. The resulting advected particle patterns are shown through the opacity of the free surface in order to observe their overall trajectories.
It is observed that as the foil gets closer to the free surface, both the resulting elevation on the free surface and the particle wake patterns are more highly influenced. This confirms that the influence calculations between particle and free surface are represented.

Downstream Bodies
Due to the modifications made to the body boundary condition, as particles get closer to the body boundary, the influence from the body enforcing the no flux condition should push the particles away from the body.
A simple test was completed in which the particles were given zero strength to test the influence from body boundaries onto the particles. The particles successfully advected downstream and were diverted to flow around the downstream body. Continued simulations are necessary to observe the relationship between the body and particle influences when the particles are created with strength and advected normally. Some preliminary testing suggested sensitivities between these influences and the time-step of the numerical simulation.

Discussion
The presented results show that the VPM captures comparable lift forces with theory and provides a convergent solution. Additionally, it is observed in section 3.1.2 that realistic flow characteristics are obtained by the particle wakes. The effects of wake rollup and downwash are accounted for, which is a shortcoming of using Aegir's pre-existing dipole wake method.
The VPM is more powerful than the dipole wake method in that it was designed to have no geometric restrictions produced by the requirement of a user defined dipole wake sheet. Consequently, the implemented Vortex Particle Method eliminates error deriving from dipole wake sheets being incorrectly generated, placed or sized.
Though the current focus on the implementation of the described VPM is to enhance the lift modeling calculations and capabilities in an ocean environment, work will continue to improve the computational efficiency of the VPM by implementing a fast solver. Due to the increase in data storage for the particle structures, various approaches for an accelerated solver are being investigated to determine one most compatible with Aegir's code framework including: Fast Multipole Tree Method, and pre-corrected FFT.
With continued work, the VPM's ability will be improved to model complex simulations which include free surface effects as well as bodies traveling in the wake of a hydrofoil such as the one depicted in Figure 8. In the case of this example, the wake particles of the two bow lifting bodies would advect around the two aft lifting bodies freely, improving the artificial shape of the depicted wakes. Particle interactions would contribute more details to the disturbance of the free surface than the dipole wake approach. Later applications will focus on rotating propellers.