An Open-Source Python-Based Boundary-Element Method Code for the Three-Dimensional, Zero-Froude, Inﬁnite-Depth, Water-Wave Diﬀraction-Radiation Problem

In this paper, a new open-source implementation of the lower-order, 3-D Boundary Element Method (BEM) of solution to the deep-water, zero Froude-number wave-body interaction problem is described. A validation case for OMHyD, the new open-source package, is included, where the outputs are compared to results obtained using the widely used frequency-domain hydrodynamic analysis package ANSYS-AQWA. The theory behind the solution to the diﬀraction-radiation problem is re-visited using the Green function method. The Hess and Smith panel method is then extended to the case of a ﬂoating object using the image-source to impose the wall condition at the free-surface, and a wavy Green function component to account for presence of free-surface waves. An algorithm for computer implementation of the procedure is developed and subsequently implemented in PYTHON. The wavy part of the Green function is determined using a veriﬁed and published FORTRAN code by Teleste and Noblesse, wrapped for PYTHON using the Fortran to Python (F2PY) interface. Results are presented for the various stages of implementation viz. panelling, body in inﬁnite ﬂuid domain, eﬀect of the free-surface, and eﬀect of surface-waves. The hydrodynamic coeﬃcients obtained from this preliminary frequency-domain analysis are shown to be in satisfactory agreement with ANSYS-AQWA results. Conclusions are drawn based on the performance of the code, followed by suggestions for further improvement by including the removal of irregular-frequencies, multi-body interactions, and bottom interference, which are not considered in the present implementation.


Introduction
In offshore applications, there is often a need to simulate the dynamic response of ships and floaters to the environmental loads due to the wind, waves, and currents, and to the operational loads from cranes, station keeping systems, risers, pipe-lay systems etc.
The problem of a ship in waves is commonly referred to as the sea-keeping problem (Fossen, 2011, pp. 11-12), and involves frequency dependent hydrodynamic parameters. Determination of these parameters requires numerical solutions to the wave-body interaction problem. The theory for this is well developed and is implemented in commercial frequency-domain hydrodynamic analysis packages like WAMIT, ANSYS-AQWA etc.
The multiphysical simulation of marine operations can be done in a component-oriented modeling approach where models representing each component of the system are interconnected, as demonstrated in (Viswanathan and Holden, 2019). The interconnection of the models can be done in a co-simulation arrangement, or as one integrated model. The interconnections can be implemented using either commercial MODEL-ICA environments, or the open-source OpenModelica environment.
Commonly used hydrodynamic analysis packages like WAMIT and ANSYS-AQWA do not render themselves to the open-source philosophy of OpenModelica. This makes the use of component models requiring inputs from such software cumbersome. This is especially true if there is a need to recalculate the hydrodynamic parameters, e.g. due to a large change in the mean position of the wetted surface, based on which the hydrodynamic parameters are calculated. A few open-source, frequency domain, hydrodynamic analysis packages like NeMOH exist, but with limitations in the number of nodes, documentation, and simulation times as discussed in (Penalba et al., 2017).
On this note, we decided to develop an opensource software package for frequency domain hydrodynamic analysis, which makes it possible to calculate the required hydrodynamic parameters for use in a component-oriented modeling-and-simulation package based on OpenModelica. The new hydrodynamic analysis package, OMHyD, is a three-dimensional, lowerorder boundary-element method program, utilizing the source formulation, and the zero-speed, infinite-depth, free-surface Green function. This paper presents the theory and implementation details behind the development of OMHyD. The program package is implemented in PYTHON, and the developed code is made public at github.com/Savin-Viswanathan/OMHyD-PA.
The paper includes a validation of the program package for a cuboidal body where the parameters from OMHyD are compared to the parameters generated for the same body in ANSYS-AQWA.
Lagrange introduced the concept of the velocity potential in the 1780's. Lamb (1879, Ch. 3) applied Green's second identity to express the velocity potential, in the case of a body in an infinite fluid domain, as the effect of a distribution of simple and double sources over the body boundary. John (1949) obtained expressions for the wave function in the case of regular incident waves interacting with bodies satisfying certain geometric assumptions. Hess and Smith (1967) pioneered a method to calculate the flow around arbitrary, non-lifting, 3D bodies, in an infinite fluid domain, based on simple source distributions on flat quadrilateral panels approximating the body surface, also referred to as the panel-method. Hess and Wilcox (1969) reported the progresses in the extension of the panel-method in evaluating the velocity potential by a small oscillating body in the presence of a free surface. Wehausen (1971) formulated a potential flow method for wave-induced motions of free-floating bodies, and demonstrated the agreement of numerical calculations based on the diffraction-radiation formulation with experimental results. Newman and Sclavounos (1988) presented the WAMIT software package, which was based on the panel method, for analysis of water-wave radiation and diffraction.
Further, (Newman, 1977, Ch. 6) and (Faltinsen, 1990, Ch. 3) discussed the wave-body interaction problem within the confines of potential theory. The Green's theorem and the distribution of singularities is discussed in (Newman, 1977, Ch. 4). Garrison (1978) described the use of the Green function to represent the source potential in calculating wave-loads on large floating-bodies. Numerical implementation of the panel-method was described in (Faltinsen, 1990, Ch. 4). Telste and Noblesse (1986) presented a method for numerically evaluating the free surface Green function and its gradient, in deep-water, and zero forward-speed conditions, and also made available the print-out of the associated FORTRAN subroutine in their appendix. McTaggart (2002) described a method for computing three-dimensional hydrodynamic coefficients in the frequency domain, using the Green function given in (Telste and Noblesse, 1986). Guha (2012) described the development of a three-dimensional panel-method using the Green function, as described by Telste and Noblesse (1986), and based on the overall approach by McTaggart (2002). The software implementations of McTaggart (2002) and Guha (2012) have not been made publicly available.
In implementing the PYTHON code for OMHyD, we follow the general framework presented in (McTaggart, 2002) and (Guha, 2012).
The paper proceeds with a brief description of the theory behind the wave-body interaction problem, and the Green function method. A discussion on the extension of the Hess and Smith panel method to the case of an object in the presence of free-surface waves is then presented, along with analytical expressions for the source potential, and its derivatives. The method of determination of the wavy part of the Green function, and its derivatives is also discussed. Subsequently, the algorithm for OMHyD is described, and results of the implementation discussed thereafter.

Theory
In the discussions that follow, x represents a vector, andx represents a unit vector.u represents the time derivative of u.

The Wave-Body Interaction Problem
Considering the physical effects of a floating object interacting with incident waves, we observe the following (Wehausen, 1971): 1. The fluid pressure force on the body changes as a result of the incident waves, and this causes the body to oscillate about its mean position.
2. The presence of the body scatters the incident wave in all directions, and these scattered waves exert associated fluid pressure forces on the body, affecting its motion.
3. The motion of the body generates radiation waves moving out in all directions, and these waves inturn exert associated fluid pressure forces on the body.
Within assumptions of potential flow and linearized boundary conditions, considering the steady-state interaction of a floating object with a regular-wave propagating in deep water, on an infinite free-surface, the total velocity potential in the vicinity of the floating object in regular waves can be expressed as (Mc-Taggart, 2002), (Wehausen, 1971) and (Faltinsen and Michelsen, 1975): Here, ϕ 0 is the spatial component of incident wave potential, ϕ 7 is the scattered wave potential, ϕ j=1,...,6 represent the radiation potential for unit amplitude oscillation in the j th DoF (mode) of the body, and η j is the amplitude of oscillation for mode j. i is the imaginary unit. x is the position vector xî + y + zk of any field point, and γ is the angle of incidence of the wave with respect to the positive direction of the x axis of a right-handed coordinate system with its origin lying on the water plane directly above the centre of gravity of the floating object. The notation is illustrated in Figure 1.
For deep-water, the spatial component of the incident wave potential is expressed as ϕ 0 ( x, γ, ω, ζ a ) = i gζ a ω e kz e −ik(x cos γ−y sin γ) . (2) where, g [m/s 2 ] is the acceleration due to gravity, and h [m] is the water depth. Further, ϕ j , j ∈ {1, . . . , 7} must satisfy: In addition, ϕ j η j , j ∈ {1, . . . , 6} and ϕ 7 must satisfy the linearized kinematic body-boundary condition Here, ∂ ∂n is the normal derivative in the direction of the unit outward normaln to the submerged surface of the body in its mean position, S 0 . The scalars n j are the components of the generalized normal vectorŝ n = (n 1 , n 2 , n 3 ), x ×n = (n 4 , n 5 , n 6 ).
For further details see (Wehausen, 1971), (Faltinsen and Michelsen, 1975), and (Milgram, 2003, pp. 238, 241, and 266 ). It is noted that waves originating as a consequence of incident waves being scattered by the body surface, and those originating due to the motions of the body, must be outgoing, and should have proper amplitude behavior at infinity. This imposes an additional radiation condition on ϕ j , j ∈ {1, . . . , 7}. Details can be found in (Wehausen, 1971) and (John, 1949).
In determining ϕ j=1,...,7 , we make use of the Green function method as described in Sec. 2.2.

The Bernoulli equation
The Bernoulli equation relates the fluid pressure and the velocity potential in the fluid domain.
From the linearized form of the Bernoulli equation, (Faltinsen, 1990, (2.4)), we get With reference to (1), the total velocity potential Φ is expressed as the sum of the diffraction potential Φ D , and the radiation potential Φ R . Thus, where, Consider the wetted body surface S( ξ), ξ = (ξ, η, ζ).
Here the ξηζ frame is the body coordinate frame, which coincides with the global xyz coordinate frame. At the mean position, let the wetted body surface be denoted by S 0 ( ξ). To determine the linear terms of the pressure loads, the hydrostatic pressure should be integrated over the instantaneous position of the body, while the hydrodynamic pressure should be integrated along the mean wetted surface. Thus, from (10) Hydrostatic pressure load , k = 1, . . . , 6. (14) Here n k is defined in (8) and (9). The hydrostatic pressure loads give the mean buoyancy. The component of the hydrodynamic pressure load due to the diffraction potential gives the waveexcitation loads and that due to the radiation potential gives the wave-radiation loads on the oscillatory body.

The diffraction problem: The Froude-Kryloff and diffraction loads
With reference to (12), the diffraction potential at the body surface is written as With reference to the first term on the R.H.S. of (14), the wave-excitation load along the k th DoF is given by Here, the loads due to the incident wave potential Φ 0 are commonly referred to as the Froude-Kryloff loads, and those associated with the scattered wave potential Φ 7 , are referred to as the diffraction loads (Faltinsen, 1990, (3.36)).

The radiation problem: Added-mass and damping
With reference to (13), the radiation potential at the body surface is written as where η j is the amplitude of motion in mode j, and ϕ j is the spatial component of the complex potential due to unit amplitude of motion in the j th DoF. With reference to the first term on the R.H.S. of (14), the wave-radiation load along the k th DoF is given by The wave-radiation load along the k th DoF due to body oscillation along the j th DoF is now expressed as Here, R and I represent the real and imaginary components, respectively, and the dot and double-dot represent time derivatives. The load component proportional to the acceleration is referred to as the added mass, and that proportional to the velocity is referred to as damping. Thus, where A kj is referred to as the added-mass, and B kj is referred to as the radiation-damping, in the k th direction due to wave radiation in the j th mode of oscillation. For details, see (Faltinsen and Michelsen, 1975), and (Guha et al., 2016).

The free-space Green function
Consider a body placed inside an unbounded fluid domain. The velocity potential at any field point Q in the fluid domain is given by (Lamb, 1879, p. 60) i.e., as the resultant potential due to the distribution of simple sources of strength ∂ϕ ∂n − ∂ϕ ∂n on the body surface, and where ϕ and ϕ correspond to potentials outside and inside the closed body surface S, respectively. Here, r is the distance between the field point Q( x), and the source located at P ( ξ), on the body surface, where x = (x, y, z) is the position vector of the field point, and ξ = (ξ, η, ζ) is the position vector of the source point on the body surface. The next step is to define the free-space Green function G * ( x, ξ) = −1 4πr , which is analogous to the sourcepotential, and the source density σ( ξ) = ∂ϕ ∂n − ∂ϕ ∂n . We may then express the potential at any field point as

The free-surface Green function
In cases where the body may move in a fluid domain bounded by other boundaries such as the freesurface, the fluid bottom, or lateral boundaries, additional boundary conditions are imposed on the problem.
The Green function should satisfy the Laplace equation in the fluid domain, and all boundary conditions satisfied by ϕ, except the body boundary condition.
Determination of such a Green function provides an explicit solution to the potential. However, this type of source potential is not known except for some very simple body geometries (Newman, 1977, pp. 138-137). (Wehausen and Laitone, 1960) gives expressions for the finite-depth Green function. Kim (2008) gives a detailed derivation for the deep-water Green function. (Telste and Noblesse, 1986) gives expressions for the infinite-depth free-surface Green function and its derivatives, with the velocity potential expressed as R{G( x, ξ; f )e −iωt }. These expressions, modified for a temporal component of e iωt of the velocity potential, are given as where,  Considering the R.H.S. of (24), and the definition sketch in Figure 2, we make the interpretations that the first term represents the sum of the source potential at a field point x = (x, y, z) due to a unit source at ξ = (ξ, η, ζ), and the image-source, having coordinates (ξ, η, −ζ). This image source (Newman, 1977, p. 160) accounts for the interaction between the body and the linearized free surface which behaves like a rigid boundary. The second term G( x, ξ, f ) is the oscillating potential at the field point due to the oscillation of the source strength, to account for the oscillating flow due to the presence of waves, and therefore referred to as the wavy part of the Green function.
The wavy part of the Green function given by (25) has partial derivatives: where, Here, J 0 (h) and J 1 (h) are usual Bessel functions of the first kind, and R 0 (h, v) and R 1 (h, v) are real functions, which are evaluated by the FORTRAN subroutine GRADIF, the printout of which, is available in the appendix of (Telste and Noblesse, 1986). Numerous works cite this method and it has been validated in (Chakrabarti, 2001).

Determination of the source density distribution
Once the Green function is determined, the solution for ϕ j , j ∈ {1, . . . , 7} can be expressed as given in (Faltinsen and Michelsen, 1975) by where G is the deep-water Green Function given by (24). The application of the body boundary condition given by (7) to the above gives Here, the integral over the body surface is to be interpreted in the Cauchy principal value sense, since it is singular when x = ξ. The limiting process of the Cauchy principal-value gives (Malenica and Chen, 1998): Thus, (38) gives From this the source density distributions associated with the diffraction and radiation problems can be determined, provided we have knowledge about the incident potential, the generalized normals to the body surface, and the normal derivative of the associated Green function. The numerical procedure for the determination of the source density distribution is presented in the following section.

Numerical Solution to the
Diffraction-Radiation Problem

The Hess and Smith panel-method
We refer to (Hess and Smith, 1967), and consider a body of arbitrary shape in an infinite fluid domain, where the body surface is approximated as being built up of N flat panels, each with an associated constant source density distribution σ p . This would allow for the representation of the field point velocity potential as where, ξ p is the position vector of the centroid of the p th panel with surface S p , and G * p is the associated free-space Green function.

Extending the Hess and Smith panel-method
to the free-surface problem Extending the panelling philosophy to the case of a body in the presence of a free surface with waves, the free-space Green function, is now replaced by the deep-water, free-surface Green function given by (24), where, r and r represent the distance of the field point from the source panel and the image panel, as shown in Figure 3. The field point velocity potential is then given by and the component velocities are given by In view of (24), the double integrals in (42)-(45) can be written as Applying the body-boundary condition given by (38) to the body-surface, now approximated by N panels, each having respective constant source-densities, for each mode j of the body, where j ∈ {1, . . . , 6} corresponds to the six radiation modes, and j = 7 corresponds to the diffraction mode, Thus, the Hess and Smith panel-method transforms the body-integral equation (38) to a set of linear algebraic equations, which can be solved to determine the source densities σ jp associated with each mode j, and for each panel p.
The integral of the normal derivative of the Green function in (50) can be determined from the derivatives along the x, y, z directions as given by (63)-(65), and the velocity potential φ j associated with each mode can be determined from equations of the form of (42). Once the velocity potentials associated with each mode is known, the hydrodynamic parameters can be determined.
In determining the integrals of the Green function, and its derivatives, at any field point x, due to a constant source distribution over a flat panel p with centroid at ξ p , as indicated in (46)-(49), we need to solve integrals of the form Sp 1 r dS, Sp GdS, and their derivatives. The method for evaluation of such integrals follow.

Analytical expressions for the integral of the
source-potential, and its derivatives (Katz and Plotkin, 2001, pp. 245-247) refer to (Hess and Smith, 1967) and give analytical expressions for the integrals of the source-potential and its derivatives, as given by equations (62)-(65). The expressions are for a unit source density distribution on a rectangular panel with vertices V k = (x k , y k , 0), k ∈ {1, 2, 3, 4}, and evaluated at a field point Q(x, y, z), with respect to a panel co-ordinate system (PCS) with its origin at the centroid of the panel as shown in Figure 4.
Q(x,y,z) r r 1 r 3 r 4 dS r 2 Figure 4: The panel co-ordinate system.
The following are defined where k = 1, 2, 3, 4. Now, We note that direct use of the expression for S 1 r dS as given in (Katz and Plotkin, 2001, pp. 245-246) appears to be for cases where the field point z-coordinate, in the PCS, is always positive. Hence, we modify the term z in (62), based on discussions in Appendix A, to suit our methodology where the field point may have a negative z co-ordinate.
It is noted that the values given by (63) and (64) become infinite when the field point is on the element edges, and go to zero when the field point is at the centroid (Katz and Plotkin, 2001, p. 247). This is not a problem in the present development, as the computation at the singular field points are not used.
The value given by (65) is singular at the centroid and tends to ∓2π, depending on the direction of approach towards z = 0. I.e., the normal velocity induced by a source panel at its centroid is +σ 2 on the positive side, and −σ 2 on the negative side. This accounts for the limiting procedure described in the formulation of (40).
We avoid this singularity by considering the source panels s to lie at an infinitesimal distance below the body panels p as shown in Figure 5 (Guha, 2012). The point where the body boundary condition is to be satisfied now lies at the centroid of the body panel, at a very small distance from the source panel, and thus the singularity is avoided in the calculations based on (65).

Actual body contour
Edges of flat body panels approximating the body surface Source panels lying below body panels body panel null point collocation distance Figure 5: Source panels beneath body panels.

Numerical evaluation of the integrals of the wavy Green function, and its derivatives
The wavy part of the Green function, G p , and its derivatives along the x, y, and z directions, can be evaluated using the GRADIF subroutine of Telste and Noblesse (1986), as discussed in Sec. 2.2.2. To implement the GRADIF subroutine in our PYTHON code, we reconstruct the subroutine in FORTRAN and wrap it to PYTHON as discussed in Sec. 3.2 and 3.2.4. In evaluating the integrals of the wavy part of the Green function, and of its derivatives, we make use of the consideration that these terms are regular throughout the fluid domain and vary spatially with the wave length, which is generally large compared to the dimension of the immersed surface panel. Hence, G and ∂ G/∂n can be considered constant over a panel, and a valid approximation to the integral is to evaluate the integrands at the centroid of a panel and multiply it by the associated panel area. See (McTaggart, 2002, Sec. 5.3), and Guha and Falzarano (2013). Thus, where ∆S p is the surface area of panel p.

Determination of the source densities -The α matrix
When the expressions for the Green function have been computed, the next step is to compute the source densities. To this end, let α be an N × N matrix where each term α ps represents the induced normal velocity at the centroid of body panel p, p ∈ {1, . . . , N }, due to unit source density distributions on source panel s, s ∈ {1, . . . , N }. In the presence of a free surface, this term would also include the induced normal velocity at the body panel centroids due to the image of the source panels about the free surface. In addition, if there is a wave present on the free surface, then this term would also include the integral of the normal derivative of the wavy part of the Green function due to an oscillating source of unit magnitude on the source panel s, oscillating at a frequency equal to the wave frequency ω.
Thus, (68) is the integral equation (40) expressed as a set of algebraic equations in the matrix form. The source density distribution associated with each panel s for the j th mode, corresponding to the frequency ω, σ ω js , may now be determined from (68).

Determination of the velocity potentials -The β matrix
Next, the velocity potentials are computed. Let β be an N × N matrix where each term β ω ps represents the velocity potential at the centroid of panel p, p ∈ {1, . . . , N }, due to unit source distributions on each panel s, s ∈ {1, . . . , N }, subject to the same conditions as above for the presence of a free surface, and waves on the free surface. Then, we have We note that (69) is the matrix equivalent of the set of algebraic equations representing the integral equation (37) when evaluated at the panel centroids. The solution of (69) gives the velocity potentials ϕ ω jp for j ∈ 1, . . . , 7, for each incident wave frequency ω, evaluated at the centroid of each panel p.

Determination of the Froude-Kryloff and diffraction loads
Having determined the diffraction potentials ϕ ω 7p from (69), the excitation force along any DoF k, for any incident frequency ω may now be expressed with reference to (16) as Here, n kp is the k th component of the generalized unit normal vector given by (8), and (9), evaluated at the centroid of body panel p, and ∆S p is the area of body panel p.

Determination of the radiation loads
Having determined the radiation potentials ϕ ω jp from (69), the added-mass and damping loads along the k th DoF due to body oscillation along the j th DoF with frequency ω, can now be expressed with reference to (19) as Here, j, k ∈ {1, . . . , 6}, n kp is the k th component of the generalized unit normal vector given by (8), and (9), evaluated at the centroid of body panel p, j is the radiation mode of the body, and ∆S p is the area of body panel p.

The aims
To develop a methodology for computer implementation of the panel method, we list out the desired results from the implementation: • body visualization.
• mean free surface visualization.
• image body visualization in case of the presence of a free surface.
• visualization of panel diagonals and normals to check correct orientation of panels.
• field point velocities along the xy and xz planes for visualization of the effects of the source part, image part, and wavy part of the Green function.
• graphical representation of the source density distribution associated with each panel.
• determine added-mass, damping and excitation forces for specified degrees of freedom.
• output storage as a .txt file with simulation parameters specified for identification.
• graphical representation of specific hydrodynamic results from the results stored in the .txt file, and generation of .csv files to enable presentation of results in a L A T E X document or similar environment.

The implementation
From the flow-chart shown in Figure 6, we note that the implementation has the following structure: SFI OMHyD V0R1 : The file containing the below listed code components.
3D Diffraction.py : The front end PYTHON code where the analysis parameters and options are specified.
diff 3d obj func.py : The PYTHON code containing the various functions used in carrying out the analysis.

Hyd mstr( ):
The master function which interfaces the other miscellaneous functions and solves the various sub-problems. am.txt : The text file that contains the results of the analysis.
Plotter.py : PYTHON code for generating plots of the outputs stored in the generated text file.
plt file.csv : The comma separated value file for easy plotting of hydrodynamic parameters in L A T E X or similar package.
The description of the implementation is given below.
3.2.1. Plotting the body, free surface, image body,determination of panel parameters, and calculation of hydrostatics The environment, body, panel, DoF, and plot parameters are to be specified in the 3D Diffraction.py, which is the front-end of the code. Details of inputs are indicated by comments in the code.
The above parameters are passed on as arguments to the function Hyd mstr.py, which is present inside the file diff 3d obj func.py, and a message is displayed in the shell indicating this transfer.
The unit conversion of the user-defined data to units for use during computations is performed, e.g., the wave direction is specified in degrees by the user, and this is converted to radians for use in the code.
The next step is to check if the type of body specified by the user is defined in the code. An error message is displayed if the body type is not defined, and the program execution is terminated. This step is not shown in the flow-chart.
If the body type is defined, then the next step is to set the plot definitions for displaying the body and the free surface, if any. This is followed by a block of code to mitigate conflicting scenarios, as described by the comments in the code. This step is not shown in the flow chart. The next step is to determine the parameters for plotting the body, image, and free surface as required. It is effected by calling the function cube param. Details of the arguments passed, and values returned can be found in the code comments. To summarize, this function takes the passed arguments and generates parameters for plotting (i) the body, and the complete image body, in the case of a submerged body, in the presence of a free surface, (ii) the body surface, excluding the top surface, and image of all surfaces except the top surface, for a floating body, and (iii) only the body, if there is no free-surface. The next step is the plotting of the body/image/free surface based on the return values from the previous step. This is effected by a call to the plt cube function which contains the PYTHON code for plotting the required surfaces.
Once this stage of execution is reached, a progress message is displayed in the shell window. This is not depicted in the flow chart.  The next step is to determine the panel vertices of the source and image body. The order of storage of the panel vertices is important in the view of determining the diagonals and the surface normals. This is effected by calling the function pan vert cub. The process is easily understood from the code comments. The return values are vectors 'vert', 'verti', 'z bp', and 'z tp', which contains the vertex coordinates of the body panels, the image panels, the z-coordinate of the bottom surface, and the z-coordinate of the top surface, respectively.
The next step is the determination and storage of panel parameters as described below.
If the coordinates of the four vertices of a quadrilateral panel are V i (x i , y i , z i ), i ∈ {1, . . . , 4}, then the diagonals D 1 and D 2 are given as Here,î,,k are unit vectors along x, y and z directions. The normal n and the unit normaln to the panel are then given as where n 1 , n 2 , and n 3 are the x, y and z components of the unit surface normal. Similarly, defining L lying on the surface of the panel as The unit vector along L iŝ Another vector lying on the surface of the panel and perpendicular to bothl andn would then be given by The unit vector along P is given as, The postion vector of the centroid of the panel is where and The panel coordinates (x i , y i , z i ), i ∈ {1, . . . , 4}, can be transformed from the body coordinate system defined by (î,,k) to the local coordinate system defined by (l,p,n) as (x l i , y l i , 0), i ∈ {1, . . . , 4}, using The centroid of the source panel (c x p , c y p , c z p ), lying at a distance specified by parameter c p from the body panel, is given by The panel parameters calculated above for each source and image panel are stored as a vector of vectors for use in later calculations. Details of the storage can be understood from the code.
If the diagonals, normals and centroids of the source panels are required to be plotted, functions plt diag, plt norm, and plt src are called.
The function hydrostatics determines the hydrostatic parameters of the body as described below.
Thus, the volume bounded by N panels is given by Here, n 3,j is the z component of the unit normal to panel j, and z j is the z-coordinate of the centroid of panel j.
The vertical and longitudinal centres of buoyancy are given by Here, x j is the x -coordinate of the centroid of panel j.
The longitudinal centre of flotation is given as where A WP is the waterplane area given by The moment of the waterplane about the x and y axes are given by The hydrostatic stiffness terms are then given by They may be non-dimensionalized as where L nd is the characteristic length. This is similar to the non-dimensionalization in WAMIT (Newman and Lee, 2013, Sec. 3.1).
Once the hydrostatics are determined, they are displayed in the shell output along with a progress message.

The source, and image-source Green functions
The next step is to determine the non-wavy part(s) of the Green function at the body panel centroids due to unit source distribution on the source panels/image panels. This is effected by calling the function non wav green func, which returns the following: The IV matrix: An N × N matrix, each element of which represents the induced velocity at the centroid of panel p, by unit source density distribution on source panel s. Each element iv ps is itself a vector [v x ps , v y ps , v z ps ], which represents the x, y, z components of the induced velocity. In determining the IV and IVI matrices, we make use of (63), (64), and (65).
In determining the DNR and DNRI matrices, we make use of the relation v n ps = v ps ·n p , where v n ps is the induced normal velocity at the centroid of panel p due to source density distribution on panel s. v ps = v x psî + v y ps + v z psk is the induced velocity vector at the centroid of panel p due to source density distribution on panel s, andn p is the unit positive surface normal at the centroid of panel p given by (76).
In determining the R INV and RI INV matrices, we make use of (62).

Body in steady, uniform incident flow
In the presence of an incident flow, assuming the body to be fixed, the problem to be solved is a diffraction problem, and hence the mode j = 7 for (68). Each term α ps of the α matrix represents the induced normal velocity at the centroid of the p th body panel due to unit source density distribution on the s th source panel, and on the s th image-source panel in addition, if there is a free surface . Thus, each term α ps is given as Here, p, s, s ∈ {1, . . . , N }. With reference to the definition of the DNR, and DNRI matrices, the α matrix can be expressed as Once the α matrix is determined, the source density distribution associated with each panel may now be determined as Here, in the sole presence of an incident flow defined by , v n 7p is the negative of the component of the incident flow velocity in the direction of the positive surface normal of the p th panel, as given in the description of (68), and is expressed as Here, v ∞ is specified in the inputs, andn p is determined using (76). Once the density distribution associated with each source panel is determined, the null point velocities at the panel centroids may be determined as These are then summed up with the incident flow at the centroid of each panel and plotted as a quiver plot.

Body in the presence of free-surface waves
If waves are present, and the body is free to move along the specified DoFs, then the diffraction-radiation problem is to be solved. We start by determining the generalized normals given by (8) and (9).
The presence of waves necessitates the evaluation of the integral of the wavy part of the Green function and its derivatives. This is accomplished by calling the function wav green func, which returns, for each wave frequency specified in the input: The SG0 matrix: An N × N matrix, each element of which represents the integral of the wavy part of the Green function at the centroid of panel p due to an oscillating unit density source distribution on panel s, calculated based on (25) and (66).
The SDG0 N matrix: An N × N matrix, each element of which represents the integral of the normal derivative of the wavy part of the Green function at the centroid of panel p due to an oscillating unit density source distribution on panel s, calculated based on (32)-(34), (67), and (105).
In determining the R 0 (h, v) and R 1 (h, v) in (32)-(34), we make use of the gradif function which is originally written in FORTRAN, the printout of which, is given as an appendix to (Telste and Noblesse, 1986). We use F2PY (Peterson, 2005) to wrap the FOR-TRAN code for PYTHON. This wrapping generates the gradif... .pyd and the libgradif... .dll files in Windows OS. Instead of .dll files, .so files are generated by F2PY on Linux/Mac OS. The gradif function may now be called into PYTHON, as any other regular PYTHON function.
With reference to (68), the α matrix corresponding to wave frequency ω is now defined as matrix α ω , where each element α ω ps represents the normal component of the induced velocity at the centroid of body panel p due to 1. unit source density distribution on source panel s, and, 2. unit source density distribution on the image of source panel s, referred to as s , in case of the presence of a free surface, and, 3. oscillating source density distribution of unit magnitude on panel s.
With reference to the definitions of the DNR, DNRI, and the SDG0 N matrices, we see that Again, in determining the normal components from the (x, y, z) components given by (43), (44), (45), and (67), we make use of (105).
Once the α matrix is determined, we may determine the source density distribution σ ω js associated with each source panel s, for a given frequency of oscillation ω, for a given mode of oscillation j from equations of the form of (108), with v n,ω jp 's corresponding to the normal velocity at the panel centroids due to oscillations with unit amplitude corresponding to each mode j, as indicated by (7), (8), and (9). Thus, With reference to (69), the β matrix corresponding to wave frequency ω, is now defined as matrix β ω where each element β ω ps represents the induced velocity potential at the centroid of body panel p due to 1. unit source density distribution on source panel s, and, 2. unit source density distribution on image of source panel s, referred to as s , in case of the presence of a free surface, and, 3. oscillating source density distribution of unit magnitude on panel s.
With reference to the definitions of the R INV, RI INV, and the SG0 matrices, we see that The complex spatial component of the velocity potential at the centroid of body panel p, corresponding to the j th mode, with frequency ω, can now be expressed, with reference to (69), as Once the radiation potentials are known, the corresponding added mass and damping components in the k th DoF, due to body motion along the j th DoF, can be determined by using (71) and (72).
Mode j = 7 in this case, is the diffraction problem with the body held fixed at its mean position in incident waves of specified frequencies.
With reference to (68), Where ϕ ω 0 is the velocity potential given by (2), ξ p (ξ p , η p , ζ p ) is the centroid of submerged body panel p.
The wave induced water-particle velocity at the centroid of the p th panel is indicated by u 0p (u 0p , v 0p , w 0p ), where Here, γ [rad] is the wave direction, and Z a = 1 m is the wave amplitude. See, (Guha, 2012, pp. 31-32). Also, with reference to (105), wheren p is the unit outward normal to the p th body panel given by (76). The α matrix being given by (112), an equation of the form (108) may now be solved to determine the source densities σ ω 7s corresponding to each incident wave frequency ω.
The matrix of diffraction potentials ϕ ω 7p associated with each body panel p, for each incident wave frequency ω, is obtained from a relation similar to (116).
The excitation force along any DoF, k, may now be determined with reference to (16) as The added-mass, damping and excitation force matrices are now written to an output file am.txt, along with identification details, and placed inside the SFI OMHyD folder.
If field point velocities are to be plotted, the function field point vel is called. This function calculates the field point velocities in a similar manner to the calculation of the body panel null point velocities described earlier. The only difference being that the IV and IVI matrices are determined for the field points and source/image panels, in this case.
In plotting the field velocities to illustrate the effect of the radiated waves, only the IV matrix is considered, since the effect of the image source is inherent in the calculation of the source strengths corresponding to the radiation modes.
If source densities associated with each panel are to be displayed, they are plotted as colormapped spheres at the body panel centroids by calling the function src str.
A separate post processing code Plotter.py enables the user to generate required plots of data contained within the .txt file. The code has options to specify the k, j components of the added-mass and damping terms, as given by (20) and (21), to be plotted. The plot for the excitation force in the k th direction is also generated. Options are also available to plot dimensional as well as non dimensional values. In addition to this, the code also generates a .csv file for the data contained in the generated plots to make it readily available for use in L A T E Xdocuments by using the pgfplots package.

Results
We endeavour to discuss the results generated using the OMHyD code, as far as possible, by using the screenshot of the graphical output as obtained. Wherever possible, we compare OMHyD and ANSYS-AQWA results.
The required parameters for the various cases presented are specified in the front-end files identified by the figure number in this paper. Thus, for each case presented, there is a different front-end file available in the download, but they all use the same functions contained inside the file named diff 3d obj func.py. More instructions on how to run the different scenarios are given in the README.txt file of the download.
For better comprehension of some of the scenarios, it is recommended to run the corresponding front-end code to get the three-dimensional graphical output. Figure 7a shows the body panels generated by the code for a cubical body of side 10 m discretized into square panels of side 5 m, while Figure 7b shows the normals and diagonals of the body panels. The indication of the surface normals and diagonals helps us in ascertaining that the vertices are stored in the required order during the function call pan vert cub.

Plotting the body and image panels
If the free surface is present, the effect of the image panels are also to be considered in the determination of the Green function. Hence, for a submerged body, in the presence of a free surface, the plot is as shown in Figure 8a. The image body is represented by the grey dashed lines, and the free surface by the red grid at z = 0. The diagonals and normals are not shown to avoid cluttering the image. However, these may be switched on, if required. In this case the body has dimensions 10 × 10 × 5 m and the panel side is 2.5 m.
If the body is floating, then the top surface of the body, and the image of the top surface need not be plotted as is shown in Figure 8b.

Body in infinite fluid in the presence of a steady uniform flow: Effect of the free-space Green function
Consider a cube of side 10 m, in an unbounded fluid, in the presence of a steady flow along the positive X axis given by v f = [1, 0, 0] m/s. The free space Green function is applicable, and is given by G * = −1 4πr . Consider a field point lying at a distance, greater than the length of the panel side, from the body. The velocity potential of this field point may now be expressed as the sum of the incident velocity potential and the velocity potential due to source density distributions on the quadrilateral panels discretizing the body surface. The inability to plot velocities at field points lying closer to the body surface is due to the fact that the velocity potential increases as one approaches the element edges and goes to infinity at the edge, as indicated by (62)   locities is shown in Figure 9. The velocities of the fluid at grid points lying only in the XY and XZ planes are shown to avoid clutter. We observe that the presence of the body causes a change in the flow field, and it is no longer uniform. The flow diverges as it approaches the body and converges as it leaves the body.
The effect of a diagonal flow given by v f = [0.707, 0.707 , 0] m/s is shown in Figure 10.
The source strengths may also be indicated by a color bar as shown in Figure 11. We notice that, on the aft panels facing the incident flow, the source strengths are positive since an outward flow from the panel is required to oppose the incident flow and bring the total velocity of the fluid to zero at the panel null points, thus satisfying the boundary condition. Similarly on the forward side, the body boundary condition implies the presence of a flow into the panels such that the total fluid velocity at the null points is zero. This flow into the panel implies the presence of a sink. This combination of sources and sinks ensures the fulfillment of the conservation of mass equations in the fluid domain. The source strengths are negative on the starboard, port, bottom, and top panels, bordering the aft panels. This is to prevent the flow from separating from the body surface due to the effect of the induced velocity by the aft source panels. Similarly, the positive source strengths on the starboard, port, bottom and top panels, bordering the forward panels prevent the fluid from penetrating the body under the influence of the aft sink panels.

Body in semi-infinite fluid domain:
Effect of the image-source Green function Consider a cuboidal body of dimensions 15 × 10 × 10 m, submerged such that the top surface is 2.5 m below the free surface. Consider a steady, uniform flow along the positive X axis given by v f = [1, 0, 0] m/s in the fluid domain. Since a free-surface is present, the effect of the image-body is also to be considered, and the free-surface Green function is now given by G = −1 4π 1 r + 1 r . The field and null point velocities are shown in Figure 12a.
When the immersion is such that top surface is 0.5 m below the water surface, the field point and null point velocities are as shown in Figure 12b.
When the top surface immersion depth is 2.5 m, the incident flow deviates almost symmetrically. However, when the immersion depth is 0.5 m, the flow deviation is not symmetric. This can be observed from the null point velocities of the top and bottom starboard panels bordering the aft panels, and the second and third row    This points to the wall condition at the free surface. As the body moves towards the wall, the area available for the flow above the body decreases, and the flow deviates to pass around the body along the other available paths. An increase in the fluid velocity in the space between the free surface and the top surface of the body is also observed. The function of the image body is to enforce this wall condition as described in Sec. 2.2. If it were not for the image-source, the whole of the free surface would have to be modeled with source panels, as is done for the body. Figure 13 shows the field and null point velocities when the body is floating at a draft of 5 m. Again, we observe that the flow near the surface does not deviate in the XZ plane but deviates in the XY plane to pass around the body.

Body in semi-infinite fluid domain in
the presence of a free surface and waves: Effect of the wavy part of the Green function We consider a cuboidal body of dimensions 15×10×10 m floating at a draft of 5 m. To illustrate the effect of radiation waves, we consider the induced velocities at the field points, due to the wavy part of the Green function alone, as given by (25). The incident wave frequencies [rad/s] to be considered are specified by the user against the vector 'omega' in the 3D Diffraction.py file. We can select a particular frequency and plot the induced velocity at the field points due to the body oscillating along the required DoF by specifying the requirements against parameters under the n rad freq and the Degrees of freedom fields, respectively. To make the velocities visible, we might need to scale the velocities using the parameter sf rad. The image body is kept hidden, for aesthetic purposes.
Considering the surge motion of the body, the field point velocities are as shown in Figure 14. As the body moves in the positive x-direction, there is an increase in the fluid pressure forward of the body, and a decrease in the fluid pressure aft of the body. The fluid now circulates around the body, from the high pressure zone to the low pressure zone, as shown by the plots for the field point velocities in the XY plane. Considering the XZ plane, we observe that low-and high-pressure zones also develop along the bottom surface.
Considering the sway motion of the body, the field point velocities are shown in Figure 15. Again, we observe similar behaviour of the field velocities as in the surge case.
Considering the heave motion of the body, the field point velocities are shown in Figure 16. We observe    that the figure shows the field point velocities as the body moves down, pushing water away in a radial direction.
Considering the roll motion of the body, the field point velocities are shown in Figure 17. We observe that the figure shows the field point velocities as the body rolls to the starboard. From the XY field point velocities, we observe that as the body rolls to the starboard side, the port side rises, pushing the water away, thus creating a high pressure zone on the port side, while on the starboard side, the motion of the body causes a low pressure zone. From the Y Z field point velocities, we correlate that a high pressure zone is created near the starboard side of the bottom surface as it moves down, and a low pressure zone is created on the port side of the bottom surface as it moves up. These pressure differences causes the fluid flow.
Considering the pitch motion of the body, the field point velocities are shown in Figure 18. As the body pitches forward down, a high pressure is created near the forward part of the bottom surface and a low pressure is created at the aft part as seen in the XZ field point velocities. This causes a flow from the high pressure regions to the low pressure regions. As indicated by the XY plane field point velocities, fluid is sucked down from the free surface at the aft of the body while the fluid is pushed towards the free surface in the forward part.
Considering the yaw motion, the field point velocities are shown in Figure 19. As the body yaws towards the port, a high pressure zone forms near the forward of the port side, while a low pressure zone forms towards the aft of the port side, and vice versa at the starboard side. This sets up the flow corresponding to the yaw oscillation of the body. Only the XY plane field point velocities are shown since the other planes do not show other relevant information.
Considering the heave motion of a submerged body of dimensions 20 × 10 × 5 m, the field point velocities are shown in Figure 20. Here, the body moves up and pushes the fluid on the top while pulling the fluid at the bottom. Looking at Figures 20b and 20c, we observe that this causes a wave crest to form on top of the body. The field velocities oscillate in time, with a frequency corresponding to the radiation frequency of the body, and after half a cycle, the velocities would point in the downward direction with the same magnitude, causing a trough formation above the top surface of the body.
The crest/trough formation in both longitudinal and transverse cross-sections indicate that the wave is three-dimensional.
Another inference drawn is that the wave is a propagating wave and not a stationary one. Had it been a stationary wave, the velocities of the water particles     would have had only the vertical component. Yet another conclusion drawn is that, far away from the body, the wave induced velocities die out.
The above paragraphs, read in conjunction, indicate the presence of circular outgoing waves with reducing amplitudes, as the wave moves away from the body. Thus, we confirm from the results that the wavy part of the Green function generates waves that satisfy the radiation condition as described by the theory in Sec.

2.
These radiation waves cause a change in the hydrodynamic pressure distribution on the body surface, and as a result, the body experiences a force which is expressed mathematically by (19). The component of the force proportional to the acceleration of the body is referred to as the added mass, and the component proportional to the velocity of the body is termed as the damping. Since both velocity and acceleration of the body are frequency-dependent, the radiation addedmass and damping are also frequency-dependent.

The hydrodynamic coefficients of a deeply submerged cube.
Considering a deeply submerged object, we hypothesize: 1. Since the object is submerged beyond half the wave length of the longest wave, the excitation forces must be zero.
2. There exists an associated added mass, which must be independent of the frequency of motion for a deeply submerged body. This is because the extent of the volume of water which is accelerated by the body, as is the interpretation of added mass made by Newman in (Newman, 1977), is constant in the absence of the influence of the free surface.
3. There is no associated potential damping since no radiation waves are present to carry energy away from the body. (Sarpkaya and Isaacson, 1981) gives the analytical expression for the translational added mass of a deeply submerged cube as 0.7ρa 3 [kg]. Here ρ [kg/m 3 ] is the density of water, and a [m] is the side of the cube. Table 1 compares the added masses calculated by the code to results obtained by Guha using MDLHydroD frequency domain solver, as given in (Guha, 2012, p. 60).
The results in OMHyD are for a cube of side 1 m, sea water of density 1025 kg/m 3 , panel size of 0.125 m giving 384 diffracting panels, (or)panel size of 0.0625 m giving 1536 diffracting panels, and an immersion depth of 10000 m. In MDLHydroD, 864 panels were used, while no details are given for the immersion depth and frequencies.
The added mass, damping and excitation forces obtained from the code are plotted in Figure 21 Currently, only square panels can be used in OMHyD. Use of exactly 864 panels gives erroneous results. However, we see that Guha's results with 864 panels lie between our results using 384 and 1536 panels. We also notice that as the mesh becomes finer, the results are diverging from the analytical value, contrary to what one would normally expect.
Sec. 2.2.1 General Modelling Requirements of the ANSYS-AQWA user manual Release 14.5, 2012, requires the surface to be split at the water plane, and hence such an analysis seems implausible in AQWA. Further, there is no indication of the results of a panel convergence study in (Guha, 2012), for this case.
Notwithstanding the above deficiencies, the results presented confirm our hypothesis above. We also arrive at the conclusion that our implementation of the surface panelization, and evaluation of the integrals of the wavy Green function and its derivatives, are not wrong. Further proof of correctness of our procedure is presented in sections that follow.

Panel convergence.
We look at the heave added mass, damping and excitation forces for a cube of side 24 [m], floating at a draft of 12 [m]. Results for discretization with square panels of side 12, 6, 4, 3, and 2 m are shown in Figure 22.
We observe convergence when the panel side is 3 m.

The hydrodynamic coefficients of a cuboidal barge.
Results from MDLHydroD and WAMIT, for a cuboidal barge of dimensions 80 × 20 × 20 m, floating at a draft of 10 m, in sea water with density 1025 kg/m 3 , are given in (Guha, 2012 In OMHyD, the submerged body surface is discretized into 576 square panels of side 2.5 m, and the hydrodynamic coefficients are non-dimensionalized for a characteristic length of 40.0 m. In ANSYS-AQWA, the water depth of 10, 000 m is specified to simulate deep water conditions. ANSYS-AQWA gives dimensional hydrodynamic coefficients, and these are nondimensionalized using (123) and (124), in the Excel sheet available in the download.
The results are compared with results generated for the same barge using ANSYS-AQWA in Figures 23-28.
For higher frequencies, we encounter the irregular frequencies, the effects of which are noticeable in the vicinity of non-dimensional frequency of 2 or higher, in Figures 23-27.
Numerically, irregular frequencies are frequencies at which the α matrix becomes ill conditioned in equations of the form of (68), as described in detail in (Malenica and Chen, 1998). Physically, irregular frequencies correspond to the frequencies of the freesurface standing waves formed in the interior fluid domain bounded by the body surface, (Ohmatsu, 1975, p. 4). As there is no function to eliminate irregular frequencies in OMHyD yet, their appearance is anticipated.
The effect of irregular frequencies can be observed in the results for the same barge using MDLHydroD in (Guha, 2012). Irregular frequency effects can also be observed in analyses carried out using NeMOH, e.g., in (Penalba et al., 2017).
With reference to (Newman and Sclavounos, 1988, p. 4), and (Faltinsen and Michelsen, 1975, p. 6), these       The satisfactory agreement demonstrated here between OMHyD and ANSYS-AQWA results, and also with WAMIT and MDLHydroD, as given in (Guha, 2012), indicates the correct modelling of the wave radiation and diffraction effects, or in other words, the wavy Green function, which was reconstructed from the printout given in (Telste and Noblesse, 1986).

Conclusions
The theory behind the linearized, infinite water-depth, zero forward-speed wave-body interaction problem is refreshed to demonstrate that the determination of the diffraction and radiation potentials enable the evaluation of the first order wave excitation, added-mass, and potential-damping loads.
Lamb's use of the free-space Green function method to represent the field point velocity potential in the infinite fluid domain as the effect of a source distribution on the surface of the submerged body, and its extension to the wave-body interaction problem through the use of the appropriate free-surface Green function, is then presented in brief.
The infinite water-depth, zero Froude-number freesurface Green function formulated in (Telste and Noblesse, 1986) is discussed, and the break-down of the Green function to its source, image-source, and pulsating-source components is presented.
Analytical representations presented for the source, and image-source Green functions, and its derivatives, are based on the lower-order quadrilateral boundary panel method given in (Hess and Smith, 1967), and those for the wavy Green function and its derivatives, are based on the FORTRAN code sourced from (Telste and Noblesse, 1986).
A detailed description of the extension of the Hess and Smith method to the wave-body interaction problem, in order to determine the source density distribution associated with each body panel, from the matrix formulation of the set of algebraic equations representing the linearized kinematic body boundary condition, is presented. Calculation of the panel null-point velocity potentials, for the corresponding frequencies and radiation modes, is then discussed. This is followed by the discussion on the method of calculation of the added-mass and damping loads, for the corresponding frequencies and radiation modes. The method of calculation of the wave excitation loads for the corresponding frequencies along the various DoFs are then discussed.
A detailed description of the implementation of the above method of solution using PYTHON is then pre-sented, and the code for the frequency domain hydrodynamic analysis package OMHyD is made available for public access.
A comprehensive discussion of the outputs from OMHyD is presented, to demonstrate the effects of the source, image-source, and pulsating-source part of the Green function, in order to gain better perspectives on the effect of a submerged/floating body on incident fluid flow, the effect of the free surface, and of radiation loads.
The hydrodynamic analysis results for a cuboidal barge, obtained using OMHyD, is then benchmarked with results obtained using the commercial software ANSYS-AQWA. From other sources discussed, it was noticed that there is satisfactory agreement between OMHyD results, and results from WAMIT, MDLHy-droD, and ANSYS-AQWA.
A shortcoming of the present implementation lies in the presence of irregular frequencies. Such effects are observed in MDLHydroD and NeMOH implementations also. The ability to remove irregular frequency effects is one major feature that sets apart commercial software like ANSYS-AQWA and WAMIT from the other software discussed above.
It is thus concluded that the work succeeded in: 1. Building up a compendium on the source formulation of three dimensional, lower-order boundary element method of solution, for the first order, infinite water-depth, zero-Froude number wave-body interaction problem, using the free-surface Green function.
2. Demonstrating the computer implementation of the theory discussed in the item above, using PYTHON, to develop the open-source frequency domain hydrodynamic analysis software, OMHyD.
3. Demonstrating satisfactory agreement between OMHyD results, and results from other commonly used commercial and proprietary software.
The simplicity of the PYTHON syntax is particularly advantageous for those interested in gaining a better understanding of the working of panel methods. Further, the comprehension of the implementation algorithm would enable code development using compiled software like C or FORTRAN.

Considerations for improvement
The following may be considered for improving the present work: 1. A global approximation to the deep-water freesurface Green function is given in (Wu et al., 2017), and is much simpler than the approximations given in (Telste and Noblesse, 1986). Replacing the wav green func function with the global approximation is expected to reduce the associated computation time.
2. The panel parameters in OMHyD are computed solely based on the vertex coordinates of the quadrilateral panels discretizing the body surface. All subsequent computations are thereafter based on these panel parameters. Noting that triangular panels are special cases of quadrilateral panels, with two coincident vertices, as mentioned in (Newman and Sclavounos, 1988, p. 5), extending the present capabilities of OMHyD to accommodate arbitrary shapes is relatively easy. One may develop a customized panelling procedure, or import panel vertex co-ordinates generated using external software.
3. For objects with a symmetric submerged surface, one can take advantage of the reflection relationship as described in (Hess and Smith, 1967, p. 25), (Hess and Smith, 1962, p. 65), and p. 30 of the lecture notes about sink-source methods and wave induced loads by O.M Faltinsen, prescribed as course material for course MR8300 (Hydrodynamic Aspects of Marine Structures-1), Department of Marine Technology, NTNU. Most commercial software take advantage of symmetry (Newman and Lee, 2013, Sec. 6.1). The advantage is that only the non-redundant portion of the body needs to be discretized, effecting a significant reduction in computation time, as can be inferred from (Newman, 1992).
4. The induced velocity potential at a point lying far from the quadrilateral source panel of area A, of given source density distribution σ, approximates the velocity potential due to a point source of strength Aσ placed at the centroid of the panel, as can be inferred from (Garrison, 1978, p. 106). One may make use of this fact to reduce the associated computational load.
5. Use of compiled languages like FORTRAN or C can improve computation time. OMHyD may entirely be implemented in FORTRAN, wrapped for PYTHON using F2PY (Peterson, 2005), and called into MODELICA taking advantage of interoperability, as described in (OpenModelica, 2020). Or, OMHyD may also be coded in C, and called directly into MODELICA using the interoperability features.
6. Lau and Hearn (1989) present a review of different methods for removal or irregular frequencies. In the 'extended boundary condition method' formulated by Ohmatsu (1975), and Kleinman (1982), a rigid lid is placed on the interior free surface to suppress the internal sloshing modes. Guha (2012) refers to the results of the numerical implementation of this method presented in (Zhu, 1994), including the algorithm for automatic lid generation, and opines that this method is, by far, the most efficient method.
A. Analytical expression for the velocity potential due to a plane source panel Table 2 gives the comparison for S 1 r dS evaluated using the different approximations above for a square panel with side 5 m. We notice that using (126) gives erroneous values for z < 0, as indicated by the values enclosed in the blue box. Hence, we use the corrected expression given by (62) in OMHyD.