系泊系统的设计 摘要 本文对系泊系统的设计问题进行了分析，给出合理的假设，建立优化模型，巧妙地解决了题目中所提出的问题。
针对问题一，首先采用集中质量的多边形近似法对单点系泊系统进行静力学分析，结合单点浮标系统特性，建立单点浮标的静力学模型，并对其算法进行改进，使算法能够迭代修正浮标受到的浮力。其次通过适当的假设列出平衡方程并求解，得出锚链各节点处张力的递推公式，利用MATLAB软件迭代验证，最后得到了较为准确的结果。
针对问题二，基于问题一建立的模型，得出在题设条件下，浮标系统已不处于最优工作状态，须通过改变重物球来对系统进行调节。计算出临界条件下重物球的质量，利用MATLAB拟合得到的重物球重力与钢桶倾斜角度之间的关系曲线，得出对重物球进行调整的范围。
针对问题三，首先求得极端环境条件下钢桶倾角仍满足约束条件时候的重物球质量，然后通过合理的假设，在问题一建立的模型基础上，改变算法的迭代约束条件，从而得出不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。
关键词：系泊系统；
集中质量的多边形近似法；
MATLAB；
迭代 一、问题的重述 1.1问题的背景 随着各国不断加大对海洋事业的投入以及不断深入对海底观测领域的探索，各海洋研究机构和海洋管理部门都相继建立了符合自身业务需求的海洋观测系统，其中浮标系统、系泊系统和水声通讯系统组成了近浅海观测网的传输节点。而影响其系泊系统工作效果的因素很多，例如水流力、海风和水深等。系泊系统的设计问题就是根据这些影响因素确定锚链的型号、长度和重物球的质量，使得系泊系统处于最佳工作状态。从国家海洋资源战略角度来讲，研究各因素对系泊系统的不同影响显得尤为重要。
1.2问题的提出 问题一：将该型传输节点布放在水深18m、海床平坦、海水密度为1.025×103kg/m3的静止海域时，选用II型电焊锚链22.05m，质量为1200kg的重物球，分别计算海面风速为12m/s和24m/s时钢桶和各节钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。
问题二：在问题一的假设下，计算海面风速为36m/s时钢桶和各节钢管的倾斜角度、锚链形状和浮标的游动区域。要求通过调节重物球的质量，使得钢桶的倾斜角度不超过5°，锚链在锚点与海床的夹角不超过16°。
问题三：当布放海域的实测水深介于16m~20m之间，布放点的海水速度最大可达1.5m/s、风速最大可达36m/s时，在综合考虑风力、水流力和水深等因素的情况下对系泊系统进行设计，并分析不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。
二、问题分析 2.1问题一的分析 通过集中质量的多边形近似法[1]来对单点系泊系统进行静力学分析，得出钢桶和锚链之间的作用力，这是问题的关键点。分析该作用力后，结合由静力学平衡方程得出的锚链各节点间张力大小及其方向的递推公式，利用MATLAB软件进行数据计算和曲线拟合，便可以求得锚链的具体形状。求解钢桶和锚链之间的作用力时需要对钢桶、重物球、钢管、浮标组成的整体进行受力分析，通过MATLAB软件迭代计算可以求得较为精确的钢桶和锚链之间的作用力、钢桶的倾斜角度及浮标受到的浮力，由浮力进一步得到浮标的吃水深度。同时，由于各个钢管倾斜角度都非常小，对钢管之间的相互作用力进行近似计算可以得出各个钢管的倾斜角度以及钢桶距离海底的垂直高度，这为后续判定锚链末端是否拖地提供了约束条件。经分析，浮标的游动区域应为一个圆面。当锚链形状和钢管倾斜角度已知时，计算出它们在海底的投影长度即为该圆面的半径大小，游动区域亦求得。
2.2问题二的分析 首先采用问题一的静力学分析方法，得到当前情况下各节钢管的倾斜角度、锚链形状和浮标的游动区域；
然后比较该条件下钢桶的倾斜角度和锚链末端与锚链接处的切线方向与海床的夹角是否符合要求；
最后反推计算出在临界条件下重物球的质量，结合拟合出的Gq-曲线，指导系泊系统的正常运行。
2.3问题三的分析 此题与前两问题相比增加了近海水流力等因素的影响，使问题更加复杂化。本题延续了前两问题的部分求解思想，先求出极端情况下重物球的质量作为中间值，再结合问题一、二的结论，以此展开对多变量影响结果的计算，便可分析水深、水流速度、风速等对钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域的影响。基于问题一的模型，可以求得不同锚链型号条件下所需锚链的长度，完成系泊系统的设计。
三、模型假设 1. 在不同海水深度中，海中水流力和水流速度恒定，与深度无关；
2. 波浪对系统载荷无影响；
3. 浮标所受的重力、浮力以及海风对它产生作用力的作用点近似为同一点；
4. 浮标受到海风和水流力作用时不发生倾斜；
5. 假定海流为平面流，在垂直方向上无分量，将分析简化为二维问题；
6. 假设近海风负载的方向与海平面水平；
7. 由于电焊锚链的特殊结构，海水对锚链的作用力忽略不计；
8. 假定重物球所受到的浮力忽略不计。
四、符号说明 表4.1 相关符号说明 符号 说明 符号 说明 H 布点水深 L6 钢桶与锚链顶端的作用力 ρ 海水密度 Fw 近海水流力 Gi 锚链第i个节点的重力 h0 浮标吃水深度 mi 锚链单位长度的质量 φi Li与竖直方向上的夹角 g 重力加速度 Ft 钢桶浮力 Ti 锚链第i个节点处的张力 Gt 钢桶和重物球总重力 Fi 锚链第i个节点处的浮力 xi 锚链第i个节点处的x坐标 αi Ti与水平方向的夹角 yi 锚链第i个节点处的y坐标 F 浮标所受浮力 S风 物体在风向法平面的投影面积 G 浮标重力 S水 物体在水流速度法平面的投影面积 Ff 近海风荷载 Gq 重物球重力 L1 第一根钢管与浮标的作用力 mq 重物球质量 Li 第i-1根钢管与第i根钢管的作用力(i=2,3,4) β 浮标与海平面的倾角 L5 第四根钢管与钢桶的作用力 V风/水 风/水速 五、模型的建立与求解 5.1模型建立前的准备 5.1.1锚链某节点受力分析 分析锚链线从顶端处开始第i个节点的受力情况，有重力Gi，浮力Fi，锚链线两端张力Ti和Ti-1。如图5.1所示：
图5.1锚链线某节点处受力分析 根据锚链单位长度的质量，得单位节点的重力：
（1） 由于电焊锚链的特殊结构，计算时可简化认为：
（2） 锚链线处于平衡时，可列静力平衡方程：
（3） （4） 其中，——与水平方向的夹角；
——与水平方向的夹角。
求解上述静力平衡方程，可以得出锚链每个节点处下端的张力以及水平夹角；
（5） （6） 综上可知，只要已知第i个节点处前一个节点的张力，就可以求得该节点处的张力，即。再通过MATLAB软件编程，递推出锚链线各个节点处的张力。
本文选定锚链线最高节点处张力L6为突破点进行计算，L6的大小和方向可通过对钢桶及以上部分进行局部受力分析求得。
5.1.2浮标受力分析 分析浮标倾角为零的情况下的受力情况，有浮力F，重力G，近海风荷载Ff，近海水流力Fw（问题一、二中可忽略不计），以及第一根钢管对浮标的作用力L1。如图5.2所示：
图5.2浮标不倾斜时受力分析 分析得如下平衡方程：
（7） （8） （9） （10） 其中，——物体在水流速度法平面的投影面积（m2）；
——水速(m/s)。
——物体在风向法平面的投影面积（m2）；
——风速(m/s)。
5.1.3钢管受力分析 分析四根钢管受力情况。以图5.3中a管为例进行分析，有钢管自身的重力G0，受到海水的浮力F0，浮标对钢管施加的作用力，第二根钢管对第一根钢管的作用力。如图5.3所示：
图5.3钢管受力分析 分析得如下平衡方程：
（11） （12） （13） （14） （15） （16） （17） （18） 5.1.4钢桶受力分析 如图5.4所示，为钢桶受力情况分析，有钢桶受到的浮力Ft，钢桶和重物球的总重力Gt，第四根钢管对钢桶的作用力,锚链最上端的一节对钢桶的作用力。
图5.4 钢桶受力分析 分析得如下平衡方程：
（19） （20） 5.2问题一的模型建立与求解 采用基于集中质量的多边形近似法来研究单点系泊系统的静力学问题。在处理分析系泊系统受力情况时，将其负载的各仪器部件单独作为一个节点进行质量集中处理；
在将锚链线进行分段处理时，假设每段锚链线的质量都集中到该段锚链线的几何中心。因忽略锚链在水中的受力而产生的误差可通过优化模型算法来减小。
设锚链顶端对钢桶的作用力初始值为零，以此为条件对钢桶、钢管、浮标进行受力分析，可以求得浮力的初始估值；
联立平衡方程（7）~（20）可得：
（21） 将式（21）带入式（5）和式（6）进行递推计算，如果得到的锚链线方程不与海床接触，则说明浮力初始估值过小，此时应该增大浮力的取值，编写算法程序对浮力初值进行修正，反复迭代计算，直到锚链线方程有一点和海底相切，即达到边界条件。此时反演到实际情况，锚链线上的切点到末端这一段距离的链锁已经拖到海床上。迭代算法流程如下图5.5所示。用MATLAB软件编程实现上述过程，并且可直接绘制出锚链的形状图。
图5.5 锚链迭代算法流程图 模型求解结果如下表5.1，图5.6，图5.7所示。
表5.1 问题一的模型求解结果 问题一 风速为12m/s 风速为24m/s 钢桶倾斜角度 0.997° 3.822° 第1根钢管倾斜角度 0.974° 3.736° 第2根钢管倾斜角度 0.979° 3.753° 第3根钢管倾斜角度 0.985° 3.775° 第4根钢管倾斜角度 0.991° 3.799° 浮标吃水深度/m 0.7365 0.7482 锚链末端与锚连接处切向方向与海床的夹角 0° 0° 浮标游动区域 圆心与锚在同一垂线位置，且半径为13.7m的圆面 圆心与锚在同一垂线位置，且半径为18.1m的圆面 图5.6 风速为12m/s时的锚链姿态示意图 图5.7 风速为24m/s时的锚链姿态示意图 当风速为12m/s时，从顶端计数，第143节（共计210节）链环往后全部平铺在海床上；
钢桶倾斜角度为0.997°；
各个钢管自上而下倾斜角度分别为0.974°、0.979°、0.985°、0.991°；
浮标吃水深度为0.7365m；
浮标游动区域为一圆心与锚在同一垂线位置且半径为13.7m的圆面。
当风速为24m/s时，从顶端计数，锚链第204节（共计210节）链环往后全部平铺在海床上；
钢桶倾斜角度为3.822°；
各个钢管自上而下倾斜角度分别为3.736°、3.753°、3.775°、3.799°；
浮标吃水深度为0.7482m；
浮标游动区域为一圆心与锚在同一垂线位置且半径为18.1m的圆面。
5.3问题二的模型建立与求解 通过问题一构建的数学模型，将相关数据带入，可求得如下表5.2所示的结果：
表5.2 问题二的模型求解结果 问题二 风速为36m/s 钢桶倾斜角度 8.027° 第1根钢管倾斜角度 7.855° 第2根钢管倾斜角度 7.895° 第3根钢管倾斜角度 7.941° 第4根钢管倾斜角度 7.981° 浮标吃水深度/m 0.7688 锚链末端与锚连接处切向方向与海床的夹角 16.993° 浮标游动区域 圆心与锚在同一垂线位置，且半径为18.8m的圆面 图5.8 风速为36m/s时的锚链姿态示意图 由表5.2及图5.8可知，在问题一的假设下，当海面风速达到36m/s时，各个钢管自上而下倾斜角度分别为7.855°、7.895°、7.941°、7.981°；
吃水深度为0.7688m；
浮标游动区域为一圆心与锚在同一垂线位置且半径为18.8m的圆面。其中，钢桶倾斜角度为8.027°，锚链末端与锚的链接处的切线方向与海床的夹角为16.993°，均超出了系统有效工作范畴，故此时应该改变重物球的质量来调节系统，使得锚不会被拖行，同时水声通讯系统的工作效果最佳。
结合静力学分析，在问题一编写的算法基础上，通过改进算法可求得重物球重力Gq和钢桶倾斜角度的函数关系。
如图5.9为重物球重力Gq和钢桶倾斜角度的函数关系图像，拟合得出曲线方程近似为。当纵坐标°时，横坐标N，即kg。说明在题设条件下，只有当附加在钢桶上的重物球质量达到1975.5kg时，钢桶的倾斜角度才能被控制在5°以内。同时，将这一临界重物球挂在题设中的钢桶上，改变算法程序中重物球的重力参数，通过MATLAB迭代拟合可得出新的锚链姿态示意图，如图5.10所示。
图5.9 Gq-关系图图5.10临界状态时锚链的姿态示意图 综上，在问题一的假设基础上，海面风速为36m/s时，系泊系统不在最优工作状态，只有钢桶下连接的重物球质量超过1975.5kg，钢桶的倾斜角度才可能在要求范围之内，由图5.10可知，此时锚链末端与锚链接处的切线方向与海床的夹角也满足最佳工作条件。
5.4问题三的模型建立与求解 此题与前两问题相比增加了近海水流力等因素的影响，使问题更加复杂化。此时式（8）中的水流力不可忽略，同时对问题一中的算法进行优化，使其适用于问题三的求解。
考虑极端情况，取海水速度为最大值1.5m/s，风速为最大值36m/s时，通过MATLAB编程，求解出使得钢桶倾角保持在5°时候的重物球的质量，具体算法流程图见图5.11所示。
图5.11重物球迭代算法流程图 分析得，在极端环境情况下，上述求得的重物球依旧可以维持系泊系统较为正常的工作，那么在其余环境条件下，若悬挂同样质量的重物球，钢桶倾角均符合约束条件。所以，在研究多变量对钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域的影响时，可以选定该重物球作为中间定值进行其余变量的影响分析。
表5.3 极端条件下不同锚链型号的系统特性 问题三 型号Ⅰ 型号Ⅱ 型号Ⅲ 型号Ⅳ 型号Ⅴ 钢桶倾斜角度 4.733° 4.733° 4.733° 4.733° 4.733° 第1根钢管倾斜角度 3.299° 3.299° 3.299° 3.299° 3.299° 第2根钢管倾斜角度 4.703° 4.703° 4.703° 4.703° 4.703° 第3根钢管倾斜角度 4.711° 4.711° 4.711° 4.711° 4.711° 第4根钢管倾斜角度 4.722° 4.722° 4.722° 4.722° 4.722° 浮标吃水深度/m 1.4 1.4 1.4 1.4 1.4 锚链长度/m 11.856 11.865 11.88 12.00 12.06 浮标游动区域半径/m 1.884 1.907 1.930 1.971 2.015 重物球重量/N 33797 33797 33797 33797 33797 图5.12 不同型号锚链的姿态示意图 将问题三转化为问题一以后，用同样的方法进行算法迭代计算。表5.3列出了极端条件下不同锚链型号的系统特性，图5.12展示了极端条件下不同型号锚链的姿态示意图，在进行系泊系统的设计时可以参考上述参数。从表5.3中可以看出，当重物球重量超过33797N时，无论该系统处于多么恶劣的环境中，钢桶倾角都在约束条件范围内，此时系泊系统均可以正常工作；
当选定某一重物球不变时，改变锚链的型号，钢桶和各钢管的倾角以及浮标吃水深度都不改变，只有锚链长度和浮标游动区域发生变化。
六、灵敏度与稳定性分析 6.1灵敏度分析 表6.1 钢桶倾角与风速的变化情况 钢桶倾角/(度) 0.997 3.822 8.027 风速/(m/s) 12 24 36 图6.1Gq-关系图 风速和水流速度是极易变化的参数，重物球质量也是影响钢桶倾斜角度的变量，应该考察这几个参数的灵敏度。从表6.1及图6.1可知，当风速和水流速度有微小变化时，即对钢桶倾角有明显影响。所以风速、水流速度及重物球质量关于钢桶倾斜角度的灵敏度较大。所以要密切关注风速、水流速度的变化，一旦超越临界值，需要立即调节重物球的质量，以保证系泊系统的正常运行。
6.2稳定性分析 在问题三模型的结论中，得到极端条件下不同锚链型号的系统特征，锚链的型号与钢桶倾斜角度、各钢管倾斜角度、浮标吃水深度、重物球重量无关；
不同型号锚链使用时，锚链长度变化微小，浮标游动区域半径变化微小。由表5.3可知，锚链的型号对系泊系统稳定性和参数影响不大。
七、模型的评价与推广 7.1模型的优点 （1）在对物体进行受力分析时，采用了基于集中质量的多边形近似法来进行研究，简化了模型的复杂程度。
（2）绘制锚链的姿态曲线时，加入浮力修正算法对模型进行改进，经初始值估算，反复迭代求得最佳锚链姿态。
7.2模型的缺点 （1）由于锚链姿态事先无法确定，导致循环步长大小无法准确确定。在利用MATLAB编程解决问题时，得到的结果存在一定的误差。
（2）在实际情况中，浮标受到近海风荷载和近海水流力作用时，会发生倾斜。虽然影响非常小，但模型中忽略此情况，也会引起一定误差。
7.3模型的推广及改进 首先，本题所建立的模型与实际情况较为符合，具有一定的参考价值。同时可以针对模型进行更深层次的分析，如考虑到浮标受近海风荷载和近海水流力作用发生倾斜时，假设倾角为β，其受力情况如图7.1所示。
图7.1浮标倾斜时受力情况 近海风荷载计算公式（10）：
（22） 将这个新自变量添加到算法程序中，使得本题中的算法模型更加符合实际情况，结果更精准。
同时在具体的实施中可以根据不同的实际背景对模型进行修改，在参数的设定上也可以根据实际情况而变化，如按照不同的要求，根据具体的环境情况对力学模型进行完善，更为方便、准确地解决问题。
八、参考文献 [1]王磊.单点系泊系统的动力学研究[D].中国海洋大学硕士论文,2012. [2]王文波,数学模型及其基础知识详解[M].武汉:武汉大学出版社,2006. [3]张志涌,精通MATLAB R2011a[M].北京:北京航空航天大学出版社,2011. 附录：
1.问题一的MATLAB算法程序 %暂时忽略锚链对桶的拉力，估算浮力以及钢桶、钢管角度 Gq=1200*9.8;Gt=980;Gg=40*9.8;Gf=9800;Ffg=789;v=12;Hw=18; Ff=Gq+Gt+Gg+Gf-Ffg; h=Ff*4/(1.025*1000*9.8*pi*4); Fw=0.625*(2-h)*2*v^2; datafai=[]; fai=atan(Fw/(Ff-Gf)); datafai=[datafai;fai]; T=Fw/sin(fai); for j=1:1:4 %对每一节钢管进行计算 fai=atan(Fw/(T*cos(fai)+19.72-98)); T=Fw/sin(fai); datafai=[datafai;fai]; end fai=atan(Fw/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=Fw/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; %按照所估浮力值计算锚链的形状和受力情况 datax=[];datay=[];dataT=[];datatheta=[]; G=7.203; theta=pi/2-fai; x=0; y=0; for i=1:1:210 %计算每一节锚链的受力和角度 datax=[datax;x]; datay=[datay;y]; dataT=[dataT;T]; datatheta=[datatheta;theta]; if y>H-Hw %判断是否到底 x=x-0.105*cos(theta); y=y-0.105*sin(theta); else x=x-0.105; y=H-Hw; theta=0; end T=sqrt((T*sin(theta)-G)^2+(T*cos(theta))^2); if theta<0.01; theta=0; else theta=atan((T*sin(theta)-G)/(T*cos(theta))); end end %若锚链不与海床接触，对浮力进行迭代补偿 while min(datay)>-12 Ff=Ff+1; h=Ff*4/(1.025*1000*9.8*pi*4); Fw=0.625*(2-h)*2*v^2; datafai=[]; fai=atan(Fw/(Ff-Gf)); datafai=[datafai;fai]; T=Fw/sin(fai); for j=1:1:4 fai=atan(Fw/(T*cos(fai)+19.72-98)); T=Fw/sin(fai); datafai=[datafai;fai]; end fai=atan(Fw/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=Fw/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; datax=[];datay=[];dataT=[];datatheta=[]; G=7.203; theta=pi/2-fai; x=0; y=0; fori=1:1:210 datax=[datax;x]; datay=[datay;y]; dataT=[dataT;T]; datatheta=[datatheta;theta]; if y>H-Hw x=x-0.105*cos(theta); y=y-0.105*sin(theta); else x=x-0.105; y=H-Hw; theta=0; end T=sqrt((T*sin(theta)-G)^2+(T*cos(theta))^2); if theta<0.01; theta=0; else theta=atan((T*sin(theta)-G)/(T*cos(theta))); end end end plot(datax,datay); legend('锚链姿态示意图'); grid on Xz=X-datax(210); %求出浮标游动半径 2.问题二的MATLAB算法程序 dataGq=[]; datafai5=[]; Gq=1200*9.8; %Gq初始值为1200*9.8N dataGq=[dataGq;Gq;]; Gt=980;Gg=40*9.8;Gf=9800;Ffg=789;v=36;Hw=18; Ff=Gq+Gt+Gg+Gf-Ffg; h=Ff*4/(1.025*1000*9.8*pi*4); Fw=0.625*(2-h)*2*v^2; datafai=[]; fai=atan(Fw/(Ff-Gf)); datafai=[datafai;fai]; T=Fw/sin(fai); for j=1:1:4 fai=atan(Fw/(T*cos(fai)+19.72-98)); T=Fw/sin(fai); datafai=[datafai;fai]; end datafai5=[datafai5;datafai(5)]; fai=atan(Fw/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=Fw/sin(fai); %对Gq进行迭代计算，使钢桶倾斜的角度小于5° whiledatafai(5)>0.0875 Gq=Gq+100; dataGq=[dataGq;Gq;]; Gt=980;Gg=40*9.8;Gf=9800;Ffg=789;v=36;Hw=18; Ff=Gq+Gt+Gg+Gf-Ffg; h=Ff*4/(1.025*1000*9.8*pi*4); Fw=0.625*(2-h)*2*v^2; datafai=[]; fai=atan(Fw/(Ff-Gf)); datafai=[datafai;fai]; T=Fw/sin(fai); for j=1:1:4 fai=atan(Fw/(T*cos(fai)+19.72-98)); T=Fw/sin(fai); datafai=[datafai;fai]; end datafai5=[datafai5;datafai(5)]; fai=atan(Fw/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=Fw/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; end %依据所确定的Gp算出锚链的形状 datax=[];datay=[];dataT=[];datatheta=[]; G=7.203; theta=pi/2-fai; x=0; y=0; fori=1:1:210 datax=[datax;x]; datay=[datay;y]; dataT=[dataT;T]; datatheta=[datatheta;theta]; if y>H-Hw x=x-0.105*cos(theta); y=y-0.105*sin(theta); else x=x-0.105; y=H-Hw; theta=0; end T=sqrt((T*sin(theta)-G)^2+(T*cos(theta))^2); if theta<0.01; theta=0; else theta=atan((T*sin(theta)-G)/(T*cos(theta))); end end plot(datax,datay); %画出锚链的姿态示意图 legend('锚链姿态示意图'); grid on Xz=X-datax(210); 3.问题二的MATLAB算法程序 l=0.18;G=49.60; %选定锚链型号，确定每节链环的长度和单位长度的重量 h=0.9;%h初始值设为0.9 Ff=10045*pi*h;Gt=980;Gg=40*9.8;Gf=9800;Ffg=789;v=36;Hw=18;vv=1.5; Gq= Ff-(Gt+Gg+Gf-Ffg); Fw=0.625*(2-h)*2*v^2; Fs=374*(2-h)*2*vv^2; Fh=Fw+Fs; datafai=[]; fai=atan(Fh/(Ff-Gf)); datafai=[datafai;fai]; T=Fh/sin(fai); for j=1:1:4 Fss=841.5*cos(fai); fai=atan((Fh+Fss)/(T*cos(fai)+19.72-98)); T=(Fh+Fss)/sin(fai); datafai=[datafai;fai]; end Fss=841.5*cos(fai); fai=atan((Fh+Fss)/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=(Fh+Fss)/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; %对h进行迭代计算，使钢桶倾斜的角度小于5° whiledatafai(5)>0.087 h=h+0.1; Ff=10045*pi*h; Gq= Ff-(Gt+Gg+Gf-Ffg); Fw=0.625*(2-h)*2*v^2; Fs=374*(2-h)*2*vv^2; Fh=Fw+Fs; datafai=[]; fai=atan(Fh/(Ff-Gf)); datafai=[datafai;fai]; T=Fh/sin(fai); for j=1:1:4 Fss=841.5*cos(fai); fai=atan((Fh+Fss)/(T*cos(fai)+19.72-98)); T=(Fh+Fss)/sin(fai); datafai=[datafai;fai]; end Fss=841.5*cos(fai); fai=atan((Fh+Fss)/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=(Fh+Fss)/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; end %依据所确定的h算出锚链的形状 datax=[];datay=[];dataT=[];datatheta=[]; theta=pi/2-fai; x=0; y=0; fori=1:1:500 datax=[datax;x]; datay=[datay;y]; dataT=[dataT;T]; datatheta=[datatheta;theta]; if y>H-Hw x=x-l*cos(theta); y=y-l*sin(theta); else x=x-l; y=H-Hw; theta=0; end T=sqrt((T*sin(theta)-G)^2+(T*cos(theta))^2); if theta<0.01; theta=0; else theta=atan((T*sin(theta)-G)/(T*cos(theta))); end end while min(datay)>-11 %若锚链不与海床接触，对浮力进行迭代补偿 Ff=Ff+1; h=Ff*4/(1.025*1000*9.8*pi*4); Fw=0.625*(2-h)*2*v^2; datafai=[]; fai=atan(Fw/(Ff-Gf)); datafai=[datafai;fai]; T=Fw/sin(fai); for j=1:1:4 fai=atan(Fw/(T*cos(fai)+19.72-98)); T=Fw/sin(fai); datafai=[datafai;fai]; end fai=atan(Fw/(T*cos(fai)+710.039-1300*9.8)); datafai=[datafai;fai]; T=Fw/sin(fai); X=0; H=0; for k=1:1:5 X=sin(datafai(k))+X; H=cos(datafai(k))+H; end H=H+h; datax=[];datay=[];dataT=[];datatheta=[]; theta=pi/2-fai; x=0; y=0; fori=1:1:500 datax=[datax;x]; datay=[datay;y]; dataT=[dataT;T]; datatheta=[datatheta;theta]; if y>H-Hw x=x-l*cos(theta); y=y-l*sin(theta); else x=x-l; y=H-Hw; theta=0; end T=sqrt((T*sin(theta)-G)^2+(T*cos(theta))^2); if theta<0.01; theta=0; else theta=atan((T*sin(theta)-G)/(T*cos(theta))); end end end i=find(datay<-11.6); xx=datax(1:min(i)); yy=datay(1:min(i)); plot(xx,yy) Xz=X-min(xx); %算出游动半径 %假设风力最大时刚好没有拖地部分，此时的L为可取的最小锚链长度 L=0; fori=1:1:500 ifdatatheta(i)>0.05 L=L+l; end end