P3_水面方箱对波浪的反射和透射

W_星烨

已于 2024-05-17 20:43:26 修改

阅读量285

点赞数 10

分类专栏： A01_波浪对海上建筑物的作用文章标签：线性代数算法矩阵

于 2024-05-17 20:38:52 首次发布

本文链接：https://blog.csdn.net/w1418169121/article/details/138921103

版权

A01_波浪对海上建筑物的作用专栏收录该内容

3 篇文章 0 订阅

订阅专栏

1、简介

右向波入射时，可简化为对称波入射和反对称波入射的叠加情况，叠加后，左向波因相位不同而抵消掉，右向波相当于两次入射，故而速度势要除以2，如下

此时得到可得到对称、反对称问题中的反射系数Rs和Ra，用其可以求出反射系数R。

依然可以从给出的四个方程组联合求解，稍微麻烦点的是系数的计算，以及计算波浪力时，速度势中求和部分的积分计算，我是逐项计算再求和，没想什么其他的方法。

注意下图中的速度势，对称势和反对称势是分段函数，别看错了分段，Fx应该用φ2a计算，之前误用φ1a计算区别还挺大的。

2、方程建模及求解

clc,clear,%close all
d=20;
T_all=[0.25*d,0.5*d,0.75*d];
w_all= 0.05:0.025:1.4 ;
A=1;
g=9.81;
m=1;n=m;

for T_n=1:size(T_all,2)
    T=T_all(T_n);
    S=d-T;
    B=d;

    for k=1:size(w_all,2)
        w=w_all(k);
        k0=w_k(w,d);
        lambda0=0;
        km=dispersion(k0,d,m);
        lambdan=pi/S*(1:n);

        %求系数f
        f=zeros(n+1,1);
        %Z0=@(z)cosh(k0*(z+d))/cosh(k0*d);
        Y0=sqrt(2)/2;
        Z0Y0=@(z)cosh(k0*(z+d))/cosh(k0*d)*Y0;
        f0_val=2/S*integral(Z0Y0,-d,-T);
        f1_val=zeros(n,1);
        for i=1:n
            %Yn=@(z)cos(lambdan(i)*(z+d));
            Z0Yn=@(z)cosh(k0*(z+d))/cosh(k0*d).*cos(lambdan(i)*(z+d));
            f1_val(i)=2/S*integral(Z0Yn,-d,-T);
        end
        f(1,1)=f0_val;f(2:n+1,1)=f1_val;

        %求系数矩阵A
        a=zeros(n+1,m+1);
        a0m_val=zeros(1,m);
        for i=1:m
            %Zm=@(z)cos(km(i)*(z+d))/cos(km(i)*d);
            ZmY0=@(z)cos(km(i)*(z+d))/cos(km(i)*d)*Y0;
            a0m_val(i)=2/S*integral(ZmY0,-d,-T);
        end

        a_val=zeros(n,m);
        for i=1:n
            for j=1:m
                %Zm=@(z)cos(km(j)*(z+d))/cos(km(j)*d);
                %Yn=@(z)cos(lambdan(i)*(z+d));
                ZmYn=@(z)cos(km(j)*(z+d))/cos(km(j)*d).*cos(lambdan(i)*(z+d));
                a_val(i,j)=2/S*integral(ZmYn,-d,-T);
            end
        end
        a(:,1)=f;
        a(1,2:m+1)=a0m_val;
        a(2:n+1,2:m+1)=a_val;

        %求第二个方程组矩阵系数
        b0n=zeros(1,n);
        for i=1:n
            %Yn=@(z)cos(lambdan(i)*(z+d));
            Z0Yn=@(z)cosh(k0*(z+d))/cosh(k0*d).*cos(lambdan(i)*(z+d));
            Z0Z0=@(z)cosh(k0*(z+d))/cosh(k0*d).*cosh(k0*(z+d))/cosh(k0*d);
            b0n(i)=(-lambdan(i)/k0*tanh(lambdan(i)*B)*integral(Z0Yn,-d,-T)/integral(Z0Z0,-d,0));
        end

        bm0=zeros(m+1,1);

        bmn=zeros(m,n);
        for i=1:m
            for j=1:n
                %Zm=@(z)cos(km(i)*(z+d))/cos(km(i)*d);
                %Yn=@(z)cos(lambdan(i)*(z+d));
                ZmYn=@(z)cos(km(i)*(z+d))/cos(km(i)*d).*cos(lambdan(j)*(z+d));
                ZmZm=@(z)(cos(km(i)*(z+d))/cos(km(i)*d)).^2;
                bmn(i,j)=-lambdan(j)/km(i)*tanh(lambdan(j)*B)*integral(ZmYn,-d,-T)/integral(ZmZm,-d,0);
            end
        end

        b=zeros(m+1,n+1);
        b(1,2:n+1)=b0n*1i;
        b(1:m+1,1)=0;
        b(2:m+1,2:n+1)=bmn;

        %c1=w/k0*0.5*(1+2*k0*d/sinh(2*k0*d));
        %c2=w/lambda0*0.5*(1+2*lambda0*d2/sinh(2*lambda0*d2));

        e=zeros(m+1,1);
        e(1)=-1;

        y=[f;e];
        AA=[-1*a,eye(m+1);(-1*eye(m+1)),b];
        x=AA\y;
        x_r=real(x);
        x_i=imag(x);
        Asm=x_r(1:m+1);
        Bsn=x_r(m+2:end);
        AAA_s(k,1+5*(T_n-1))=w;
        AAA_s(k,2+5*(T_n-1))=k0;
        AAA_s(k,3+5*(T_n-1))=Asm(1,1);
        AAA_s(k,4+5*(T_n-1))=Bsn(1,1);
        AAA_s(k,5+5*(T_n-1))=T;
        igA=(1i*g*A/w);
        Fzm(1)=1i*w*B*(-igA) * (Bsn(1)*Y0);
        for i_fz=1:m
            phiasm=@(x,z)(-igA)*( Bsn(1+i_fz)*cosh(lambdan(i_fz)*x)/cosh(lambdan(i_fz)*B) )*cos(lambdan(i_fz)*(z+d));
            Fzm(1+i_fz)=1i*w*integral(@(x)phiasm(x,-T),-B,0);
        end
        Fz(k,T_n)=(sum(Fzm));
    end
end


%求反对称系数
for T_n=1:size(T_all,2)
    T=T_all(T_n);
    S=d-T;
    B=d;

    for k=1:size(w_all,2)
        w=w_all(k);
        k0=w_k(w,d);
        lambda0=0;
        km=dispersion(k0,d,m);
        lambdan=pi/S*(1:n);

        %求系数f
        f=zeros(n+1,1);
        %Z0=@(z)cosh(k0*(z+d))/cosh(k0*d);
        Y0=sqrt(2)/2;
        Z0Y0=@(z)cosh(k0*(z+d))/cosh(k0*d)*Y0;
        f0_val=-2/S*integral(Z0Y0,-d,-T);
        f1_val=zeros(n,1);
        for i=1:n
            %Yn=@(z)cos(lambdan(i)*(z+d));
            Z0Yn=@(z)cosh(k0*(z+d))/cosh(k0*d).*cos(lambdan(i)*(z+d));
            f1_val(i)=-2/S*integral(Z0Yn,-d,-T);
        end
        f(1,1)=f0_val;f(2:n+1,1)=f1_val;

        %求系数矩阵A
        a=zeros(n+1,m+1);
        a0m_val=zeros(1,m);
        for i=1:m
            %Zm=@(z)cos(km(i)*(z+d))/cos(km(i)*d);
            ZmY0=@(z)cos(km(i)*(z+d))/cos(km(i)*d)*Y0;
            a0m_val(i)=-2/S*integral(ZmY0,-d,-T);
        end

        a_val=zeros(n,m);
        for i=1:n
            for j=1:m
                %Zm=@(z)cos(km(j)*(z+d))/cos(km(j)*d);
                %Yn=@(z)cos(lambdan(i)*(z+d));
                ZmYn=@(z)cos(km(j)*(z+d))/cos(km(j)*d).*cos(lambdan(i)*(z+d));
                a_val(i,j)=-2/S*integral(ZmYn,-d,-T);
            end
        end
        a(:,1)=f;
        a(1,2:m+1)=a0m_val;
        a(2:n+1,2:m+1)=a_val;

        %求第二个方程组矩阵系数
        Z0Y0=@(z)cosh(k0*(z+d))/cosh(k0*d)*Y0;
        Z0Z0=@(z)cosh(k0*(z+d))/cosh(k0*d).*cosh(k0*(z+d))/cosh(k0*d);
        b00=1/k0/B*integral(Z0Y0,-d,-T)/integral(Z0Z0,-d,0);

        b0n=zeros(1,n);
        for i=1:n
            %Yn=@(z)cos(lambdan(i)*(z+d));
            Z0Yn=@(z)cosh(k0*(z+d))/cosh(k0*d).*cos(lambdan(i)*(z+d));
            Z0Z0=@(z)cosh(k0*(z+d))/cosh(k0*d).*cosh(k0*(z+d))/cosh(k0*d);
            b0n(i)=(lambdan(i)/k0/tanh(lambdan(i)*B)*integral(Z0Yn,-d,-T)/integral(Z0Z0,-d,0));
        end

        bm0=zeros(m,1);
        for i=1:m
            %Yn=@(z)cos(lambdan(i)*(z+d));
            ZmY0=@(z)cos(km(i)*(z+d))/cos(km(i)*d).*Y0;
            ZmZm=@(z)(cos(km(i)*(z+d))/cos(km(i)*d)).^2;
            bm0(i)=(1/km(i)/B*integral(Z0Yn,-d,-T)/integral(Z0Z0,-d,0));
        end


        bmn=zeros(m,n);
        for i=1:m
            for j=1:n
                %Zm=@(z)cos(km(i)*(z+d))/cos(km(i)*d);
                %Yn=@(z)cos(lambdan(i)*(z+d));
                ZmYn=@(z)cos(km(i)*(z+d))/cos(km(i)*d).*cos(lambdan(j)*(z+d));
                ZmZm=@(z)(cos(km(i)*(z+d))/cos(km(i)*d)).^2;
                bmn(i,j)=lambdan(j)/km(i)/tanh(lambdan(j)*B)*integral(ZmYn,-d,-T)/integral(ZmZm,-d,0);
            end
        end

        b=zeros(m+1,n+1);
        b(1,1)=b00*1i;
        b(1,2:n+1)=b0n*1i;
        b(2:m+1,1)=bm0;
        b(2:m+1,2:n+1)=bmn;

        %c1=w/k0*0.5*(1+2*k0*d/sinh(2*k0*d));
        %c2=w/lambda0*0.5*(1+2*lambda0*d2/sinh(2*lambda0*d2));

        e=zeros(m+1,1);
        e(1)=-1;

        y=[f;e];
        AA=[-1*a,eye(m+1);(-1*eye(m+1)),b];
        x=AA\y;
        x_r=real(x);
        x_i=imag(x);
        Asm=x_r(1:m+1);
        Bsn=x_r(m+2:end);
        AAA_a(k,1+5*(T_n-1))=w;
        AAA_a(k,2+5*(T_n-1))=k0;
        AAA_a(k,3+5*(T_n-1))=Asm(1,1);
        AAA_a(k,4+5*(T_n-1))=Bsn(1,1);
        AAA_a(k,5+5*(T_n-1))=T;
        phiaa=@(z)(-1i*g*A/w)*( (1+Asm(1))*cosh(k0*(z+d))/cosh(k0*d) );
        Fx(k,T_n)=1i*w*integral(phiaa,-T,0);

    end
end

run B_As1_figure.m

3、计算结果

和书中对比还是存在一定误差，但趋势一致，原因不明；

方箱对长波的反射差，对短波的反射好，透射反之；

水平波浪力先增后减，垂向波浪力一直降低。