SMS:CGWAVE Math Details: Boundary Conditions

Domains on which the elliptic eq. (5) is solved are enclosed by closed boundaries (represented by coastlines and surface-penetrating structures like pier walls or pier legs, breakwaters, seawalls, etc.) and open boundaries (which represent an artificial boundary between the area being modelled and the sea region outside. A separation between the model domain and an outer water area from where no waves enter the model domain (e.g. a creek or tributary at the backbay or down wave end of the domain) may be considered to be a fully-absorbing closed boundary. An open boundary is considered to be one where an incident wave is specified (and may contain other radiated waves). Along these boundaries, appropriate conditions must be specified to solve (1); however, even in the best of circumstances, only approximate boundary conditions can be developed (e.g. see Dingemaans, 1997).

Closed Boundary Conditions

Along coastline and surface-protruding structures, the following boundary condition has traditionally been used (e.g. Berkhoff 1976; Tsay & Liu 1983; Tsay et al. 1989; Oliveira and Anastasiou 1998; Li 1994a)

${\frac {\partial \phi }{\partial h}}=\alpha \phi$ ⁽⁶⁾

where n is the outward normal to the boundary and a is related to a user-specified reflection coefficient as follows:

$\alpha =ik{\frac {1-K_{r}}{1+K_{r}}}$ ⁽⁷⁾

K_r varies between 0 and 1 and specific values for different types of reflecting surfaces have been compiled by Thompson et al. (1996).

It may be verified that (6) is strictly valid only for fully-reflecting boundaries (i.e. K_r= 1). For partially reflecting boundaries, it is valid only if waves approach the boundary normally. For other conditions, (6) is approximate, and may produce distortions in the model solutions. These limitations may be eliminated by describing the solution at the boundary more fully as a sum of incident and reflected waves:

$\phi =A\{\exp[ik(n\,\cos \theta +s\,\sin \theta )]+K,\ \exp[ik(-n\,\cos \theta +s\,\sin \theta +\beta )]\}$ ⁽⁸⁾

where A is the amplitude of the approaching waves, θ is the direction at which they intersect the boundary (θ = 0 for normally incident waves), s is the coordinate along the tangent to the boundary, and β is a phase shift between the incident and the reflected wave. (8) leads to the following boundary condition.)

${\frac {\partial \phi }{\partial n}}=ik\cos \theta {\frac {1-K_{r}\exp(ik\beta )}{1+K_{r}\exp(ik\beta )}}\phi$ ⁽⁹⁾

Unfortunately, θ and β are not known a priori inside the model domain, and must be estimated by approximation. For fully absorbing boundaries K_r = 0), Li and Anastasiou (1992) and Li et al. (1993) have used (9) after estimating θ from Snell's Law and the deep-water incident wave angle. Alternatively, Isaacson and Qu (1990) estimated θ as follows:

$\theta =\arctan\{(\partial \chi /\partial s)/(\partial \chi /\partial n)\}$ ⁽¹⁰⁾

where χ is the argument of the complex quantity θ (i.e. the phase of θ). For implementation, they first used (6) as a boundary condition, obtained χ from the results, determined θ from (10), used (9) as a boundary condition to perform a second iteration of the model, recalculated χ and θ, performed a third model iteration using (9), and so on. Like Pos (1985), they assumed β = 0 while using (9), based on limited numerical tests that showed little sensitivity to β. Clearly, like the Snell’s Law approach, (10) is valid only for K_r = 0 (although problems with non-zero K_r were also considered). To include the effect of the reflected waves (i.e the second term in the right hand side of (8)), Isaacson et al. (1993) suggested estimating θ as follows:

$\theta =(1/k)\arcsin\{\partial \chi /\partial s\}$ ⁽¹¹⁾

Again an iterative method with repeated model calculations were needed. Steward and Panchang (2000) analyzed these methods and noted difficulties with convergence of the above iterative methods and with the quality of the solutions obtained with (10) and (11). They were able to eliminate these difficulties by estimating θ from the following expression:

$\tan \theta ={\frac {\partial \chi }{\partial s}}{\Bigg /}\left({\frac {\partial \chi }{\partial n}}+{\frac {2K_{r}k(\cos(k\beta )+K_{r}}{1+2K_{r}\cos(k\beta )+K_{r}^{2}}}\cos(\theta )\right)$ ⁽¹²⁾

(12) is a generalization of (10) that allows non-zero K_r and β. For a detailed comparison of results, see Steward and Panchang (2000). Fig. 2 shows a simulation with of waves propagating into a rectangular harbor area. Clearly, solutions obtained with (12) are qualitatively superior to those obtained with (10).

Despite the increasing sophistication seen progressively in (6) and (9) and in the various ways of estimating θ, some fundamental problems remain. The most important one is inherent in (8), i.e. the assumption that the total wave field near the boundary can be represented either by one set of plane waves (in the case of (10) or the Snell's Law approach of Li et al. (1993)) or by two sets of plane waves (in the case of (11) and (12)) propagating in constant depth. In domains of complex shapes (as in Fig. 1) with arbitrary bathymetry and boundaries with varying reflectivities, a complex pattern of waves can result; simple wave trains are not easily discernible and, as noted by Isaacson and Qu (1990), the definition of a single θ (and β) can become meaningless. Further, even when there is a well-defined train of waves near the boundary (justifying the use of the above methods), precise estimation of K_r and β is still problematic. Values of K_r provided by Thompson et al. (1996) certainly do not cover full range of reflecting surfaces that the modeller encounters, nor do they cover the dependence of these parameters on the incident wave frequency. Efforts to incorporate the work of Dickson et al. (1995) and Sutherland and O’Donoghue (1998) pertaining to β in models such as the one described here are lacking.

In some ways it may be best to recognize these difficulties at the outset and use the simplest expressions (6) and (7) by combining all the uncertainties noted above into a single parameter α, which may be regarded as a tuning parameter. This is the approach followed in CGWAVE.

Open Boundary Conditions

Along the open boundary, an incident wave φ_i must be specified. Along this boundary, however, waves backscattered from within the domain will also exist, and their magnitude is generally not known. In the context of simple rectangular domain models, with one side (aligned, say, in the y direction) constituting the open boundary, Panchang et al. (1988, 1991), Li (1994a, b), and Oliveira and Anastasiou (1998) have used the following condition

${\frac {\partial \phi }{\partial x}}=ik(2\phi _{i}-\phi )$ ⁽¹³⁾

(13) is obtained by assuming that the incident and backscattered components along this boundary can be described by $\phi _{i}=A_{i}\exp(ikx)$ and $\phi =B\exp(-ikx)$ respectively (where A_i is the (specified) amplitude of the incoming wave and B is an unknown), adding the two components, and differentiating. Obviously, this is valid only if the incident and backscattered waves near the boundary are plane waves propagating in the +/- x direction.

For more complex domains involving multidirectional scattering, (13) is inappropriate. Harbor applications generally use model domains such as that described in Fig. 3, where the semicircle is used to separate the model area from the open sea. In the exterior domain Ω' the potential φ is comprised of three components:

$\phi =\phi _{i}+\phi _{r}+\phi _{s}^{}$ ⁽¹⁴⁾

where φ_i = the incident wave that must be specified to force the model, φ_r = a reflected wave that would exist in the absence of the harbor, and φ_s = a scattered wave that emanates as a consequence of the harbor and must satisfy the Sommerfeld radiation condition. With appropriate descriptions for these components, a boundary condition can be developed along the semicircle.

In traditional harbor models (Mei 1983; Tsay and Liu, 1983; Thompson et al. 1996; Chen and Houston, 1987; Xu and Panchang, 1993; Demirbilek and Panchang, 1998), the exterior wave conditions are described as follows:

$\phi _{i}=A_{i}[ikr\cos(\theta -\theta _{i}^{})]$ ⁽¹⁵⁾, which is the specified input

$\phi _{r}=A_{i}[ikr\cos(\theta -\theta _{i}^{})]$ ⁽¹⁶⁾

$\phi _{s}=\sum _{n=0}^{\phi }H_{n}(kr)(A_{n}\cos n\theta +B_{n}\sin n\theta )$ ⁽¹⁷⁾

where (r, φ) denotes the location of a point in polar coordinates, H_n is the Hankel function of the first kind and order n, and A_n and B_n are unknown coefficients.

For the specified incident wave field given by (15), equations (16) and (17) result from the solution of the relevant eigenvalue problem in the traditional method. As demonstrated by Xu et al. (1996), however, this eigenvalue problem, in which φ_s and φ_r are coupled, may be solved only under the following conditions:

(i) The exterior region must have a constant depth,

(ii) the exterior coastlines A1D1 and A2D2 must be fully reflecting and collinear.

These requirements usually cannot be met in practice where the exterior geometry varies arbitrarily, and the unrealistic bathymetric representation used perforce by the modeller invariably has an adverse influence on the solution. In field applications, the exterior bathymetry is irregular and the depth generally increases in the x-direction. Condition (i) is thus violated, causing two problems as demonstrated by Panchang et al. (2000). First, the modeller must arbitrarily select a representative “constant” depth and test the sensitivity of the solutions to these depths. This can be extremely time-consuming. Second, the effect of reflections from the sloping exterior bathymetry is ignored. These effects are often significant, especially for long periods that are of interest in harbor resonance studies. Condition (ii) is also problematic. Exterior coastlines are not always fully reflecting for all wave conditions, and imposing full reflection in such cases yields extremely large amplification factors and rapid variations in the wave pattern in the outer regions of the domain. (See examples in Xu et al. (1996), Demirbilek et al. (1996), and in Thompson et al. 1996)). One may of course enlarge the interior region in the hope that these effects do not contaminate the results in the area of interest; however, there is no guarantee that these effects are confined to specific regions. In addition, the extra memory requirements and grid-generation for a larger domain are usually exceedingly demanding. Thus, while (16)-(17) constitute rigorous solutions of the eigenvalue problem, their use renders the application of harbor wave models problematic in practice. (One consequence of the above is that many of the models in this category cannot correctly simulate fairly simple phenomena like waves approaching a sloping beach. Investing confidence in model results when applied to field situations is therefore difficult.)

An effective alternative, followed in CGWAVE, is to use a “parabolic approximation” to describe $\phi _{s}$:

${\frac {\partial \theta _{s}}{\partial _{r}}}=-p\phi _{s}-q{\frac {\partial ^{2}\phi _{s}}{\partial \theta ^{2}}}$ ⁽¹⁸⁾

where

$p=ik_{0}+{\frac {1}{2r}}-{\frac {i}{8k_{0}r^{2}}}$ , ${\frac {d}{dx}}(CC_{g}{\frac {d\psi }{dx}})+kCC_{g}(k\cos ^{2}\theta +iW)\psi =0$ ^(18a)

where r and θ represent the polar coordinates of a point on the open boundary. (p and q are not unique, and alternative forms, each obtained with an appropriate rationale, have been investigated by Givoli (1991), Xu et al. (1996), Panchang et al. (2000)). The parabolic approximation (18) allows the scattered waves to exit only through a limited aperture around the radial direction. Unlike (17), therefore, it does not rigorously satisfy the Sommerfeld radiation condition. However, using this formulation decouples φ_s from the other components. These components (φ_i and Φ_r) may be obtained by making a compromise between a detailed exterior bathymetric representation (which as noted earlier, is difficult) and the constant depth representation (which is unrealistic). A one-dimensional representation, where the depths vary in the cross-shore direction only (Fig. 1), may be selected. This is reasonable, since in general, this is often the direction in which the depths vary the most. If natural variations do not permit the representation of the exterior depths by only one section, a second one-dimensional section, shown as transect 2 in Fig. 1, may be constructed. For transects 1 and 2 with varying depths, no simple analytical expression (such as (16)) can be found for the reflected wave (since &phi_i and &phi_r are coupled). However, the quantity:

$\phi _{0}=\phi _{i}+\phi _{r}^{}$ ⁽¹⁹⁾

may be obtained by the solution of the one-dimensional version of (5), since the depths along these transects vary in one direction only. This one-dimensional equation is (Schaffer and Jonsson, 1992; Panchang et al. 2000):

${\frac {d}{dx}}(CC_{g}{\frac {d\psi }{dx}})+kCC_{g}(k\cos ^{2}\theta +iW)\psi =0$ ⁽²⁰⁾

where, for one-dimensional geometry,

$\phi _{0}^{}=\psi (x)\exp(iky\sin \theta )$ ⁽²¹⁾

(20) is an elliptic ordinary differential equation requiring two boundary conditions. It may easily be solved via a simple finite-difference scheme. (For the present, the dissipation factor W is considered to be prespecified). Assuming that transect 1 extends out to a region of constant depth (or deep water), a condition at P₁ may be obtained by combining a specified incident wave

$\phi _{i}^{}(P_{1})=A_{i}\exp(ikx\cos \theta _{i}+iky\sin \theta _{i})$ ⁽²²⁾

(where $A_{i}$ = a given input wave amplitude) and an unknown reflected wave:

$\phi _{r}^{}(P_{1})=B\exp(-ikx\cos \theta _{i}+iky\sin \theta _{i})$ ⁽²³⁾

Without loss of generality, the point P₁ may be located at x = 0, which allows elimination of B to yield

${\frac {\partial \psi }{\partial x}}=ik\cos \theta _{i}(2A_{i}-\psi )$ ⁽²⁴⁾

At the coastal boundary point Q₁, the partial reflection boundary condition (9) may be used in the following form:

${\frac {\partial \psi }{\partial x}}={\frac {i{\sqrt {k^{2}-k^{2}\sin ^{2}\theta }}(1-K_{r})}{1+Kr}}\psi$ ⁽²⁵⁾

where K_r is the reflection coefficient for the exterior coastline (i.e. near Q₁) and $k\sin \theta$ is constant for the one-dimensional problem.

The solution of (20) using boundary conditions (24) and (25) along with (21) produces φ₀ along transects 1 and 2. These solutions are denoted by φ₀₁ and φ_02}. The desired φ₀ along the semicircle may be obtained by laterally translating φ₀₁ and φ₀₂ via interpolation between transects 1 and 2 as follows:

$\phi _{0}=(1-m)\phi _{01}^{}\exp(-ik(r-y)\sin \theta )+m\phi _{02}\exp(ik(r+y)\sin \theta )$ ⁽²⁶⁾

where we have set y = 0 at the center of semicircle, the interpolation function $m=(r-y)/2r$, r is the radius of the semicircle, y is the lateral coordinate of the open boundary node relative to the origin of semicircle (Fig.3).

The boundary condition for φ along the semicircle Γ may be obtained by using the continuity of the potential (equations 14 and 19) and its derivative along with (18) and (26):

${\frac {\partial \psi }{\partial x}}={\frac {i{\sqrt {k^{2}-k^{2}\sin ^{2}\theta }}(1-K_{r})}{1+K_{r}}}\psi$ ⁽²⁷⁾

Thus, the solution of (20) provides φ₀ along the 1-dimensional transects. These values can be translated laterally and substituted into (27) to obtain the open boundary condition for the two-dimensional equation (1). Zhao et al. (2000) and Panchang et al. (2000) have demonstrated that this procedure provides extremely satisfactory solutions for a large number of test cases.

In CGWAVE, if non-linear breaking is implemented, a version of (21) is solved by iteration since W depends on the (unknown) amplitude.

Related Topics

• CGWAVE
• Governing Equations
• Numerical Solution