这次我们的目标如下:
可以看到项目里面最大的挑战是控制船的倾翻角度。
设计步骤:
下面我们将按照这一设计步骤一步一步地来完成。
计算重心,计算空载曲线
把这两步合并在一起的原因是:船的重心其实就等于空载曲线下面排开水的部分的质心,所以当我们算出船的空载曲线时,船的重心也就出来了。
可以看我的这篇文章来了解浮心的算法。
我们计划的船型是一个100X30X20(cm)的船。
为了方便计算与拼装,我们不妨设船舶的每一片切面的形状都是相同的,只是大小发生变化。
于是我们设计出船舶的切片的标准形状:
你可以认为船的每一片切面都是由这个“标准切面”等比例缩放而来的
而这个等比例缩放在x轴满足下面的分布:
有了这几点前提我们的空载曲线就很好算了:
首先定义出水平面:
%这条是空载水线z=-4
d = -4;
theta = 0;
y = mesh.ys;
z = tand(theta) .* y + d;%每个水线对应的方程
为什么空载曲线是z=-4?
其实是我们利用二分的思想一次次试出来的,
我们计划船舶的重量是5800g左右,排开水的重量随排水线的高度单调增加,所以我们二分排水线高度,当排开水的体积等一重量时,这条水线就是我们的空载水线。
下面是计算浮心和重量的代码:
%这条是空载水线z=-4
d = -4;
theta = 0;
y = mesh.ys;
z = tand(theta) .* y + d;%每个水线对应的方程
water = mesh.zgrid < z;
bluemap = [1,1,1; 0,0,1];
plotMatrix(water,mesh,bluemap);
%计算水线与切面的公共面积
sub_region = hull & water;%&位与运算(都是1才得1)
purplemap = [1,1,1; 1,0,1];
plotMatrix(sub_region,mesh,purplemap);%绘制图像
mass=matrixSum(sub_region)
COB = centerOfMass2(sub_region,mesh);
V=0;
%计算浮心2D
for i=1:30 %X从0到100切30次每次长度为3
Zjudge=Judge(i*3)
bili=(10-Zjudge)/20.0%相关的缩放比例
if i>=5 && i<=25
bili=1;
end
yyy=COB(1);%标准面的重心的y坐标
zzz=COB(2);%标准面的重心的x坐标
ZMass(end+1)=Zjudge+(zzz+10)*bili
YMass(end+1)=yyy*bili
XMass(end+1)=(i*3)
MMass(end+1)=mass*bili
tot_mass=tot_mass+mass*bili;%用来算总质心的权重比
V=V+mass*3*bili;%总体积
end
ANS=centerOfMass3(XMass,YMass,ZMass,MMass,tot_mass);%计算三维浮心
M=V*1%总重量
算出来在这条水线下的排开水的体积等于4275g(还会加上电机电池与两个载重)与我们船的质量5800g相近,我们就得出了船的空载水线z=-4;
但是算质心时,我们应该再加上电机和电池的重心和载重物品的重心:
于是:
算出裸船的重心在(48,0,-5.6)处。
再加上两个电机一个电池(800g),重心在(65,0,-6)
ANS(1)=ANS(1,1)*M/(M+800)+65*800/(M+800)
ANs(2)=ANS(1,2)*M/(M+800)+0*800/(M+800)
Ans(3)=ANS(1,3)*M/(M+800)+(-6)*800/(M+800)
ANS
M=M+800;
再加上载重(800g),重心再在(50,0,6)
ANS(1)=ANS(1)*M/(M+800)+50*800/(M+800)
ANs(2)=ANS(2)*M/(M+800)+0*800/(M+800)
Ans(3)=ANS(3)*M/(M+800)+(6)*800/(M+800)
ANS
M=M+800
得出船在正常情况载重为800g时的重心:
(48.4,0,-5.6)、
船重:5800g
计算130度水线的浮心
我们只要把上面的水线的代码稍作改动,使得它的角度为130度排开水的体积是船的重量即可,经一番测试,我们得出水线为:
d = 9;
theta = 50;%取>的部分故反过来为50度
y = mesh.ys;
z = tand(theta) .* y + d;%每个水线对应的方程
时,截取的船截面为:
计算这种水线下的浮心为:
于是我们得出了130度时的浮心位置。
计算扶正力矩:
给出计算图:
会发现在130度的情况下船舶是具有扶正能力的,且在140度时不具备扶正能力。
后续
接下来就是按照这一模型进行soildwork建模,造出实物来验算了。
完整版代码如下:
计算浮心(单位:厘米):
clear all
boat.L = 100; %船的长度(厘米)
boat.W = 30; %船的宽度
boat.HB = boat.W / 2; %船的半宽
boat.D = 20; %船的深度
boat.HD = boat.D / 2; %船的半宽
max_area = boat.D * boat.W %沿X轴的切面(Y_Z平面切面)的最大面积
max_volume = max_area * boat.L %立方体的体积
density_water = 1; %水的密度 (克/立方厘米)
max_mass = max_volume * density_water %排开水的最大质量
从Y,Z轴切片(取微元)
dy = 1; % meters(微元_Y轴)
dz = 1; % meters(微元_X轴)
mesh.ys = -1.0*(-boat.HB:dy:boat.HB); % meters Y轴切片
mesh.zs = -1.0*(-boat.HD:dz:boat.HD); % meters Z轴切片
[mesh.ygrid,mesh.zgrid] = meshgrid(mesh.ys,mesh.zs);
total_area = boat.W * boat.D %平方英尺
mesh.dA = total_area / numel(mesh.ygrid) %微元面积
%XMass,YMass,ZMass%记录的是每个切片的浮心(2D)
ZMass=[];
XMass=[];
YMass=[];
MMass=[];%每个浮心点的权重(质量)
tot_mass=0;%总质量
对每个x所对应的点切片:
ZZ=[];
for i=0 : 30 %讨论y从[-10,10]
t=caculate(i);
ZZ=[ZZ t];
end
y = mesh.ys;
figure(1);%绘制不同的图
%axis xy;
hull = mesh.zgrid > ZZ%在所计算的Z的值之上的就是船的截面
redmap = [1,1,1;1,0,0];
colormap(redmap);
image(mesh.ys,mesh.zs,flipud(hull),'AlphaData',0.5);
描述水线 y=kx+d
%这条是空载水线z=-4
d = -4;
theta = 0;
y = mesh.ys;
z = tand(theta) .* y + d;%每个水线对应的方程
water = mesh.zgrid < z;
bluemap = [1,1,1; 0,0,1];
plotMatrix(water,mesh,bluemap);
%计算水线与切面的公共面积
sub_region = hull & water;%&位与运算(都是1才得1)
purplemap = [1,1,1; 1,0,1];
plotMatrix(sub_region,mesh,purplemap);%绘制图像
mass=matrixSum(sub_region)
COB = centerOfMass2(sub_region,mesh);
V=0;
%计算浮心2D
for i=1:30 %X从0到100切30次每次长度为3
Zjudge=Judge(i*3)
bili=(10-Zjudge)/20.0%相关的缩放比例
if i>=5 && i<=25
bili=1;
end
yyy=COB(1);%标准面的重心的y坐标
zzz=COB(2);%标准面的重心的x坐标
ZMass(end+1)=Zjudge+(zzz+10)*bili
YMass(end+1)=yyy*bili
XMass(end+1)=(i*3)
MMass(end+1)=mass*bili
tot_mass=tot_mass+mass*bili;%用来算总质心的权重比
V=V+mass*3*bili;%总体积
end
ANS=centerOfMass3(XMass,YMass,ZMass,MMass,tot_mass);%计算三维浮心
M=V*1%总重量
算出裸船的重心在(48,0,-5.6)处。
再加上两个电机一个电池(800g),重心在(65,0,-6)
ANS(1)=ANS(1,1)*M/(M+800)+65*800/(M+800)
ANs(2)=ANS(1,2)*M/(M+800)+0*800/(M+800)
Ans(3)=ANS(1,3)*M/(M+800)+(-6)*800/(M+800)
ANS
M=M+800;
再加上载重(800g),重心再在(50,0,6)
ANS(1)=ANS(1)*M/(M+800)+50*800/(M+800)
ANs(2)=ANS(2)*M/(M+800)+0*800/(M+800)
Ans(3)=ANS(3)*M/(M+800)+(6)*800/(M+800)
ANS
M=M+800
得出船在正常情况载重为800g时的重心:
(48.4,0,-5.6)、
船重:5800g
function z=Judge(zj)
if zj<15
z=-4.0/3*zj+10;
end
if zj>75
z=4.0/3*(zj-75)-10;
end
if zj>=15 && zj<=75
z=-10;
end
end
function ANS=centerOfMass3(XMass,YMass,ZMass,MMass,m)
x=0;
y=0;
z=0;
for i=1:30
x=x+XMass(1,i)*MMass(1,i)/m;
y=y+YMass(1,i)*MMass(1,i)/m;
z=z+ZMass(1,i)*MMass(1,i)/m;
end
ANS=[x,y,z];
end
function M = matrixSum(A)
% matrixSum: returns total of all elements in the matrix
% A: matrix
% returns: scalar
% normally sum(m) computes the sums of the columns
% selecting m(:) flattens the matrix and computes the sum of all elements
% see https://stackoverflow.com/questions/1721987/what-are-the-ways-to-sum-matrix-elements-in-matlab
M = sum(A(:));
end
function COM = centerOfMass2(masses,mesh)
% centerOfMass2: computes center of mass in 2D
% masses: matrix of masses
% mesh: structure containing ygrid and zgrid
% returns: Vector [ycom,zcom]
M = matrixSum(masses);
ycom = matrixSum(masses .* mesh.ygrid) / M;
zcom = matrixSum(masses .* mesh.zgrid) / M;
COM = [ycom,zcom];
end
function plotMatrix(A,mesh,cmap)
% plotMatrix: plots a matrix using image
% A: matrix
% mesh: srtruct containing ys and zs
% cmap: Colormap
colormap(cmap);
image(mesh.ys,mesh.zs,flipud(A),'AlphaData',0.5);
end
分段函数
function z=caculate(a)
a=a-15;
if a<=-11
z=-5.0/2.0*a-55.0/2;
end
if a>=11
z=5.0/2.0*a-55.0/2;
end
if a>(-11) && a<11
z=10.0/121*(a-11)*(a+11);
end
end