Resonance in moonpool type structures

I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% end

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n

end

for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G

G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %

end

plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
endI have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

end

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

end

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end moonpool, resonance, system of linear equation, semi analytical solution, wave motion MATLAB Answers — New Questions

Cart

Cart