Wave Resonance and Dissipation in Collections of Partially Reflecting Vertical Cylinders

Wave fields near coastal features exhibit complicated patterns as incoming waves meet scattered waves to generate areas of resonance with wave amplification and areas of dissipation with reduced amplitude. A mathematical framework to study such problems was developed by MacCamy and Fuchs (1954) for one vertical cylinder that fully reflects plane waves in water of constant depth. Such periodic solutions reduce to the wave equation in two dimensions (Lamb 1879, p. 195):

This may be formulated using the Helmholtz equation of a complex function φ in the reduced form (Courant and Hilbert 1962):

Here, the wave number k is given for waves with specified period T in water depth h by the dispersion relation: ${\omega }^2=kg\text{tanh}kh$ (Airy 1845, p. 290), where $\omega =2\pi /T$ is the angular frequency. Berkhoff (1976) developed this complex function for wave interactions in water of constant depth to study wave amplitude using its modulus $|\phi |$ and wave phase using its argument $\text{arg}(\phi )$. This formulation of the problem is also important for the study of waves in acoustics (Morse and Ingard 1968, p. 728), electromagnetism (Sommerfeld 1972, p. 159), and optics (Goodman 1996, p. 55).

Simple analytic solutions provide insight into coastal interactions, such as the exact solution for plane waves that encounter a single vertical cylinder shown in Fig. 1; this fully reflective solution follows from MacCamy and Fuchs (1954), and the equations necessary to construct this figure are articulated later. Plane waves are shown moving from left to right by depicting the wave field’s amplitude $|\phi |$ (with thicker lines at the scaled amplitude 1 of the plane waves) and the cosine of its phase, $\text{cos}\hspace{0.17em}[\text{arg}(\phi \left)\right]$, with wavelength equal to 4 times the radius of the element, similar to the fully reflective solution shown by Berkhoff (1976, Fig. 3.4). Note that contours illustrate amplitude at intervals of 0.1 times the incoming amplitude except for Fig. 11, where the interval is 0.2. The wave interactions with the coast are incorporated within a coastal boundary condition adapted from Berkhoff (1976):where the reflection coefficient R is real and takes on values from 0 for full absorption by the element to 1 for full reflection. Clearly, the resonance of partially standing waves to the windward side becomes more pronounced for a reflective coast, and dissipation to the leeward side is larger for an absorbing coast. Such solutions inform the design of wave run-up around a circular island (Sarpkaya and Isaacson 1981) and the computation of wave forces (Morison et al. 1950) and breaking wave forces (Swift 1989) on a vertical cylinder.

More advanced problems exist where wave fields are shaped by reflection and absorption in collections of interacting objects. Trapped modes with resonance occurring at particular frequencies are known to occur as plane waves moving past cylinders arranged as a straight line in a waveguide (Linton and McIver 2002) and as a grating of vertical cylinders (McIver et al. 1998). Recent investigations have also examined the propagation of periodic waves through fields of aquatic vegetation (Liu et al. 2015) using vertical cylinders with both regular and nonuniform spacings (Mei et al. 2014). Study of such systems of scatterers requires the capacity to identify the specific geometric configurations through which waves either resonate or dissipate (Platts et al. 2002). Accurate depiction of wave fields is also important for the design and interpretation of monitoring networks (Pandey et al. 2014), where the location and spacing of sensors provide information about an unknown composition of partially reflecting objects.

Various numerical methods have been employed to compute wave fields near an assemblage of elements. Propagation of waves through a cylindrical lattice was examined by Guo and McIver (2011), ignoring interactions among elements for elements much smaller than the wavelength and applied to elements isolated at relatively large distances from each other. The method of matched expansions was used by Kakuno and Liu (1993) and Kakuno and Nakata (1998) to study reflection by a slit-type breakwater formed by an array of vertical cylinders, with geometry similar to the waveguides from acoustics and electromagnetism. Waves near a cylindrical element have also been studied using the finite-difference method (Britt et al. 2013), finite-element method (Calhoun et al. 2008), finite-volume method (Lemoine and Ou 2014), integral representation (Chen et al. 2009a), and the boundary integral equation method (Chen et al. 2013).

Clearly, it would be valuable to achieve solutions with the precise depiction provided by the analytic solution of a single circular element in Fig. 1 for problems where collections of circular elements collectively interact to form more complicated wave fields. This is accomplished here through the adaption of methodology from the analytic element method (AEM), which has been developed over the past 30 years (Strack 1989, 2003) to design solutions for engineering problems in hydrogeology governed by the Laplace, Poisson, heat, and modified Helmholtz equations. This AEM paradigm has evolved the capacity to achieve nearly exact solutions to complicated problems with hundreds to thousands of closely interacting features (Janković et al. 2009; Steward 2015). The AEM is formulated next to study the interactions of waves within collections of objects, representing its first application to engineering problems governed by the wave equation. Results illustrate the resonance and dissipation of waves occurring across variations in wavelength, element size, and reflection coefficients, and the errors associated with matching boundary conditions quantify the achievement of nearly exact solutions. Although this study used coastal features shaped like vertical cylinders in water of uniform depth, it provides a foundation for future extension to problems with more complicated coastal features and bathymetry.

The AEM provides a framework to study wave fields, and the following presentation is organized within the four steps of Steward and Allen (2013): (1) subdivide a domain into elements with prescribed geometries, (2) develop linear combinations of influence functions that satisfy the governing partial differential equations, (3) superimpose elements and their influence functions to achieve a comprehensive solution, and (4) evaluate these functions at control points to solve boundary conditions.

Coastal features are approximated as a collections of I circular objects inserted into a field of plane waves, similar to Sumer and Fredsøe (1998). The geometry of each vertical cylinder is prescribed using the variables in Fig. 2, where the element i is centered at $({\underset{i}{x_c}}_{},{\underset{i}{y_c}}_{})$ with radius ${\underset{i}{r_0}}_{}$, and the underscript i is used to associate the center and radius with a particular circle. A set of M control points is identified where the coastal boundary condition [Eq. (3)] will be applied later. They are evenly spaced about each element at anglesso, the control point m for element i is located at

The influence functions for analytic elements are provided as a summation of the complex function φ for plane waves moving in the x-direction with the linear combination of Bessel functions of r times trigonometric functions of θ from separation of variables (Moon and Spencer 1961). This gives a series solution to the Helmholtz equation [Eq. (2)]:where the wave number k = 2π/L is related to the wavelength L (Goodman 1996); the terms $\stackrel{\text{cos}}{c_n}$ and $\stackrel{\text{sin}}{c_n}$ are complex coefficients; and the Hankel function of the first kind satisfies the Sommerfeld (1972, p. 189) radiation condition. Summation of these functions provides the wave solutions illustrated in Fig. 1 with appropriately chosen coefficients. The gradient of this complex wave function is needed for the derivative term in the boundary condition [Eq. (3)] and is given by the chain rule for partial differentiation:withwhere the partial derivative of influence functions with respect to r is evaluated usingwhich was adapted from Abramowitz and Stegun (1972, p. 361).

A comprehensive solution for problems with wave interactions among many circular elements is obtained through linear superposition of the influence functions for each element:where these functions are evaluated at a distance from the center of element i and the angle about its center:

Note that φ varies as a continuous function and may be computed at any point in the domain to provide the amplitude and phase of the wave field, and this series approximation has been used in coastal engineering problems with one cylinder (MacCamy and Fuchs 1954) and multiple cylinders (Linton and Evans 1990, Eq. 2.9). Now, the AEM separates the functions for element i from those for all other elements:using the additional function from Steward (2015):This partitioning separates the coefficients for element i from the unknowns for all other elements, which is necessary to achieve the iterative solution methods on an element-by-element basis as presented next.

The prescribed boundary condition [Eq. (3)] is matched at control points, ($\underset{i}{x_m},\hspace{0.17em}\underset{i}{y_m}$) in Eq. (5), by adjusting the coefficients ${\underset{i}{\overset{\text{cos}}{c}}}_n$ and ${\underset{i}{\overset{\text{sin}}{c}}}_n$ in Eq. (12a). Various methods exist to match these conditions. For example, the exact solution in Fig. 1 for a single element is given bywhere J_n are Bessel functions of the first kind. This extends published solutions for fully reflecting boundary conditions (Morse 1936; Havelock 1940; MacCamy and Fuchs 1954) to boundaries with partial reflection, with coefficients found by expanding plane waves as an infinite series (Sommerfeld 1972, Eq. 21.2b) and satisfying the boundary condition [Eq. (3)] term by term for each n. A solution for multiple elements (when I > 1) was developed by Linton and Evans (1990) as a system of linear equations relating all coefficients of all elements to all boundary conditions for all elements.

The AEM is formulated here as an iterative method that sequentially finds the coefficients for each element while holding the coefficients for all other elements constant in the additional function [Eq. (12b)]. Although the details to formulate the solution via matrix multiplications are deferred to Appendix I, the equations that provide coefficients for vertical cylinder i [Eq. (26)] are repeated here:using coefficients ${\mathcal{J}}_n$ and ${\mathcal{Y}}_n$ from Eq. (18). Thus, element i is solved by computing its additional function using the coefficients for all other elements, then adjusting the coefficients for element i to match its boundary conditions. And iteration continues this process for elements i = 1, ⋯, I until all elements have been solved for the iterate. Gauss–Seidel iteration is adopted whereby the coefficients for each analytic element are updated after solving, and the solution sequences continuously through elements until convergence occurs, which is defined to occur here when the variation in φ is very small (10⁻¹²) at every control point for every element between successive iterates. The number of coefficients is truncated at N = 20 for all figures to match the accuracy used to depict the exact solution in Fig. 1.

The number of control points is chosen following the overspecification principle of Janković and Barnes (1999) with M = 62 for each circle (so the number of unknowns M is approximately 1.5 times the number of unknowns 2 N + 1). This principle matches conditions approximately at a larger number of control points than the number of unknowns by minimizing a least-squares objective function, and residual errors are calculated for each figure to demonstrate the nearly exact solutions that are achieved across all control points for all elements. Solutions were obtained using Scilab 5.4.0 on Beocat (the Beowulf High-Performance Computing Cluster at Kansas State University) and visualized by contouring values of φ computed across a grid of 401 by 401 points.

Models that accurately satisfy the boundary conditions imposed by coastal features across variations in reflection coefficient provide insight into the capacity of breakwaters and coastal structures with gaps to provide shelter (Dalrymple and Martin 1990). Three sets of model simulation were conducted to study coastal interactions. The first two sets contained a collection of regularly spaced vertical cylinders, such as those studied by Sumer and Fredsøe (1998), and illustrate the impact of wavelength and the impact of gap size versus element size. The third set contained randomly placed elements and illustrates the impact of changes in reflection on a wave field. Figures allow visual identification of zones with wave amplification (where the scaled amplitude is higher than 1) and dissipation, in addition to phase shifts and the formation of partially standing waves.

The computational model is founded on well-established and well-utilized governing equations and series expansions using Hankel and Fourier series, as shown in Fig. 2. Furthermore, the wave field must by unique because it is decomposed into regular solutions that individually satisfy the Sommerfeld radiation condition (Courant and Hilbert 1962, p. 317). Thus, the differences between the AEM and the corresponding exact solution for each problem are isolated to how well the solutions for coefficients satisfy the coastal boundary conditions [Eq. (3)] along each of the vertical cylinders. This is quantified by computing the residual error in satisfying this condition across all control points for all elements, which is reported in Tables 1–3 as the RMS error normalized by dividing by the wave number and the amplitude of the incoming wave field. Thus, the error in satisfying the condition on $\partial \phi /\partial r$ was very small for all cases, and nearly exact solutions were achieved. The capacity of the numerical method to compute coefficients was also verified by comparing the coefficients obtained using the AEM method in Eq. (14) or Eq. (26) to the coefficient for the exact solution in Eq. (13) for the configurations in Fig. 1. The differences between the absolute value of difference in coefficients are $4\times {10}^{-15}\left|{\phi }_0\right|$ or smaller for all coefficients. Therefore, the AEM method is able to accurately reproduce the exact solutions and provide a very accurate set of benchmark solutions.

A visual interpretation of how well the coastal boundary conditions are satisfied can be found by examining the amplitude and phase of the wave fields. This is accomplished by separating Eq. (3) into real and imaginary parts to provide variation of these quantities in the normal direction, following (Berkhoff 1976):Thus, lines of constant wave amplitude are mutually orthogonal to the circular boundaries. Likewise, lines of constant phase intersect fully reflecting (R = 1) boundaries orthogonally and as a function of k for absorbing and partially reflecting boundaries. An inspection of each figure reveals that these contours have been shaped to satisfy these conditions for the specific configurations in the examples. For example, the strong wave interactions off a fully reflective coast in Fig. 5 are magnified in the upper panels of Fig. 10. This illustrates the wave interactions among groups of scattering elements that differ from the exact solutions of a single element in Fig. 1, and these images are further magnified in Fig. 10 to a quadrant of the circle. Clearly, the lines of constant amplitude and phase intersect this fully reflective boundary orthogonally to satisfy its boundary conditions [Eq. (15)].

A visual demonstration of the choice used for N terms in the series [Eq. (12a)], and the need for an overspecified number of control points M [Eq. (4)], is found in Fig. 11. The left panels illustrate the wave field obtained with N = 6 terms with M = (2 N + 1), as per Linton and Evans (1990, Eq. 2.15), for four closely spaced fully reflective cylinders with wavelengths corresponding to k r₀ = 5. This system of equations provided 13 complex boundary conditions with 13 complex coefficients for each element (26 equations for real and imaginary terms), and conditions were satisfied with a very small residual error at the computer accuracy of double-precision computations. However, boundary conditions were clearly not satisfied between the 13 control points, with lines of constant amplitude and phase intersecting the boundary nonorthogonally between these visualized points. The center panels used the same number of terms N but overspecified the number of control points by a factor of approximately 1.5. Clearly, the residual error was larger because conditions were satisfied by least squares; however, overall, a better match was obtained. The right panels show the results with the values of N and M used for all other figures and illustrate a better visual match of conditions and lower error. The model was also run for N = 50 and M = 152 with RMS error = 2.0 × 10⁻¹¹ $k\left|{\phi }_0\right|$, and the wave field (not shown) almost identically matched the curves shown for N = 20.

The iterative solution process requires numerically intensive evaluation of the Hankel functions and Fourier series in Eq. (12) that increases exponentially with the number of elements N in the application. Computational efficiencies may be enhanced by computing these functions once for the coefficients ${\mathcal{J}}_n$ and ${\mathcal{Y}}_n$ in Eq. (14) using Eq. (18). Furthermore, the Hankel and trigonometric functions in the additional functions may be precomputed before solving the problem, as shown in Appendix II, thereby reducing their evaluation to matrix multiplication. This organizes the solution methods to match those used for circular elements in groundwater flow problems by Barnes and Janković (1999), which have been applied to problems with 100,000 elements (Janković et al. 2009, Fig. 1). And yet, the coastal boundary conditions, Eqs. (3) and (15), may have difficulty achieving solutions for complicated wave interactions. Highly resonant problems with full reflection generate standing waves along features with tightly focused standing waves (e.g., on the element in the middle panel of Fig. 10), where small perturbations may cause significantly different wave fields to become manifest, and the iterative solution may jump across such solutions. The boundary condition for fully absorbing coasts requires waves to travel toward the coast with a normal derivative of the argument equal to minus the wave number [Eq. (15)], whereas the direction of incoming waves may be accounted for through addition terms, as in Steward and Panchang (2001). Thus, trapped or near-trapped modes may exist at irregular frequencies (Mei 1978) that cause the problem to become ill-posed. Such numerical instability is also observed at the corresponding fictitious frequencies in BEM formulations (Chen et al. 2009b, 2012).

The results provide a foundation toward future applications to engineering monitoring networks (Pandey et al. 2014) where sensors measure wave properties (amplitude and phase) at discrete locations to infer properties of wave interactions with partially reflecting coastal features. Clearly, a sensor network would have difficulty interpreting signals in collections of reflecting objects for shorter wavelengths (Figs. 3–5) because vastly different results might be measured over very small changes in the position of sensors. Thus, longer wavelengths were examined for regularly (Figs. 6–8) and randomly spaced objects (Fig. 9). Absorbing objects cause wave dissipation with a decrease in wave amplitude and a lag in phase that is shown to be impacted by the element size. Collections of reflecting objects may have zones of amplification that occur in both regularly and randomly placed objects and are also impacted by element size. These results suggest that design of monitoring networks aimed at distinguishing the reflection coefficient of collections of objects should utilize a sensor array strategically situated to measure both amplitude and phase across variations in the wavelength of incoming waves.

The AEM provides very accurate solutions to study waves in collections of partially reflecting coastal features. The published solutions using Hankel and trigonometric functions for a single circle provide a foundation for these developments. The AEM formulates the system of equations for boundary conditions in an iterative, element-by-element scheme that employs overspecification with least squares and uses orthogonality of Fourier series to achieve a solution via matrix multiplications. Solutions are developed for problems with many interacting elements, to enable computation of coefficients that satisfy a Robin boundary condition [Eq. (3)], along each element. Particular applications study the impact on wave fields resulting from variation in wavelength across elements, variation in the size of objects and the size of the gap between elements, and variation in the reflection coefficient for both regularly and randomly located elements. Water run-up is visualized for each figure by showing the contours of the amplification factors of the incoming plane waves.

This foundation enables identification of conditions under which resonance and dissipation of waves occur. Resonance exists in collections of coastal features with highly reflective boundaries, and dissipation occurs with more absorbing boundaries; additionally, both resonance and dissipation increase as objects become more tightly spaced. Phase shifts, where waves travel more slowly through the collection of objects, occur more readily for objects with absorbing boundaries. Errors in the computational implementation are quantified in Tables 1–3 and demonstrate the capacity of the AEM to provide a robust, nearly exact solution method for wave fields near circular objects.

Extension of the AEM to a broader set of coastal engineering problems than those embodied within the coastal boundary condition [Eq. (3)] could be envisioned, for example, by incorporate a phase shift occurring at the boundary through use of a complex reflection coefficient with real and imaginary parts (Berkhoff 1976). Such extensions could be used to study wave dissipation resulting from the uprush and downrush of waves in the swash zone of coastal structures (Kobayashi et al. 1987), the impact of nearshore vegetation (Kobayashi et al. 1993), and the direction at which waves intersect partially reflecting boundaries (Steward and Panchang 2001). Methods are extensible to other design problems with similar geometries, for example, vertical breakwaters surrounded by a rigid cylinder (Darwiche et al. 1994), vertical cylinders used with perforated breakwaters (Fugazza and Natale 1992; Kakuno and Liu 1993), and problems with porous coastal structures (Liu et al. 1999). These methods also provide a foundation for the study of circular shoals (Homma 1950; Liu et al. 2012), scour pits (Liu et al. 2013), and the formation of scour pockets around vertical cylinders (Sumer et al. 1992).

The details necessary to compute the coefficients of a coastal feature with matrix multiplication follow. Although the solution is presented for element i, a notation to simplify expressions is adopted whereby it is implied that coefficients and functions are evaluated for element i without explicitly specifying these indices. First, the complex wave function in Eq. (12a) with complex coefficients ($\stackrel{{}_{\text{cos}}}{c_n}=\Re \stackrel{{}_{\text{cos}}}{c_n}+\mathrm{i}\Im {\underset{\text{cos}}{c}}_n$ and ${\stackrel{\text{sin}}{c}}_n=\Re {\stackrel{\text{sin}}{c}}_n+\mathrm{i}\Im {\stackrel{\text{sin}}{c}}_n$) and Hankel function [$H_n^{\left(1\right)}\left(r\right)=J_n\left(r\right)+\mathrm{i}{Y}_n\left(r\right)$] is separated into real and imaginary parts:where the additional function $\stackrel{\text{add}}{\phi }$ in Eq. (12b) contains the contributions from all other elements. The boundary condition [Eq. (3)], is evaluated at control point m where r = r₀ and θ = θ_m and separated into real parts:and imaginary parts:using a condensed notation:

These partial derivatives of the additional function in the radial direction are obtained from the x- and y-gradient terms [Eq. (8)], with

The $m=1\cdots M$ conditions in Eqs. (17a and b) may be organized as two systems of equations in the matriceswhere A contains Fourier coefficientsand the unknown coefficients and additional functions are organized in the vectorsand

The least-squares solutiontakes on a simple form for series when the locations of control points are specified at constant intervals using θ_m in Eq. (4). In this case, the orthogonality of the Fourier series results in a diagonal matrix D as in (Steward 2015):with terms equal to M for the first diagonal entry and M/2 for the others. This gives a solution where the unknown coefficients may be computed directly fromwhich may be written in summation form asand

The coefficients are rearranged by solving sets of two equations with two unknowns for each n to give the real and imaginary components:and

These are added to give the complex coefficients:that are used to obtain the solutions with coefficients from Eq. (18).

Computation of coefficients for element i in Eq. (14) requires evaluation of the additional function [Eq. (12b)] at each point $m=1,\dots ,M$ along its boundary:

The computations associated with the additional functions may be organized as matrix multiplications within the matrices:

This provides a direct relation between the coefficients of each analytic element and their influence on the boundary condition of the other elements:

This is numerically advantageous because the matrices $\underset{\mathit{ij}}{\overset{\hspace{0.17em}\text{cos}\hspace{0.17em}}{\mathbf{F}}}$ and $\underset{\mathit{ij}}{\overset{\hspace{0.17em}\text{sin}\hspace{0.17em}}{\mathbf{F}}}$ are computed once, and changes in φ at element i as a result of adjusting coefficients in element j are computed with matrix multiplication rather than evaluating Bessel and Fourier functions. Similar matrices follow directly to provide the derivative of the additional function in Eq. (14) using matrix multiplication.

This study of reflection and adsorption contributes toward research of electromagnetic waves near wheat plants for the USDA/NIFA Award 2017-67007-25943 and interpretation of electrical resistivity in wetted soils for the USDA/ARS Ogallala Aquifer Program.