代码合集
2024年8月8日大约 7 分钟
2016-A
% 2016 CUMCM problem A - Optimal Design of Mooring System
% -------------------------------------------------------------------------
% Question 1
%
Lc = 22.05; chain = 2; M = 1200; depth = 18;
[tilt,elev,xsbed,xbuoy,~] = moor(Lc, chain, 12, 0, M, depth, 1);
[tilt,elev,xsbed,xbuoy,f] = moor(Lc, chain, 24, 0, M, depth, 1);
% -------------------------------------------------------------------------
% Question 2
%
% 2.1
Lc = 22.05; chain = 2; M = 1200; depth = 18;
[tilt,elev,xsbed,xbuoy,~] = moor(Lc, chain, 36, 0, M, depth, 1);
Mi = 1200:10:4500; tilti = []; elevi = [];
for mi = Mi
[tilt,elev,xsbed,xbuoy,f] = moor(Lc,chain,36,0,mi, depth);
tilti = [tilti, tilt]; elevi = [elevi, elev];
end
figure('name','problem 2');
plot(Mi, tilti,'r', Mi, elevi, 'b');
xi = [min(Mi), max(Mi)];
hold on;
plot(xi,[5,5],'r--',xi,[16,16],'b--');
set(gca,'ytick',[0,5,10,16,20]); grid on
legend('Tilt angle of the drum','Elevation angle of the chain at the anchor')
xlabel('Mass of the heavy ball (kg)'); ylabel('Angle (degree)')
idx = find( tilti<=5 & elevi<=16);
[tilt,elev,xbuoy,f] = moor(Lc, chain, 36, 0, min(Mi(idx)), depth, 1);
% -------------------------------------------------------------------------
% Question 3
Lc = 15.84; chain = 5; vw = 36; vs = 1.5; M = 4020; depth = 16;
[tilt,elev,xsbed,xbuoy,f] = moor(Lc,chain,vw,vs,M,depth,1);
Lc = 20.88; chain = 5; vw = 36; vs = 1.5; M = 4000; depth = 20;
[tilt,elev,xsbed,xbuoy,f] = moor(Lc,chain,vw,vs,M,depth,1);function [tilt,elev,xsbed,xbuoy,f] = moor(Lc,chain,vw,vs,M,depth,isplot)
% MOOR 2016 CUMCM Problem A - Optimal Design of Mooring System
% Reference: http://canuck.seos.uvic.ca/rkd/mooring/moordyn.php
%
% USAGE: [tiltdrum,elevanch,xsbed,xbuoy,f] = moor(Lc,chain,v,M,isplot)
%
% tilt = 滚筒沿着竖直方向倾斜的角度
% elev = 与海床的夹角
% xsbed = 倾斜在海床上的链条的长度
% xbuoy = swimming range of the buoy
% f = immersion ratio of the of the buoy
% Lc = 当前使用的锚链的长度
% ls = 锚链长度
% vw = 风速
% vs = 海水表面海水的流速
% M = 重物的质量
% isplot = is plot?
% g = 重力加速度
% rho = 海水密度
% rhoFe = 铁的密度
% cdwin = 风中阻力系数
%
if nargin==0; % nargin函数输入变量数
Lc=22.05;chain=2;vw=12;vs=0;M=1200;depth=18;isplot=1;
end
if nargin==6;
isplot=0;
end
g = 9.81; % g:重力加速度 [m/s^2]
rho = 1.025e3; % rho:海水密度 [kg/m^3]
rhoFe = 7.9e3; % rhoFe:铁的密度 [kg/m^3]
cdwin = 0.625; % cdwin:风中阻力系数
[lc, mc, dc] = chainpara(chain);
nc = round(Lc/lc); % 链条单元数
% Lc是当前使用的锚链的总长度。
% lc是单个链条单元的长度，这个值是通过调用chainpara函数并传入链条类型chain得到的。
% [浮标, 四个钢管, 圆柱形滚筒, 一系列焊接链条]
h = [ 2, 1*ones(1,4), 1, lc*ones(1,nc)]; % m
% 2：浮标的高
% ones(1,4)会创建一个1行4列的矩阵，每个元素都是1
d = [ 2, 5e-2*ones(1,4), 0.3, dc*ones(1,nc)]; % m
% 2：浮标的底面直径
m = [1000, 10*ones(1,4), 100, mc*ones(1,nc)]; % kg
% 1000：浮标的质量
phi = zeros(1,length(h)); % phi : 用来存储各个组件的倾斜角度或方向角
% 浮标的浮力
Fb = pi*(d/2).^2.*h*rho*g - m*g; % Turn masses/buoyancies into forces将质量/浮力转化为力
% 二分法求解浮标吃水比例
fmin = 0; fmax = 1;
while fmax-fmin>1e-10
% f : 浮标的浸没比例
f = (fmax+fmin)/2;
% Fb : 浮标的浮力(1)
Fb(1) = f*pi*(d(1)/2).^2.*h(1) * rho * g - m(1)*g;
% 风载
Fw = cdwin * (1-f)*h(1).*d(1) * vw.^2; % wind load
% 水载
Fs = waterload(vs, h, d, phi, depth, f);% water load
% 每个杆的倾斜角
phi = solvequileq(Fb, Fw, Fs, M, f);
% 组件在水平方向上的投影长度
x = h.*sin(phi);
% 组件在垂直方向上的投影长度
z = h.*cos(phi); % projected length
% 浮标底部到海平面的垂直距离
zw = sum(z(2:end)) + h(1)*f; % waterline
if zw>depth;
fmax = f;
else;
fmin = f;
end
end
% ①fliplr(x)和fliplr(z)用于将向量x和z左右翻转；
% 原来向右的分量现在向左
% 原来向下的分量现在向上
% ②[0 fliplr(x)]和[0 fliplr(z)]是将翻转后的向量前加上一个0，形成新的向量。
% 这个0代表锚泊系统的起始点，通常是在海床上。
% ③cumsum是MATLAB函数，用于计算输入向量的累积和，
% 对于每个元素，cumsum返回该元素及其之前所有元素的总和。
% 从锚开始，到浮标的水平距离
x = cumsum([0 fliplr(x)]);
% 从锚开始，到浮标的垂直距离
z = cumsum([0 fliplr(z)]);
tilt = phi(6)*180/pi; % tilt angle of the drum滚筒倾斜角度
elev = 90- phi(end)*180/pi; % elevation angle of chain at anchor锚处链条的仰角
xsbed = max(x(z<1e-10)); % 海床上链条的长度
xbuoy = x(end-1); % swimming range of the buoy浮标的游动范围
if isplot;
plotmoor(x,z,h,phi,Lc,chain,vw,vs,M,depth,tilt,elev,xbuoy,f);
end
% ------------------------------------------------------------------------
function Fs = waterload(vs, h, d, phi, depth, f)
cd = 374; % 海水阻力系数
z = h.*cos(phi);
zi = fliplr(cumsum(fliplr(z))) - h.*cos(phi)/2;
zi(1) = depth - f*h(1)/2;
vsi = vs./sqrt(depth)*sqrt(zi); % 水在不同深度的流速
Fs = cd * h.*d.*cos(phi) .* vsi.^2; % water load
Fs(1) = Fs(1)*f;
% ------------------------------------------------------------------------
% 解决锚泊系统的平衡方程，以确定系统中每个组件的倾斜角度phi
function phi = solvequileq(Fb, Fw, Fs, M, f)
g = 9.81; % acceleration of gravity [m/s^2]
rho = 1.025e3; % density of seawater [kg/m^3]
rhoFe = 7.9e3; % density of ferrum [kg/m^3]
N = length(Fb);
[theta, phi, Ft] = deal(zeros(1,N));
for i = 1:N-1
% theta：拉力与竖直方向的夹角
%
% 链在水平方向上
fx = Ft(i)*sin(theta(i)) + Fs(i);
% 如果是浮标，
if i==1; fx = fx + Fw; end
% 链在竖直方向上
fz = Fb(i) + Ft(i)*cos(theta(i));
% 加上重物球 = 重力 - 浮力
if i==6; fz = fz - M*g + rho*(M/rhoFe)*g; end
% 竖直方向与水平方向上的合力
Ft(i+1) = sqrt(fx^2+fz^2);
% 得到theta角
theta(i+1) = acos(fz/Ft(i+1)); % the tilt of Ft
if theta(i+1)>pi/2; theta(i+1) = pi/2; end
end
phi = atan2( Ft.*sin(theta)+Fs/2, Ft.*cos(theta)+Fb/2);
phi(phi>pi/2) = pi/2; % the tilt of each element
phi(1) = atan2( Fs(1)*f/2+Fw(1)*(f+(1-f)/2), Fb(1)*f/2 );
% ------------------------------------------------------------------------
function [lc, mc, dc] = chainpara(typeid)
% The types and parameters of the welded chains
% 锚链的类型和参数
rhoFe = 7.9e3; % density of ferrum [kg/m^3];铁的密度
rho = [3.2 7.0 12.5 19.5 28.12]; % 海水密度
% length(lc) and mass(mc) of one unit of the welded chain.
lc = [ 78 105 120 150 180]*1e-3; % m 锚链的长度
mc = rho.*lc; % kg
lc = lc(typeid);
mc = mc(typeid);
% 根据typeid的值选择对应的链条长度
% one unit of the welded chain was assumed as a cylinder.
% 锚链的一个单元被假定为圆柱体。
dc = 2*sqrt(rho(typeid)/rhoFe/pi); % diameter of the welded chain焊接链条的直径
% -------------------------------------------------------------------------
function plotmoor(x,z,h,phi,Lc,chain,vw,vs,M,depth,tilt,elev,xbuoy,f)
% 创建一个新的图形窗口，并计算向量x的长度，将其存储在变量N中，用于后续的循环和索引
figure; N = length(x);
% 绘制浮标
x1 = x(N) -cos(phi(1)); x2 = x(N) +cos(phi(1)); % x(N)表示取向量x的第N个元素
x3 = x(N-1)+cos(phi(1)); x4 = x(N-1)-cos(phi(1));
z1 = z(N) +sin(phi(1)); z2 = z(N) -sin(phi(1));
z3 = z(N-1)-sin(phi(1)); z4 = z(N-1)+sin(phi(1));
fill([x1 x2 x3 x4], [z1 z2 z3 z4], 'r');
hold on
text(x(N-1)+1.3, z(N-1), sprintf('%2.1f m',f*h(1)));
% 绘制铁管
plot(x(N-5:N-1),z(N-5:N-1),'o-','linewidth',1, 'markersize',2);
% 绘制铁管，并添加倾斜角度标注
for i = 5:-1:2 % 从第5根钢管（索引为5，即最后一根钢管）到第2根钢管
xi = x(N-i+1);
zi = z(N-i+1);
angle_deg = 90 - phi(N-i+1)*180/pi; % 转换为度
plot(xi, zi, 'o-', 'linewidth', 1, 'markersize', 2);
text(xi+0.1, zi, sprintf('%3.1f^\\circ', angle_deg), 'HorizontalAlignment', 'left');
end
% 绘制铁桶
x1 = x(N-5)-0.15*cos(phi(6)); x2 = x(N-5)+0.15*cos(phi(6));
x3 = x(N-6)+0.15*cos(phi(6)); x4 = x(N-6)-0.15*cos(phi(6));
z1 = z(N-5)+0.15*sin(phi(6)); z2 = z(N-5)-0.15*sin(phi(6));
z3 = z(N-6)-0.15*sin(phi(6)); z4 = z(N-6)+0.15*sin(phi(6));
fill([x1 x2 x3 x4], [z1 z2 z3 z4], 'r');
hold on
text(x(N-6)-2,z(N-6)+0.5,sprintf('%3.1f^\\circ',tilt));
% 绘制锚链
plot(x(1:N-6),z(1:N-6),'-m')
text(0.5,1.5,sprintf('%3.1f^\\circ',elev));
types = {'I','II','III','IV','V'};
text(x(ceil(N/2))+1, z(ceil(N/2)), sprintf('%s:%5.1f m',types{chain}, Lc));
% 绘制重物球
plot([x(N-6) x(N-6)],[z(N-6) z(N-6)-4],'-r');
plot(x(N-6),z(N-6)-4,'.r','markersize',50)
text(x(N-6)-1,z(N-6)-5,sprintf('%d kg',M));
axismax = max(depth+2,xbuoy+4);
% 绘制水平面
plot([0,axismax],[depth,depth],'--');
% 风速
quiver(x(N)-8,depth+1,2,0,'MaxHeadSize',0.5,'color','r');
text(x(N)-5.5,depth+1,sprintf('%3.1f m/s', vw),'color','r');
% 水速
quiver(x(N)-8,depth-1,2,0,'MaxHeadSize',0.5,'color','b');
text(x(N)-5.5,depth-1,sprintf('%3.1f m/s', vs),'color','b');
box on; axis image;
set(gca,'xtick',[max(x(z<1e-10)),round(xbuoy*10)/10],'xgrid','on')
axis([0,axismax,0,axismax]);
xlabel('x'); ylabel('z')