矩形方箱垂荡程序.zip

clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% 基本参数定义 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I=50; J=50; for k=1:22 period=4+k*0.5; %定义波周期 H=1.0; %定义波高 OMEGA=0.3+k*0.05; % 定义了作用于结构上的波浪频率 V=22.15*0.5144;%设计航速 beta=pi/6*0;%浪向角 p1=1025; %定义密度 D=100; % 定义水深 T=10; %定义吃水 S=D-T; g=9.82; % 定义重力加速度 L=g*period^2/2/pi; % 定义波长 comp_k=1.4; %定义空气压缩系数 Height=15; %定义结构的高度 HL=100; %定义结构的长度 HB=20; %定义结构的宽度 freeboard=Height-T; %定义结构的干舷高度 WG=4000; HW=WG/HL/HB; %定义浮力中心 P0=101325; %大气压强 density=1.025e3; %定义密度 PCS=P0+HW*density*9.8; %定义筒内气压 VCS=HB*HL*(HW+freeboard); %定义筒内气体体积 miu=OMEGA^2/g; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 定义常数 alpha % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% alpha=comp_k*PCS*HB*HL/density/9.8/VCS; %%%%%%%%%%%%%%%%%%%%% 求解波数K0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fun1=@(x)x.*tanh(D.*x)-miu; K0=abs(fzero(fun1,0,1)); fun2=@(x)x.*tan(x)/D+miu; for n=0:1:300 x1=n*pi+pi/2+0.001; y1=n*pi+pi; aa(n+1)=fzero(fun2,[x1,y1]); end aa=aa/D; %%%%%%%%%%%%%%%%%%%%% 求的所需要的k %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:length(aa)+1 if j==1 kk(j)=K0; else kk(j)=aa(j-1); end end %%%%%%%%%%%%%%%%%%%%% 求解波数B0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fun3=@(y)y.*tanh(S.*y)-miu; B0=abs(fzero(fun3,0,1)); fun4=@(y)y.*tan(y)/S+miu; for n=0:1:300 x2=n*pi+pi/2+0.001; y2=n*pi+pi; cc(n+1)=fzero(fun4,[x2,y2]); end cc=cc/S; %%%%%%%%%%%%%%%%%%%%% 求的所需要的b %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:length(cc)+1 if j==1 bb(j)=B0; else bb(j)=cc(j-1); end end %%%%%%%%%%%%%%%%%%%%% 定义常数N和M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:J if j==1 N(j)=1/cosh(kk(j)*D)^2*(D/2+sinh(2*kk(j)*D)/4/kk(j)); M(j)=1/cosh(bb(j)*S)^2*(S/2+sinh(2*bb(j)*S)/4/bb(j)); else N(j)=1/cos(kk(j)*D)^2*(D/2+sin(2*kk(j)*D)/4/kk(j)); M(j)=1/cos(bb(j)*S)^2*(S/2+sin(2*bb(j)*S)/4/bb(j)); end end %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% 速度连续条件 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% PF\PA\PP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:I for j=1:J if i==1 if j==1 cnmjian(j,i)=1/cosh(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sinh(kk(i)*S)-bb(j)*sinh(bb(j)*S)*cosh(kk(i)*S))/(kk(i)^2-bb(j)^2); else cnmjian(i,j)=1/cosh(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sinh(kk(i)*S)+bb(j)*sin(bb(j)*S)*cosh(kk(i)*S))/(kk(i)^2+bb(j)^2); end else if j==1 cnmjian(j,i)=1/cos(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sin(kk(i)*S)+bb(j)*sinh(bb(j)*S)*cos(kk(i)*S))/(kk(i)^2+bb(j)^2); else cnmjian(i,j)=1/cos(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sin(kk(i)*S)-bb(j)*sin(bb(j)*S)*cos(kk(i)*S))/(kk(i)^2-bb(j)^2); end end end end for i=1:I for j=1:J if i==1 if j==1 dnnjian(i,j)=exp(1i*bb(i)*HB/2)*M(i)+alpha/miu*(1-exp(1i*bb(i)*HB/2))/HB*tanh(bb(j)*S)/bb(j)*1i*tanh(bb(i)*S); else dnnjian(i,j)=alpha/miu*(1-exp(1i*bb(i)*HB/2))/HB*tan(bb(j)*S)/bb(j)*1i*tanh(bb(i)*S); end else if j==1 dnnjian(i,j)=alpha/miu*(1-exp(1i*bb(i)*HB/2))/HB*tanh(bb(j)*S)/bb(j)*1i*tan(bb(i)*S); else if j==i dnnjian(i,j)=exp(1i*bb(i)*HB/2)*M(i)+alpha/miu*(1-exp(1i*bb(i)*HB/2))/HB*tan(bb(j)*S)/bb(j)*1i*tan(bb(i)*S); else dnnjian(i,j)=alpha/miu*(1-exp(1i*bb(i)*HB/2))/HB*tan(bb(j)*S)/bb(j)*1i*tan(bb(i)*S); end end end end end for i=1:I for j=1:J if i==1 if j==1 pnnjian(i,j)=alpha/miu*1i*tanh(bb(i)*S)/HB*(1-exp(1i*bb(i)*HB/2))*tanh(bb(j)*S)/bb(j); else pnnjian(i,j)=alpha/miu*1i*tanh(bb(i)*S)/HB*(1-exp(1i*bb(i)*HB/2))*tan(bb(j)*S)/bb(j); end else if j==1 pnnjian(i,j)=alpha/miu*1i*tan(bb(i)*S)/HB*(1-exp(1i*bb(i)*HB/2))*tanh(bb(j)*S)/bb(j); else pnnjian(i,j)=alpha/miu*1i*tan(bb(i)*S)/HB*(1-exp(1i*bb(i)*HB/2))*tan(bb(j)*S)/bb(j); end end end end for i=1:I if i==1 fnjian(i)=-alpha/miu*tanh(bb(i)*S)/bb(i); else fnjian(i)=-alpha/miu*tan(bb(i)*S)/bb(i); end end for i=1:I for j=1:J if i==1 if j==1 cmmjian(i,j)=-1i*kk(i)*N(i); else cmmjian(i,j)=0; end else if j==i cmmjian(i,j)=-1i*kk(i)*N(i); else cmmjian(i,j)=-0; end end end end for i=1:I for j=1:J if i==1 if j==1 dmnjian(i,j)=-1i*bb(j)*exp(1i*bb(j)*HB/2)/cosh(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sinh(kk(i)*S)-bb(j)*cosh(kk(i)*S)*sinh(bb(j)*S))/(kk(i)^2-bb(j)^2); else dmnjian(i,j)=-1i*bb(j)*exp(1i*bb(j)*HB/2)/cosh(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sinh(kk(i)*S)+bb(j)*cosh(kk(i)*S)*sin(bb(j)*S))/(kk(i)^2+bb(j)^2); disp(dmnjian(i,j)); end else if j==1 dmnjian(i,j)=-1i*bb(j)*exp(1i*bb(j)*HB/2)/cos(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sin(kk(i)*S)+bb(j)*cos(kk(i)*S)*sinh(bb(j)*S))/(kk(i)^2+bb(j)^2); disp(dmnjian(i,j)); else dmnjian(i,j)=-1i*bb(j)*exp(1i*bb(j)*HB/2)/cos(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sin(kk(i)*S)-bb(j)*cos(kk(i)*S)*sin(bb(j)*S))/(kk(i)^2-bb(j)^2); end end end end %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% 速度连续条件 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% PH\PB %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:I for j=1:J if i==1 if j==1 enmjia(j,i)=1/cosh(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sinh(kk(i)*S)-bb(j)*sinh(bb(j)*S)*cosh(kk(i)*S))/(kk(i)^2-bb(j)^2); else enmjia(i,j)=1/cosh(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sinh(kk(i)*S)+bb(j)*sin(bb(j)*S)*cosh(kk(i)*S))/(kk(i)^2+bb(j)^2); end else if j==1 enmjia(j,i)=1/cos(kk(i)*D)/cosh(bb(j)*S)*(kk(i)*cosh(bb(j)*S)*sin(kk(i)*S)+bb(j)*sinh(bb(j)*S)*cos(kk(i)*S))/(kk(i)^2+bb(j)^2); else enmjia(i,j)=1/cos(kk(i)*D)/cos(bb(j)*S)*(kk(i)*cos(bb(j)*S)*sin(kk(i)*S)-bb(j)*sin(bb(j)*S)*cos(kk(i)*S))/(kk(i)^2-bb(j)^2); end end end end for i=1:I for j=1:J if i==1 if j==1 pnnjia(i,j)=exp(1i*bb(i)*HB/2)*M(i)+alpha