船舶强度与结构设计大作业（一）船海1301 禹宗昕U201312263一．数据与函数准备1.主要数据：计算船长L=200m水密度ρ=1t/m³重力加速度g=9.8m/s²2.原始资料（见附录）（1）全船重量分布（2）全船邦戎曲线数据3.静水平衡参数总重量W=122248 *9.8=120030.4 kN水线面面积A=4800m平均吃水dm=3.9m纵稳心半径R=220m漂心纵向坐标xf=4.3m4.基本函数简介（1）function [ Area ] = SAREA( DD )% 作用：已知各站吃水求各站对应的水下面积,核心方法：三次样条曲线差值法% D：已知水线% br：原始邦戎数据，21*21矩阵，(因数据过大在下面已省略)br(i,:)代表第i行每站在各吃水下的面积，% DD：参数，一维矩阵，各站实际吃水% Area：输出，一维矩阵，各站在DD吃水情况下对应的横截面积D=[0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12;]br=[…..];for i=1:length(DD);Area(i)=spline(D,br(i,:),DD(i));end(2 ) function [ D ] = draft( df,da,n )% 作用：已知首尾吃水求各站吃水% df：参数，double，首吃水% da：参数，double，尾吃水% n:站数% D：输出，一维矩阵，n站对应吃水c=(df-da)/(n-1);for i=1:1:nD(i)=da+(i-1)*c;end;end（3）function [ V] = displacement( area,zhanju )% 已知各站面积求排水，核心：辛普森面积法% area：参数，一维矩阵，各站排水面积% zhanju：参数，double，站距% V：输出，double，排水体积S=0;for i=1:2:length(area)-2;S=S+(area(i)+4*area(i+1)+area(i+2))*2/6;V=S*zhanju;end;（4）function [ xb1] = flocenter( area,zhanju )% 作用：求浮心位置；核心：displacement函数拓展，141辛普森面积法% area：参数，一维矩阵，各站排水面积% zhanju：参数，一维矩阵，站距% xb1：输出，double，浮心到尾垂线的距离S=0;for i=1:2:length(area)-2;S=S+(area(i)+4*area(i+1)+area(i+2))*2/6;V=S*zhanju;end;%mb：各站体积对尾垂线的矩%下述两个for循环，是奇数站和偶数站对应的辛普森系数，第一站和最后一站的系数为0.5，其余为1，2交叉mb=area(1)*0.5*0+area(length(area))*0.5*(length(area)-1);for i=2:2:length(area)-1;mb=mb+area(i)*2*(i-1);end;for i=3:2:length(area)-2;mb=mb+area(i)*1*(i-1);end;xb1= mb*zhanju/V/3*2*zhanju;end二． 浮态计算1.总重与重心计算公式：重心距离尾垂线xg1 = ∑M i 201∑W i 201= E23/C24*10 = 93.2136 m 2.浮态计算（1）函数说明 function [ d] = floatingstate(dm,L,A,R,xf,W,xg,zhanju)%作用：通过迭代计算浮态，% dm：正浮时平均吃水，% A：正浮时水线面面积% R：纵稳心半径% xf：漂心% W：总重% xg：浮心距尾垂线距离% zhanju：站距% d：数据类型：类，% d.xg重心距舯m% d.xb浮心距舯m% d.weight总重t% d.displacement排水t% d.index1=index1=abs((W-V)/W);判断指标1，当误差小于0.5%时停止迭代% d.index2=index2=abs((xg-xb)/L);判断指标2，当误差小于0.1%时停止迭代% d.k 迭代次数k=0;index1=1;index2=1;%赋初始值，第一次取平均吃水，进行第一次排水与浮心计算d.df=dm; d.da=dm;D=draft(d.df,d.da,(L/zhanju+1));area=SAREA(D);V=displacement(area,zhanju);xb=flocenter(area,zhanju);%while循环迭代，当误差小于0.5%和0.1%时停止迭代，迭代步长：0.3 while (index1>0.005 )||(index2>0.001)&&(k<20)d.df=d.df+((W-V)/(p*A)+(L/2-xf)*(xg-xb)/R)*0.3;d.da=d.da+((W-V)/(p*A)-(L/2+xf)*(xg-xb)/R)*0.3;D=draft(d.df,d.da,(L/zhanju+1));area=SAREA(D);V=displacement(area,zhanju);xb=flocenter(area,zhanju);d.xg=xg-100;d.xb=xb-100;d.weight=W;d.displacement=V;index1=abs((W-V)/W);index2=abs((xg-xb)/L);d.index1=index1;d.index2=index2;k=k+1;d.k=k;end;end（2）运行结果d = floatingstate (3.9,200,4800,220,-4.3,12248,93.214,10)d =df: 3.0666da: 4.6647xg: -6.7860xb: -6.6036weight: 12248displacement: 1.2240e+004index1: 6.5099e-004index2: 9.1211e-004k: 3浮态分析：da>df 尾倾精度分析：index1=abs((W-V)/W)= 0.065% < 0.5%Index2=abs((xg-xb)/L)=0.091% < 0.1%迭代次数：k=3三．剪力与弯矩计算1.初次计算（1）.载荷分布，剪力，弯矩的计算与储存% D：一维矩阵，已修正浮态下各站吃水，% area：一维矩阵，已修正浮态下各站面积% earea：一维矩阵，利用area和三次样条曲线差值，求0.5,1.5..站面积% eachv：各站的体积，*10为站距% qx：载荷分布，重力-浮力，qx1为集中力到分布力的换算% N11：剪力函数，M11：弯矩函数% A，B11，C11：分别为qx1，N11，M11对应的多项式系数矩阵D=draft(d.df,d.da,21)；area=SAREA(D)；earea=spline([0:1:20],area,[0.5:1:19.5])；eachv=10*earea；qx=eweight-eachv；qx1=qx/10*9.8;x=[5:10:195];A11=polyfit(x,qx1,4);syms x q1;q1=poly2sym(A,x);syms N11 M11;N11=int(q1);M11=int(N11);B11=sym2poly(N11); C11=sym2poly(M11);x=[0:10:200];ezplot(M11,x) ;hold on;ezplot(100*N11,x);（2）图形校核原理：1.多项式A的根为N的极值点，B的根为M的极值点，a=roots(A)2.将极值点带入函数求极值，如polyval(B,a(2))3.求第20站的剪力，弯矩值，如polyval(B,20)4.计算k1，k2指标，k1=Nmax/N（20）；k2=Mmax/M（20）；5.k1&&k2<0.05时可进行第二步修正a=roots(A11)a =20.875014.2251-7.78004.0257>> polyval(B11,a(2))ans =-1.6803e+004>> polyval(B11,a(4))ans =1.9640e+004>> polyval(B11,20)ans =6.8987e+001>> k1=ans/ 1.9640e+004k1 =0.0035>> b=roots(B)b =-11.665421.661619.97159.2145>> k2=polyval(C,20)/polyval(C,b(5))k2 =0.0470（3）结果分析k1=Nmax/N（20）=0.0035 <0.05k2=Mmax/M（20）=0.0470 <0.05由图可见，20站处N和M均不为0，即此积分方法存在误差，但均在误差允许范围内，因此可进行二次修正。

2.修正原理：（1.）N1=N-x*N（20）/20 ；M1= M-x*M（20）/20（2）N1对应B11为五次多项式，因而一次项为第五项，M1对应C11为六次多项式，因而一次项为第六项，B22=B11;C22=C11;B22(5)= B22(5)- polyval(B11,200)/200;C22(6)=C22(6)- polyval(C11,200)/200; syms x N1 M1;N1=poly2sym(B22,x);M1=poly2sym(C22,x);x=[0:10:200];ezplot(M1,x) ;hold on;ezplot(100*N1,x);x=[0:1:20];ezplot(M1,x) ;hold on;ezplot(100*N1,x);分析：利用matlab自带的图像分析功能，易得修正后，M1max=1.164*106 KN*mN1max=2*104 KN四．心得体会船舶强度与结构设计已经开课三周，给我最深的感触是，它突破了原先我们结构设计的简单认识，以为设计出来就万事大吉。

课程里，老师无数次强调反复修正的重要性，然而这门课的大作业也是我们上大学以来第一次面临计算的迭代问题。

因为在计算工具上，我选择了matlab。

之前在静力学课设上，同学们普遍使用AutoCAD绘制邦戎曲线，然后描点，量尺寸，求每一站的面积，然而如今面临迭代，手工求法显然不再使用。