View source for Boundary Element Method for a Fixed Body in Finite Depth

{{complete pages}}

== Introduction ==

We show how to solve for the wave scattering by a rigid body in constant [[Finite Depth]] using
the [[:Category:Boundary Element Method|Boundary Element Method]].

== Equations ==

{{boundary value problem for a fixed body}}

{{incident plane wave}}

{{sommerfeld radiation condition two dimensions}}

== Solution Method ==

We divide the domain into three regions, <math>x<-l</math>,
<math>x>r</math>, and <math>-l<x<r</math> so that the body surface is entirely in the
finite region. 

=== Solution in the finite region ===

We use the [[:Category:Boundary Element Method|Boundary Element Method]]
in the finite region.

===Solution in the Semi-infinite Domains===

We now solve Laplace's equation in the semi-infinite domains. First consider
the domain on the left so that <math>\Omega = \left\{ x<-l,\;-h\leq z\leq 0\right\} </math> 
The equations are
{{standard linear wave scattering equations without body condition}}
We have the following explicit boundary conditions
<center><math>
\phi =f\left( z\right) ,\;\;x=-\,l
</math></center>
and
<center><math>
\lim\limits_{x\rightarrow -\infty }\phi \left( x,z\right) = \phi_0
e^{-k_{0}x}
+R\phi_0
e^{k_{0}x},
</math></center>
where <math>f\left( z\right) </math> is
an arbitrary continuous function<math>.</math> Our aim is to find the outward normal
derivative of the potential on <math>x=-l</math> as a function of <math>\tilde{\phi}\left(
z\right) </math>. <math>\phi_0</math> and <math>k_0</math> will be defined shortly
when we separate variables, but are equivalent to the [[Sommerfeld Radiation Condition]].

Since the water depth is constant in these regions we can solve Laplace's
equation by separation of variables. 

{{separation of variables for a free surface}}

=== Expansion of the potential ===

This gives us the following expression for the potential in the region <math>
\Omega </math> 
<center><math>
\phi \left( x,z\right)  = \phi_0(z)e^{-k_0 x}
+\sum_{m=0}^{\infty } a_m  \phi _{m}\left( z\right)
e^{k_{m}\left( x+l\right) }
</math></center>
where <math>a_m</math> can be found from <math>\left. \phi \right| _{x=-l} = f(z)</math>
from the orthogonality of the vertical eigenfunctions. Therefore
<center><math>
a_m +\delta_{m0}e^{k_0 l} = \frac{1}{A_m} \int_{-h}^{0} f(z) \phi_m(z) \mathrm{d}z
</math></center>

The normal derivative of potential (with the normal outward from the interior region) we obtain
<center><math>
\partial_n\left. \phi \right| _{x=-l} 
= -\partial_x \left. \phi \right| _{x=-l} 
= k_0 \phi_0(z)e^{k_0 l}
- \sum_{m=0}^{\infty } k_m a_m  \phi _{m}\left( z\right)
</math></center>

If we define <math>\mathbf{Q}</math> as 
<center><math>
\mathbf{Q} \left[f\left( z\right) \right]=
-\sum_{m=0}^{\infty } \frac{k_{m}}{A_m} \left\{ \int_{-h}^{0} f\left( s\right) \phi
_{m}\left( s\right) \mathrm{d}s \right\} \phi _{m}\left( z\right) 
</math></center>
we can write
<center><math>
\partial_n \left. \phi \right| _{x=-l} =\mathbf{Q} \left[ f(z) \right] + 2k_0 \phi_0e^{k_0l},
</math></center>

Similarly, we now consider the potential in the region 
<math>\Omega = \left\{ x>r,\;-h\leq z\leq 0\right\} </math> 
which satisfies exactly the equations as before except the boundary
condition is  
<center><math>
\lim\limits_{x\rightarrow \infty }\phi \left( x,z\right) = \phi_0
T e^{-k_{0}x}
</math></center>
We solve again by separation of variables and obtain
<center><math>
\left. \partial_n \phi \right| _{x=r}=
\left. \partial_x \phi \right| _{x=r}=\mathbf{Q}\left[f\left( z\right)\right] ,
</math></center>
where the outward normal is with respect to the inner domain as before. Note that, if
the depth is different on each side then the matrix <math>\mathbf{Q}</math> will be different
(since we need to use a different depth solving the [[Dispersion Relation for a Free Surface]] etc.)

== Boundary Element Method ==

We have reduced the problem to Laplace's equation in a finite domain subject
to certain boundary conditions. These boundary
conditions give the outward normal derivative of the potential as a function
of the potential but this is not always a point-wise condition; on some
boundaries it is given by an integral equation. We must solve both Laplace's
equation and the integral equations numerically. We will solve Laplace's
equation by the [[:Category:Boundary Element Method|Boundary Element Method]] and the integral equations by
numerical integration. However, the same discretisation of the boundary will
be used for both numerical solutions.
We solve Laplace's equation by a modified constant panel method
which reduces it to the following matrix equation 
<center><math>
\frac{1}{2}\vec{\phi}=\partial_n\mathbf{G}\vec{\phi}-\mathbf{G}\partial_n\vec{\phi}.
</math></center>
where <math>\vec{\phi}\mathcal{\ }</math>and <math>\partial_n\vec{\phi}</math> 
are vectors which approximate the potential and its normal derivative
around the boundary <math>\partial \Omega </math>, and <math>\mathbf{G}</math> and <math>\mathbf{G}_{n}</math>
are matrices corresponding to the Green function and the outward normal
derivative of the Green function respectively. The method used to calculate
the elements of the matrices <math>\mathbf{G}</math> and <math>\partial_n\mathbf{G}</math> can be
found in [[:Category:Boundary Element Method|Boundary Element Method]].

The outward normal derivative of the potential, <math>\partial_n\vec{\phi},</math> and the
potential, <math>\vec{\phi},</math> are related by the conditions on the boundary <math>
\partial \Omega </math>. This can be expressed as 
<center><math>
\partial_n\vec{\phi}=\mathbf{A}\,\vec{\phi}-\vec{f}, 
</math></center>
where <math>\mathbf{A}</math> is the block diagonal matrix given by 
<center><math>
\mathbf{A}\mathbb{=}\left[ 
\begin{matrix}
\mathbf{Q} &  &  &   \\ 
& \mathbf{Q} &   &  \\ 
&    & \alpha \,\mathbf{I} &  \\ 
&    &  & 0
\end{matrix}
\right] ,  
</math></center>
where we have separated the boundaries to the semi-infinite domains, the free surface and the 
sea floor and body boundary. 
<math>\mathbf{Q}</math> is a matrix
approximation of the integral operator of the same name and <math>\vec{f}</math> is
the vector 
<center><math>
\vec{f}=\left[ 
\begin{matrix}
-2k_{0}\phi_0 \,e^{k_{0}l } \\ 
0 \\ 
\vdots \\ 
0
\end{matrix}
\right] . 
</math></center>

We obtain the following matrix equation for the
potential 
<center><math>
\left( \frac{1}{2}\ \mathbf{I}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G}
\vec{f}
</math></center>
which can be solved straightforwardly. The reflection and transmission
coefficients are calculated from the solution.

{{numerical calculation of Q}}

{{reflection and transmission coefficients BEM}}

=== Singularity Expansion Method ===
The matrix equation for the potential 
<center><math>
\left( \frac{1}{2}\ \mathbf{I}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G}
\vec{f}
</math></center>
can be rewritten as 
<center><math>
 \mathbf{M}(k) \mathbf{x}=\mathbf{F}(k)
</math></center>
where <math>\mathbf{x}</math> represents the unknown potential values in the new notation. Resonant frequencies (or scattering frequencies) correspond to zero eigenvalues of the matrix <math>\mathbf{M}(k)</math> for a set of discrete points <math>k=k_p</math> in the upper half of the complex plane (due to the time dependence). Since 
<center><math>
\mathbf{x}=\left( \mathbf{M}(k)\right)^{-1}\mathbf{F}(k)
</math></center>
then <math>\mathbf{x}</math> will be singular for these scattering frequencies. Each scattering frequency, corresponding to a zero eigenvalue of <math>\mathbf{M}(k)</math> will possess a corresponding eigenvector <math>\mathbf{u}_{k_p}</math> satisfying
<center><math>
\mathbf{M}(k_p)\mathbf{u}_{k_p}=0
</math></center>
and also a related conjugate transpose eigenvector <math>\mathbf{u}^{*}_{k_p}</math> satisfying 
<center><math>
\mathbf{M}(k_p)^{A}\mathbf{u}^*_{k_p}=0.
</math></center>
In the vicinity of a scattering frequency it can be shown that the vector of unknowns satisfies
<center><math>
\mathbf{x}(k)\approx \frac{\langle \mathbf{F}(k_p),\mathbf{u}^*_{k^*_p} \rangle}{\langle\mathbf{u}_{k_p},\mathbf{M}'(k_p)\mathbf{u}^*_{k^*_p}\rangle}\frac{1}{k-k_p}
</math></center>
and we define a corresponding mode shape as
<center><math>
\zeta_{k_p}(x,k_p) =\frac{\langle \mathbf{F}(k_p),\mathbf{u}^*_{k^*_p} \rangle}{\langle\mathbf{u}_{k_p},\mathbf{M}'(k_p)\mathbf{u}^*_{k^*_p}\rangle}U(x)
</math></center>
where <math>U(x)</math> is obtained

== Matlab Code ==

{{fixed body 2d bem code}}

=== Additional code ===

This program requires
* {{free surface dispersion equation code}}
* {{boundary element code}}

[[Category:Boundary Element Method]]
[[Category:Eigenfunction Matching Method]]