Difference between revisions of "Boundary Element Method for a Fixed Body in Finite Depth"

From WikiWaves
Jump to navigationJump to search
Line 158: Line 158:
 
potential  
 
potential  
 
<center><math>
 
<center><math>
\left( \frac{1}{2}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G}
+
\left( \frac{1}{2}\ \mathbf{I}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G}
 
\vec{f}
 
\vec{f}
 
</math></center>
 
</math></center>
Line 169: Line 169:
 
evanescent modes,
 
evanescent modes,
 
<center><math>
 
<center><math>
\mathbf{Q}\phi =\sum_{m=0}^{N}k_{m}\int_{-h}^{0} \phi
+
\mathbf{Q}\phi =\sum_{m=0}^{N}k_{m}\int_{-h}^{0} \tilde{\phi}
\left( s\right) ,\phi_{m}\left( s\right) \mathrm{d}s
+
\left( s\right) \phi_{m}\left( s\right) \mathrm{d}s
 
\frac{\phi _{m} \left( z\right)}{A_m} .
 
\frac{\phi _{m} \left( z\right)}{A_m} .
 
</math></center>
 
</math></center>

Revision as of 23:07, 16 June 2009

Introduction

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

Equations

The Standard Linear Wave Scattering Problem in Finite Depth for a fixed body is

[math]\displaystyle{ \begin{align} \Delta\phi &=0, &-h\lt z\lt 0,\,\,\mathbf{x} \in \Omega \\ \partial_z\phi &= 0, &z=-h, \\ \partial_z \phi &= \alpha \phi, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]


(note that the last expression can be obtained from combining the expressions:

[math]\displaystyle{ \begin{align} \partial_z \phi &= -\mathrm{i} \omega \zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \\ \mathrm{i} \omega \phi &= g\zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]

where [math]\displaystyle{ \alpha = \omega^2/g \, }[/math]) The body boundary condition for a rigid body is just

[math]\displaystyle{ \partial_{n}\phi=0,\ \ \mathbf{x}\in\partial\Omega_{\mathrm{B}}, }[/math]

The equation is subject to some radiation conditions at infinity. We assume the following. [math]\displaystyle{ \phi^{\mathrm{I}}\, }[/math] is a plane wave travelling in the [math]\displaystyle{ x }[/math] direction,

[math]\displaystyle{ \phi^{\mathrm{I}}(x,z)=A \phi_0(z) e^{\mathrm{i} k x} \, }[/math]

where [math]\displaystyle{ A }[/math] is the wave amplitude (in potential) [math]\displaystyle{ \mathrm{i} k }[/math] is the positive imaginary solution of the Dispersion Relation for a Free Surface (note we are assuming that the time dependence is of the form [math]\displaystyle{ \exp(-\mathrm{i}\omega t) }[/math]) and

[math]\displaystyle{ \phi_0(z) =\frac{\cosh k(z+h)}{\cosh k h} }[/math]

In two-dimensions the Sommerfeld Radiation Condition is

[math]\displaystyle{ \left( \frac{\partial}{\partial|x|} - \mathrm{i} k \right) (\phi-\phi^{\mathrm{{I}}})=0,\;\mathrm{{as\;}}|x|\rightarrow\infty\mathrm{.} }[/math]

where [math]\displaystyle{ \phi^{\mathrm{{I}}} }[/math] is the incident potential.

Solution Method

We divide the domain into three regions, [math]\displaystyle{ x\lt -l }[/math], [math]\displaystyle{ x\gt r }[/math], and [math]\displaystyle{ -l\lt x\lt r }[/math] so that the body surface is entirely in the finite region.

Solution in the finite region

We use the Boundary Element Method in the finite region.

Solution in the Semi-infinite Domains

We now solve Laplace's equation in the semi-infinite domains [math]\displaystyle{ \Omega ^{-}=\left\{ x\lt -l,\;-h\leq z\leq 0\right\} }[/math] and [math]\displaystyle{ \Omega ^{+}=\left\{ x\gt r,\;-h\leq z\leq 0\right\} }[/math]. Since the water depth is constant in these regions we can solve Laplace's equation by separation of variables. The potential in the region [math]\displaystyle{ \Omega ^{-} }[/math] satisfies the following equation

[math]\displaystyle{ \left. \begin{matrix} \nabla ^{2}\phi =0,\;\;\mathbf{x}\in \Omega ^{-}, \\ \phi _{n} = \alpha \phi =0,\;\;z=0, \\ \phi _{n}=0,\;\;z=-h, \\ \phi =\tilde{\phi}\left( z\right) ,\;\;x=-\,l, \\ \lim\limits_{x\rightarrow -\infty }\phi \left( x,z\right) = \phi_0 e^{-k_{0}x} +R\phi_0 e^{k_{0}x}, \end{matrix} \right\} }[/math]

where [math]\displaystyle{ \mathbf{x}=\left( x,z\right) \ }[/math]and [math]\displaystyle{ \tilde{\phi}\left( z\right) }[/math] is an arbitrary continuous function[math]\displaystyle{ . }[/math] Our aim is to find the outward normal derivative of the potential on [math]\displaystyle{ x=-l }[/math] as a function of [math]\displaystyle{ \tilde{\phi}\left( z\right) }[/math]. [math]\displaystyle{ \phi_0 }[/math] and [math]\displaystyle{ k_0 }[/math] will be defined shortly when we separate variables, but are equivalent to the Sommerfeld Radiation Condition.

Separation of variables for a free surface

We use separation of variables

We express the potential as

[math]\displaystyle{ \phi(x,z) = X(x)Z(z)\, }[/math]

and then Laplace's equation becomes

[math]\displaystyle{ \frac{X^{\prime\prime}}{X} = - \frac{Z^{\prime\prime}}{Z} = k^2 }[/math]

The separation of variables equation for deriving free surface eigenfunctions is as follows:

[math]\displaystyle{ Z^{\prime\prime} + k^2 Z =0. }[/math]

subject to the boundary conditions

[math]\displaystyle{ Z^{\prime}(-h) = 0 }[/math]

and

[math]\displaystyle{ Z^{\prime}(0) = \alpha Z(0) }[/math]

We can then use the boundary condition at [math]\displaystyle{ z=-h \, }[/math] to write

[math]\displaystyle{ Z = \frac{\cos k(z+h)}{\cos kh} }[/math]

where we have chosen the value of the coefficent so we have unit value at [math]\displaystyle{ z=0 }[/math]. The boundary condition at the free surface ([math]\displaystyle{ z=0 \, }[/math]) gives rise to:

[math]\displaystyle{ k\tan\left( kh\right) =-\alpha \, }[/math]

which is the Dispersion Relation for a Free Surface

The above equation is a transcendental equation. If we solve for all roots in the complex plane we find that the first root is a pair of imaginary roots. We denote the imaginary solutions of this equation by [math]\displaystyle{ k_{0}=\pm ik \, }[/math] and the positive real solutions by [math]\displaystyle{ k_{m} \, }[/math], [math]\displaystyle{ m\geq1 }[/math]. The [math]\displaystyle{ k \, }[/math] of the imaginary solution is the wavenumber. We put the imaginary roots back into the equation above and use the hyperbolic relations

[math]\displaystyle{ \cos ix = \cosh x, \quad \sin ix = i\sinh x, }[/math]

to arrive at the dispersion relation

[math]\displaystyle{ \alpha = k\tanh kh. }[/math]

We note that for a specified frequency [math]\displaystyle{ \omega \, }[/math] the equation determines the wavenumber [math]\displaystyle{ k \, }[/math].

Finally we define the function [math]\displaystyle{ Z(z) \, }[/math] as

[math]\displaystyle{ \chi_{m}\left( z\right) =\frac{\cos k_{m}(z+h)}{\cos k_{m}h},\quad m\geq0 }[/math]

as the vertical eigenfunction of the potential in the open water region. From Sturm-Liouville theory the vertical eigenfunctions are orthogonal. They can be normalised to be orthonormal, but this has no advantages for a numerical implementation. It can be shown that

[math]\displaystyle{ \int\nolimits_{-h}^{0}\chi_{m}(z)\chi_{n}(z) \mathrm{d} z=A_{n}\delta_{mn} }[/math]

where

[math]\displaystyle{ A_{n}=\frac{1}{2}\left( \frac{\cos k_{n}h\sin k_{n}h+k_{n}h}{k_{n}\cos ^{2}k_{n}h}\right). }[/math]

Expansion of the potential

This gives us the following expression for the potential in the region [math]\displaystyle{ \Omega ^{-}, }[/math]

[math]\displaystyle{ \phi \left( x,z\right) = \phi_0(z)e^{k_0 x} +\sum_{m=0}^{\infty }\left\langle \tilde{\phi}\left( z\right) ,\phi _{m}\left( z\right) \right\rangle \frac{\phi _{m}\left( z\right)}{A_m} e^{k_{m}^{\left( 1\right) }\left( x+l\right) } }[/math]

We define [math]\displaystyle{ \mathbf{Q} }[/math] as

[math]\displaystyle{ \mathbf{Q}\tilde{\phi}\left( z\right) = \sum_{m=0}^{\infty } k_m \int_{-h}^{0} \tilde{\phi}\left( s\right) \phi _{m}\left( s\right) ds \frac{\phi _{m}\left( z\right)}{A_m} e^{k_{m}^{\left( 1\right) }\left( x+l\right) } }[/math]

The normal derivative of the potential on the boundary of [math]\displaystyle{ \Omega ^{-} }[/math] and [math]\displaystyle{ \Omega }[/math] [math]\displaystyle{ \left( x=-l\right) }[/math] gives us

[math]\displaystyle{ \left. \phi _{n}\right| _{x=-l}=\mathbf{Q}\tilde{\phi}\left( z\right) -2k_0 \phi_0e^{-k_0l}, }[/math]

Similarly, we now consider the potential in the region [math]\displaystyle{ \Omega ^{+} }[/math] which satisfies

[math]\displaystyle{ \left. \begin{matrix} \nabla ^{2}\phi =0,\;\;\mathbf{x}\in \Omega ^{+}, \\ \phi _{n}=\alpha \phi,\;\;z=0, \\ \phi _{n}=0, \;\;z=-h, \\ \phi =\tilde{\phi}\left( z\right) ,\;\;x=\,r, \\ \lim\limits_{x\rightarrow -\infty }\phi \left( x,z\right) =T\phi_0(z) e^{k_{0}x}. \end{matrix} \right\} }[/math]

Solving by separation of variables as before we obtain

[math]\displaystyle{ \left. \phi _{n}\right| _{x=r}=\mathbf{Q}\tilde{\phi}\left( z\right) , }[/math]

where the outward normal is with respect to the [math]\displaystyle{ \Omega }[/math] domain.

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 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

[math]\displaystyle{ \frac{1}{2}\vec{\phi}=\mathbf{G}_{n}\vec{\phi}-\mathbf{G}\vec{\phi}_{n}. }[/math]

where [math]\displaystyle{ \vec{\phi}\mathcal{\ } }[/math]and [math]\displaystyle{ \vec{\phi} _{n} }[/math] are vectors which approximate the potential and its normal derivative around the boundary [math]\displaystyle{ \partial \Omega }[/math], and [math]\displaystyle{ \mathbf{G} }[/math] and [math]\displaystyle{ \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]\displaystyle{ \mathbf{G} }[/math] and [math]\displaystyle{ \mathbf{G}_{n} }[/math] can be found in Boundary Element Method.

The outward normal derivative of the potential, [math]\displaystyle{ \vec{\phi}_{n}, }[/math] and the potential, [math]\displaystyle{ \vec{\phi}, }[/math] are related by the conditions on the boundary [math]\displaystyle{ \partial \Omega }[/math]. This can be expressed as

[math]\displaystyle{ \vec{\phi}_{n}=\mathbf{A}\,\vec{\phi}-\vec{f}, }[/math]

where [math]\displaystyle{ \mathbf{A} }[/math] is the block diagonal matrix given by

[math]\displaystyle{ \mathbf{A}\mathbb{=}\left[ \begin{matrix} \mathbf{Q} & & & \\ & \mathbf{Q} & & \\ & & \nu \,\mathbf{I} & \\ & & & 0 \end{matrix} \right] , }[/math]

where we have separated the boundaries to the semi-infinite domains, the free surface and the sea floor and body boundary. [math]\displaystyle{ \mathbf{Q} }[/math] is a matrix approximation of the integral operator of the same name and [math]\displaystyle{ \vec{f} }[/math] is the vector

[math]\displaystyle{ \vec{f}=\left[ \begin{matrix} 2k_{0}\phi_0 \,e^{k_{0}l } \\ 0 \\ \vdots \\ 0 \end{matrix} \right] . }[/math]

We obtain the following matrix equation for the potential

[math]\displaystyle{ \left( \frac{1}{2}\ \mathbf{I}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G} \vec{f} }[/math]

which can be solved straightforwardly. The reflection and transmission coefficients are calculated from the solution.

Numerical Calculation of [math]\displaystyle{ \mathbf{Q} }[/math]

We begin by truncating to a finite number ([math]\displaystyle{ N }[/math]) of evanescent modes,

[math]\displaystyle{ \mathbf{Q}\phi =\sum_{m=0}^{N}k_{m}\int_{-h}^{0} \tilde{\phi} \left( s\right) \phi_{m}\left( s\right) \mathrm{d}s \frac{\phi _{m} \left( z\right)}{A_m} . }[/math]

We calculate the integral with the same panels as we used to approximate the integrals of the Green function and its normal derivative . Similarly, we assume that [math]\displaystyle{ \phi }[/math] is constant over each panel and integrate [math]\displaystyle{ \phi _{m}\left( z\right) }[/math] exactly. This gives us the following matrix factorisation of [math]\displaystyle{ \mathbf{Q} }[/math]

[math]\displaystyle{ \mathbf{Q}\vec{\phi}=\mathbf{S}\,\mathbf{R}\,\vec{\phi}. }[/math]

The components of the matrices [math]\displaystyle{ \mathbf{S} }[/math] and [math]\displaystyle{ \mathbf{R} }[/math] are

[math]\displaystyle{ \begin{matrix} s_{im} &=&\frac{\phi _{m}^{\left( 1\right) }\left( z_{i}\right)}{A_m} , \\ r_{mj} &=&k_{m} \int_{z_{j}-h/2}^{z_{j}+h/2}\phi _{m}\left( s\right) \mathrm{d}s \end{matrix} }[/math]

where [math]\displaystyle{ z_{j} }[/math] is the value of the [math]\displaystyle{ z }[/math] coordinate in the centre of the panel and [math]\displaystyle{ h }[/math] is the panel length. The integral operator [math]\displaystyle{ \mathbf{Q} }[/math] is approximated by a matrix in exactly the same manner.

Matlab Code

A program to solve this problem can be found here fixed_body_twod.m

Additional code

This program requires