Eigenfunction Matching for a Circular Floating Elastic Plate

Introduction

We show here a solution for a Floating Elastic Plate on Finite Depth water based on Peter, Meylan and Chung 2004. A solution for Shallow Depth was given in Zilman and Miloh 2000 and we will also show this. The solution is an extension of the Eigenfunction Matching for a Circular Dock.

Governing Equations

We begin with the Frequency Domain Problem for a Floating Elastic Plate in the non-dimensional form of Tayler 1986 (Dispersion Relation for a Floating Elastic Plate) We will use a cylindrical coordinate system, $(r,\theta,z)$, assumed to have its origin at the centre of the circular plate which has radius $a$. The water is assumed to have constant finite depth $H$ and the $z$-direction points vertically upward with the water surface at $z=0$ and the sea floor at $z=-H$. The boundary value problem can therefore be expressed as

$\Delta\phi=0, \,\, -h

$\phi_{z}=0, \,\, z=-h,$

$\phi_{z}=\alpha\phi, \,\, z=0,\,r>a,$

$(\beta\Delta^{2}+1-\alpha\gamma)\phi_{z}=\alpha\phi, \,\, z=0,\,r

where the constants $\beta$ and $\gamma$ are given by

$\beta=\frac{D}{\rho\,L^{4}g}, \gamma=\frac{\rho_{i}h}{\rho\,L}$

and $\rho_{i}$ is the density of the plate. We must also apply the edge conditions for the plate and the Sommerfeld Radiation Condition as $r\rightarrow\infty$. The subscript $z$ denotes the derivative in $z$-direction.

Solution Method

We use separation of variables in the two regions, $r and $r>a$.

The solution of the problem for the potential in finite water depth can be found by a separation ansatz,

$\phi (r,\theta,z) =: Y(r,\theta) Z(z).\,$

Substituting this into the equation for $\phi$ yields

$\frac{1}{Y(r,\theta)} \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial Y}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 Y}{\partial \theta^2} \right] = - \frac{1}{Z(z)} \frac{\mathrm{d}^2 Z}{\mathrm{d} z^2} = k^2.$

The possible separation constants $k$ will be determined by the free surface condition and the bed condition.

Separation of variables for a free surface

We use separation of variables

We express the potential as

$\phi(x,z) = X(x)Z(z)\,$

and then Laplace's equation becomes

$\frac{X^{\prime\prime}}{X} = - \frac{Z^{\prime\prime}}{Z} = k^2$

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

$Z^{\prime\prime} + k^2 Z =0.$

subject to the boundary conditions

$Z^{\prime}(-h) = 0$

and

$Z^{\prime}(0) = \alpha Z(0)$

We can then use the boundary condition at $z=-h \,$ to write

$Z = \frac{\cos k(z+h)}{\cos kh}$

where we have chosen the value of the coefficent so we have unit value at $z=0$. The boundary condition at the free surface ($z=0 \,$) gives rise to:

$k\tan\left( kh\right) =-\alpha \,$

which is the Dispersion Relation for a Free Surface

The above equation is not really the dispersion relation for a free surface, it would be better to refer to it as 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 $k_{0}=\pm ik \,$ and the positive real solutions by $k_{m} \,$, $m\geq1$. The $k \,$ of the imaginary solution is the wavenumber. We put the imaginary roots back into the equation above and use the hyperbolic relations

$\cos ix = \cosh x, \quad \sin ix = i\sinh x,$

to arrive at the dispersion relation

$\alpha = k\tanh kh.$

We note that for a specified frequency $\omega \,$ the equation determines the wavenumber $k \,$.

Finally we define the function $Z(z) \,$ as

$\chi_{m}\left( z\right) =\frac{\cos k_{m}(z+h)}{\cos k_{m}h},\quad m\geq0$

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

$\int\nolimits_{-h}^{0}\chi_{m}(z)\chi_{n}(z) \mathrm{d} z=A_{n}\delta_{mn}$

where

$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).$

Separation of variables under the Plate

$Z^{\prime\prime} + \kappa^2 Z =0.$

subject to the boundary conditions

$Z^{\prime}(-h) = 0$

and

$\left(\beta \kappa^4 + 1 - \alpha\gamma\right)Z^{\prime}(0) = \alpha Z(0)$

(the first term comes from the beam eigenvalue problem, where $\partial_x^4 X = \kappa^4 X$). We then use the boundary condition at $z=-h \,$ to write

$Z = \frac{\cos \kappa(z+h)}{\cos \kappa h}$

The boundary condition at the free surface ($z=0$) is the Dispersion Relation for a Floating Elastic Plate

$\kappa \tan{(\kappa h)}= -\frac{\alpha}{\beta \kappa^{4} + 1 - \alpha\gamma}$

Solving for $\kappa \,$ gives a pure imaginary root with positive imaginary part, two complex roots (two complex conjugate paired roots with positive imaginary part in most physical situations), an infinite number of positive real roots which approach ${n\pi}/{h} \,$ as $n$ approaches infinity, and also the negative of all these roots (Dispersion Relation for a Floating Elastic Plate) . We denote the two complex roots with positive imaginary part by $\kappa_{-2} \,$ and $\kappa_{-1} \,$, the purely imaginary root with positive imaginary part by $\kappa_{0} \,$ and the real roots with positive imaginary part by $\kappa_{n} \,$ for $n$ a positive integer. The imaginary root with positive imaginary part corresponds to a reflected travelling mode propagating along the $x$ axis. The complex roots with positive imaginary parts correspond to damped reflected travelling modes and the real roots correspond to reflected evanescent modes.

Inner product between free surface and elastic plate modes

$\int\nolimits_{-h}^{0}\phi_{n}(z)\psi_{m}(z) d z=B_{mn}$

where

$B_{mn}=\frac{k_{n}\sin k_{n}h\cos\kappa_{m}h-\kappa_{m}\cos k_{n}h\sin \kappa_{m}h}{\left( \cos k_{n}h\right) \left( \cos \kappa_{m}h\right) \left( k_{n} ^{2}-\kappa_{m}^{2}\right) }$

Separation for Cylindrical Coordinates

We now separate variables, noting that since the problem has circular symmetry we can write the potential as

$\phi(r,\theta,z)=\frac{\cos k(z+h)}{\cos kh}\sum_{n=-\infty}^{\infty}\rho_{n}(r)e^{i n \theta}$

We now solve for the function $\rho_{n}(r)$. Using Laplace's equation in polar coordinates we obtain

$\frac{\mathrm{d}^{2}\rho_{n}}{\mathrm{d}r^{2}}+\frac{1}{r} \frac{\mathrm{d}\rho_{n}}{\mathrm{d}r}-\left( \frac{n^{2}}{r^{2}}+k^{2}\right) \rho_{n}=0.$

We can convert this equation to the standard form by substituting $y=k r$ (provided that $\mu\neq 0$to obtain

$y^{2}\frac{\mathrm{d}^{2}\rho_{n}}{\mathrm{d}y^{2}}+y\frac{\mathrm{d}\rho_{n} }{\rm{d}y}-(n^{2}+y^{2})\rho_{n}=0$

The solution of this equation is a linear combination of the modified Bessel functions of order $n$, $I_{n}(y)$ and $K_{n}(y)$ (Abramowitz and Stegun 1964).

Therefore

$\rho_n(r) = C_1 I_{n}(kr) + C_2 K_{n}(kr)\,$

for some constants $C_1$ and $C_2$

Incident potential

The incident potential is a wave of amplitude $A$ in displacement travelling in the positive $x$-direction. The incident potential can therefore be written as

$\phi^{\mathrm{I}} =e^{k_{0}x}\phi_{0}\left( z\right) =\sum\limits_{n=-\infty}^{\infty} I_{n}(k_{0}r)\phi_{0}\left(z\right) e^{i n \theta}$

Expansion of the potential

Since the solution must be bounded we know that under the plate the solution will be a linear combination of $I_{n}(y)$ while outside the plate the solution will be a linear combination of $K_{n}(y)$. Therefore the potential can be expanded as

$\phi(r,\theta,z)=\sum_{n=-\infty}^{\infty}\sum_{m=0}^{\infty}a_{mn}K_{n} (k_{m}r)e^{i n\theta}\phi_{m}(z), \;\;r>a$

and

$\phi(r,\theta,z)=\sum_{n=-\infty}^{\infty}\sum_{m=-2}^{\infty}b_{mn} I_{n}(\kappa_{m}r)e^{i n\theta}\psi_{m}(z), \;\;r

where $a_{mn}$ and $b_{mn}$ are the coefficients of the potential in the open water and the plate covered region respectively.

Boundary conditions

The boundary conditions for the plate also have to be considered. The vertical force and bending moment must vanish, which can be written as

$\left[\bar{\Delta}-\frac{1-\nu}{r}\left(\frac{\partial}{\partial r} +\frac{1}{r}\frac{\partial^{2}}{\partial\theta^{2}}\right)\right] w=0\,\,\,(3)$

and

$\left[ \frac{\partial}{\partial r}\bar{\Delta}-\frac{1-\nu}{r^{2}}\left( -\frac{\partial}{\partial r}+\frac{1}{r}\right) \frac{\partial^{2}} {\partial\theta^{2}}\right] w=0 \,\,\,(4)$

where $w$ is the time-independent surface displacement, $\nu$ is Poisson's ratio, and $\bar{\Delta}$ is the polar coordinate Laplacian

$\bar{\Delta}=\frac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\theta^{2}}$

Displacement of the plate

The surface displacement and the water velocity potential at the water surface are linked through the kinematic boundary condition

$\phi_{z}=-i\sqrt{\alpha}w,\,\,\,z=0$

From the equations the potential and the surface displacement are therefore related by

$w=i\sqrt{\alpha}\phi,\quad r>a$

and

$(\beta\bar{\Delta}^{2}+1-\alpha\gamma)w=i\sqrt{\alpha}\phi,\quad r

The surface displacement can also be expanded in eigenfunctions as

$w(r,\theta)=\sum_{n=-\infty}^{\infty}\sum_{m=0}^{\infty}i\sqrt{\alpha} a_{mn}K_{n}(k_{m}r)e^{i n\theta},\;\;r>a$

and:

$w(r,\theta)= \sum_{n=-\infty}^{\infty}\sum_{m=-2}^{\infty}i\sqrt{\alpha}(\beta\kappa _{m}^{4}+1-\alpha\gamma)^{-1}b_{mn}I_{n}(\kappa_{m}r)e^{i n\theta},\; r

using the fact that

$\bar{\Delta}\left( I_{n}(\kappa_{m}r)e^{i n\theta}\right) =\kappa_{m} ^{2}I_{n}(\kappa_{m}r)e^{i n\theta}\,\,\,(5)$

An infinite dimensional system of equations

The boundary conditions (3) and (4) can be expressed in terms of the potential using (5). Since the angular modes are uncoupled the conditions apply to each mode, giving

$\sum_{m=-2}^{\infty}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left(\kappa_{m}^{2}I_{n}(\kappa_{m}a)-\frac{1-\nu}{a}\left(\kappa _{m}I_{n}^{\prime}(\kappa_{m}a)-\frac{n^{2}}{a}I_{n}(\kappa_{m}a)\right) \right) =0\,\,\,(6)$

$\sum_{m=-2}^{\infty}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left( \kappa_{m}^{3}I_{n}^{\prime}(\kappa_{m}a)+n^{2}\frac{1-\nu}{a^{2} }\left( -\kappa_{m}I_{n}^{\prime}(\kappa_{m}a)+\frac{1}{a}I_{n}(\kappa _{m}a)\right) \right) =0\,\,\,(7)$

The potential and its derivative must be continuous across the transition from open water to the plate covered region. Therefore, the potentials and their derivatives at $r=a$ have to be equal. Again we know that this must be true for each angle and we obtain

$e_{n}I_{n}(k_{0}a)\phi_{0}\left( z\right) + \sum_{m=0}^{\infty} a_{mn} K_{n}(k_{m}a)\phi_{m}\left( z\right) =\sum_{m=-2}^{\infty}b_{mn}I_{n}(\kappa_{m}a)\psi_{m}(z)$

and

$e_{n}k_{0}I_{n}^{\prime}(k_{0}a)\phi_{0}\left( z\right) +\sum _{m=0}^{\infty} a_{mn}k_{m}K_{n}^{\prime}(k_{m}a)\phi_{m}\left( z\right) =\sum_{m=-2}^{\infty}b_{mn}\kappa_{m}I_{n}^{\prime}(\kappa_{m}a)\psi _{m}(z)$

for each $n$. We solve these equations by multiplying both equations by $\phi_{l}(z)$ and integrating from $-H$ to $0$ to obtain:

$e_{n}I_{n}(k_{0}a)A_{0}\delta_{0l}+a_{ln}K_{n}(k_{l}a)A_{l} =\sum_{m=-2}^{\infty}b_{mn}I_{n}(\kappa_{m}a)B_{ml} \,\,\,(8)$

and

$e_{n}k_{0}I_{n}^{\prime}(k_{0}a)A_{0}\delta_{0l}+a_{ln}k_{l}K_{n}^{\prime }(k_{l}a)A_{l} =\sum_{m=-2}^{\infty}b_{mn}\kappa_{m}I_{n}^{\prime}(\kappa_{m} a)B_{ml} \,\,\,(9)$

Equation (8) can be solved for the open water coefficients $a_{mn}$

$a_{ln}=-e_{n}\frac{I_{n}(k_{0}a)}{K_{n}(k_{0}a)}\delta_{0l}+\sum _{m=-2}^{\infty}b_{mn}\frac{I_{n}(\kappa_{m}a)B_{ml}}{K_{n}(k_{l}a)A_{l}}$

which can then be substituted into equation (9) to give us

$\left( k_{0}I_{n}^{\prime}(k_{0}a)-k_{0}\frac{K_{n}^{\prime}(k_{0} a)}{K_{n}(k_{0}a)}I_{n}(k_{0}a)\right) e_{n}A_{0}\delta_{0l} =\sum_{m=-2}^{\infty}\left( \kappa_{m}I_{n}^{\prime}(\kappa_{m} a)-k_{l}\frac{K_{n}^{\prime}(k_{l}a)}{K_{n}(k_{l}a)}I_{n}(\kappa _{m}a)\right) B_{ml}b_{mn}\,\,\,(10)$

for each $n$. Together with equations (6) and (7) equation (10) gives the required equations to solve for the coefficients of the water velocity potential in the plate covered region.

Numerical Solution

To solve the system of equations (10) together with the boundary conditions (6 and 7) we set the upper limit of $l$ to be $M$. We also set the angular expansion to be from $n=-N$ to $N$. This gives us

$\phi(r,\theta,z)=\sum_{n=-N}^{N}\sum_{m=0}^{M}a_{mn}K_{n}(k_{m}r)e^{i n\theta }\phi_{m}(z), \;\;r>a$

and

$\phi(r,\theta,z)=\sum_{n=-N}^{N}\sum_{m=-2}^{M}b_{mn}I_{n}(\kappa _{m}r)e^{i n\theta}\psi_{m}(z), \;\;r

Since $l$ is an integer with $0\leq l\leq M$ this leads to a system of $M+1$ equations. The number of unknowns is $M+3$ and the two extra equations are obtained from the boundary conditions for the free plate (6) and (7). The equations to be solved for each $n$ are

$\left( k_{0}I_{n}^{\prime}(k_{0}a)-k_{0}\frac{K_{n}^{\prime}(k_{0} a)}{K_{n}(k_{0}a)}I_{n}(k_{0}a)\right) e_{n}A_{0}\delta_{0l} =\sum_{m=-2}^{M}\left( \kappa_{m}I_{n}^{\prime}(\kappa_{m}a)-k_{l} \frac{K_{n}^{\prime}(k_{l}a)}{K_{n}(k_{l}a)}I_{n}(\kappa_{m}a)\right) B_{ml}b_{mn}$

$\sum_{m=-2}^{M}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left( \kappa_{m}^{2}I_{n}(\kappa_{m}a)-\frac{1-\nu}{a}\left( \kappa _{m}I_{n}^{\prime}(\kappa_{m}a)-\frac{n^{2}}{a}I_{n}(\kappa_{m}a)\right) \right) =0$

and

$\sum_{m=-2}^{M}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left( \kappa_{m}^{3}I_{n}^{\prime}(\kappa_{m}a)+n^{2}\frac{1-\nu}{a^{2} }\left( \kappa_{m}I_{n}^{\prime}(\kappa_{m}a)+\frac{1}{a}I_{n}(\kappa _{m}a)\right) \right) =0$

It should be noted that the solutions for positive and negative $n$ are identical so that they do not both need to be calculated. There are some minor simplifications which are a consequence of this which are discussed in more detail in Zilman and Miloh 2000.

The Shallow Depth Theory of Zilman and Miloh 2000

The shallow water theory of Zilman and Miloh 2000 can be recovered by simply setting the depth shallow enough that the shallow water theory is valid and setting $M=0$. If the shallow water theory is valid then the first three roots of the dispersion equation for the ice will be exactly the same roots found in the shallow water theory by solving the polynomial equation. The system of equations has four unknowns (three under the plate and one in the open water) exactly as for the theory of Zilman and Miloh 2000.

A Simple Method To Calculate The Diffraction Transfer Matrix For The Case Of A Circular Plate

Let's consider an incident wave whose potential has the following expression

$\phi^\mathrm{I} (r,\theta,z) = \sum_{n=0}^{\infty} \phi_n(z) \sum_{\nu = - \infty}^{\infty} D_{n\nu} I_\nu (k_n r) \mathrm{e}^{\mathrm{i}\nu \theta}.$

Such an incident potential is found in the Kagemoto and Yue Interaction Theory, where it can be written as the sum of an ambient incident potential and the scattered potentials of the other bodies, which are interpretated as incident potentials for the studied body.

We can apply the same eigenfunction matching that previously, considering the potential and its normal derivative continuous at $r=a$. Thus the potential and its normal derivative expressed at each side of this value of the radius have to be equal. We obtain the relationships

$\sum_{m=0}^{\infty} D_{mn} I_n (k_m a) \phi_m(z) + \sum_{m=0}^{\infty} a_{mn} K_{n}(k_{m}a)\phi_{m}\left( z\right) =\sum_{m=-2}^{\infty}b_{mn}I_{n}(\kappa_{m}a)\psi_{m}(z)$

and

$\sum_{m=0}^{\infty} D_{mn} k_m I'_n (k_m a) \phi_m(z) + \sum_{m=0}^{\infty} a_{mn} k_m K'_{n}(k_{m}a)\phi_{m}\left( z\right) =\sum_{m=-2}^{\infty}b_{mn} \kappa_m I'_{n}(\kappa_{m}a)\psi_{m}(z)$

for each $n$. We solve these equations with the same method that before, by multiplying both equations by $\phi_{l}(z)$ and integrating from $-H$ to $0$ to obtain:

$D_{ln} I_n (k_l a) A_l + a_{ln} K_{n}(k_l a) A_l = \sum_{m=-2}^{\infty} b_{mn} I_{n}(\kappa_{m} a) B_{ml},\ \ \ (11)$

and

$D_{ln} k_l I'_n (k_l a) A_l + a_{ln} k_l K'_{n}(k_{l} a) A_l =\sum_{m=-2}^{\infty}b_{mn} \kappa_m I'_{n}(\kappa_{m}a) B_{ml} \ \ \ (12)$

The Diffraction Transfer Matrix maps the coefficients of the incident wave with the coefficients of the scattered wave within the open water domain. The relation which links these two coefficients can be written as follows

$a_{mn}=\sum_{l=0}^{\infty} T_{lmn} D_{ln}$

Furthermore the boundary conditions are exactly the same that before, namely

$\sum_{m=-2}^{\infty}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left(\kappa_{m}^{2}I_{n}(\kappa_{m}a)-\frac{1-\nu}{a}\left(\kappa _{m}I_{n}^{\prime}(\kappa_{m}a)-\frac{n^{2}}{a}I_{n}(\kappa_{m}a)\right) \right) =0$

$\sum_{m=-2}^{\infty}(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1}b_{mn} \left( \kappa_{m}^{3}I_{n}^{\prime}(\kappa_{m}a)+n^{2}\frac{1-\nu}{a^{2} }\left( -\kappa_{m}I_{n}^{\prime}(\kappa_{m}a)+\frac{1}{a}I_{n}(\kappa _{m}a)\right) \right) =0 .$

For the further study, let's call

$L^1_{mn}=(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1} \left(\kappa_{m}^{2}I_{n}(\kappa_{m}a)-\frac{1-\nu}{a}\left(\kappa _{m}I_{n}^{\prime}(\kappa_{m}a)-\frac{n^{2}}{a}I_{n}(\kappa_{m}a)\right) \right)$

and

$L^2_{mn}=(\beta\kappa_{m}^{4}+1-\alpha\gamma)^{-1} \left( \kappa_{m}^{3}I_{n}^{\prime}(\kappa_{m}a)+n^{2}\frac{1-\nu}{a^{2} }\left( -\kappa_{m}I_{n}^{\prime}(\kappa_{m}a)+\frac{1}{a}I_{n}(\kappa _{m}a)\right) \right)$

From the equations (11), (12) and the boundary conditions over the edges of the plate, we can write a linear system of equation, limiting the number of modes of the dispersion equation to $N$ real ones

$\begin{bmatrix} \begin{bmatrix} -A_0 K_n(k_0 a)&0 \quad \cdots&0\\ 0&&\\ \vdots&-A_l K_n(k_l a)&\vdots\\ &&0\\ 0&\cdots \quad 0 & -A_N K_n(k_N a) \end{bmatrix} & \begin{bmatrix} I_n(\kappa_{-2}a) B_{-20}&\cdots&I_n(\kappa_{N}a) B_{N0}\\ &&\\ \vdots&&\vdots\\ &&\\ I_n(\kappa_{-2}a) B_{-2N}&\cdots&I_n(\kappa_{N}a) B_{NN} \end{bmatrix} \\ \begin{bmatrix} 0&&\cdots&&0\\ 0&&\cdots&&0 \end{bmatrix} & \begin{bmatrix} L^1_{-2n}&\cdots&L^1_{Nn}\\ L^2_{-2n}&\cdots&L^2_{Nn} \end{bmatrix} \\ \begin{bmatrix} -k_0 K'_n(k_0 a) A_0&0 \quad \cdots&0\\ 0&&\\ \vdots&-k_l K'_n(k_l a) A_l&\vdots\\ &&0\\ 0&\cdots \quad 0 & -k_N K'_n(k_N a) A_N \end{bmatrix} & \begin{bmatrix} \kappa_{-2} I'_n(\kappa_{-2}a) B_{-20}&\cdots&\kappa_{N} I'_n(\kappa_{N}a) B_{N0}\\ &&\\ \vdots&&\vdots\\ &&\\ \kappa_{-2} I'_n(\kappa_{-2}a) B_{-2N}&\cdots&\kappa_{N} I'_n(\kappa_{N}a) B_{NN} \end{bmatrix} \end{bmatrix} \begin{bmatrix} a_{0n} \\ \\ \vdots \\ \\ a_{Nn} \\ b_{-2n}\\ b_{-1n}\\ b_{0n}\\ \\ \vdots \\ \\ b_{Nn} \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} I_n(k_0 a) A_0&0 \quad \cdots&0\\ 0&&\\ \vdots&I_n(k_l a) A_l&\vdots\\ &&0\\ 0&\cdots \quad 0 & I_n(k_N a) A_N \end{bmatrix} \\ \begin{bmatrix} 0&&\cdots&&0\\ 0&&\cdots&&0 \end{bmatrix} \\ \begin{bmatrix} k_0 I'_n(k_0 a) A_0&0 \quad \cdots&0\\ 0&&\\ \vdots&k_l I'_n(k_l a) A_l&\vdots\\ &&0\\ 0&\cdots \quad 0 &k_N I'_n(k_N a) A_N \end{bmatrix} \end{bmatrix} \begin{bmatrix} D_{0n} \\ \\ \vdots \\ \\ D_{Nn} \\ 0\\ 0\\ D_{0n}\\ \\ \vdots \\ \\ D_{Nn} \end{bmatrix}$

for each $n$.

Therefore we can find a Diffraction Transfer Matrix for each $n$, by setting

$\forall i \in [0, N], (D_{pn})_{p \in [0, N]} = \delta_{ip}$

Then we solve the linear system defined previously, so that we can find the coefficients $(a_{ln})_{l \in [0, N]}$ for each $i$. This vector represents exactly the $i^{th}$ column of the Diffraction Transfer Matrix, $n$ being set.

This method permits to obtain the matrix which links the coefficients of the incident and scattered potential in the free water domain. Applying this for each $n$, we finally obtain a 3-dimensional matrix for the Diffraction Transfer Matrix.

Matlab Code

A program to calculate the coefficients for circular plate problems can be found here circle_plate_matching_one_n.m Note that this problem solves only for a single n.