全国大学生数学建模大赛国家一等奖优秀论文系泊系统的设计

摘要

的设计

本文详细对系泊系统的各个机构进行了力学分析，针对系泊系统的要求，建立优化模型，求解系泊系统在多种环境下的最优解，使得浮标游动范围，吃水程度和钢桶倾斜角度尽可能的小。

针对问题一，本文对系泊系统的受力及力矩进行了分析，基于浮标倾斜的考虑，得到了平衡状态下关于受力平衡及力矩平衡的方程组。由于方程组数量较多及相互影响的特点，直接求解十分困难。因此我们考虑以浮标两边的浸水长度h1,h2为变量，利用搜索算法对方程组进行求解，并得到相应的结果。如当风速为12m/s时，钢桶的倾斜角度1.0405°，从上到下钢管的倾斜角度分别为1.0086°、1.0146°、1.0206°、1.0267°，浮标吃水深度0.735m，浮标游动区域半径14.4429m。

针对问题二，首先将风速为36m/s的情况代入问题一建立的模型中，但是得到的结果不满足题目所给定的要求。则考虑在重物球质量一定的条件下，以浮标的吃水深度和游动区域及钢桶的倾斜角为目标，建立了一个单决策变量的多目标最优系泊模型，相比于问题一，此问的变量更多，更加难于求解，故考虑将多目标转化成单目标的问题进行求解，并继续使用搜索法对问题进行求解。最后找到了三组可行解，其中最优解是重力球的质量为2102kg.

针对问题三，本文中有三个决策变量以及三个变系数，相比于前两问，无论是计算量还是计算维数，难度更大。为了求解该问，建立了一个多决策变量的多目标变系数的最优系泊系统模型，为了简便运算，我们建立了变步长的搜索算法，并最终求解得到结果，得到的一组解为:选用了III型号的锚链，重物球质量为2800kg，锚链长度为23.4m。

针对论文的实际情况，对论文的优缺点做了评价，文章最后还给出了其他的改进方向，以用于指导实际应用。

关键词：系泊系统设计；力的平移定理；多目标；优化模型；搜索算法

1．问题的重述

一个由浮标系统、系泊系统和水声通讯系统组成的近浅海观测网的传输节点。可以简化看作是一个浮于海平面的圆柱体浮标通过四节钢管链接装有通讯设备的钢桶，钢桶再通过锚链链接一个可以移动的锚，锚沉在海床上。为了不让锚被拖行要求锚链末端与锚的链接处的切线方向与海床的夹角不超过16度，并且要保证通讯设备的工作效果，钢桶与电焊锚链链接出可悬挂重物体，使得钢桶的倾斜角度（钢桶与竖直线的夹角）不超过5度。

需要建立模型讨论并解决下列的几个问题：

1. 若海水静止，在选用II型电焊锚链22.05m，重物球的质量为1200kg，布放水

深为18m，海床平坦，海水密度为1.025×103kg/m3的条件下，分别计算海面风速为12m/s和24m/s时钢桶和各节钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

2. 在问题1的假设下，计算海面风速为36m/s时钢桶和各节钢管的倾斜角度、锚

链形状和浮标的游动区域。并调节重物球的质量，使得钢桶的倾斜角度不超过5度，锚链在锚点与海床的夹角不超过16度。

3. 若布放海域的实测水深介于16m-20m之间。布放点的海水速度最大可达到 1.5m/s、风速最大可达到36m/s。给出考虑风力、水流力和水深情况下的系泊系统设计，分析不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

2．模型的假设

（1）假设两钢管用铰链链接在一起，可以自由转动。 （2）假设风是沿平行海平面的方向吹来。 （3）假设忽略锚链和重物球的浮力影响。

（4）假设忽略由于其它原因（如温度，湿度，受力等）而对各个系统产生的形变。

（5）假设浮标不会出现因为风力过大而被吹翻的情况。 （6）假设各个系统之间连接处的长度忽略不计。

3．符号说明

Fwind：浮标所受水平方向的风力大小；

FbuoyantA：浮标所受浮力大小； FbuoyantB：每节钢管所受浮力大小； FbuoyantCFB0：钢桶所受浮力大小；

：第一节钢管对浮标拉力的大小(或浮标对第一节钢管拉力的大小)；

FBi：第i节钢管下端所受拉力的大小(i1,2,3,4)；

FC(FD0)：钢桶所受锚链的拉力大小(或锚链上方第一段所受钢桶的拉力大小)；

FDj：锚链第j段(自上而下)下端所受拉力大小(j1,2,L,N，其中N为张紧状态

锚链的段数)；

Ff：锚与海底间的摩擦力；

Fmax：锚与海底间的最大静摩擦力；

FwaterA：浮标所受的水流力大小；

FwaterBi：第i节钢管所受水流力大小（i1,2,3,4）；

FwaterC：钢桶所受水流力大小；

：锚与海底间的最大静摩擦系数；

FN：海底对锚的支持力；

GA：浮标的重力；

GB：每节钢管的重力；GC：表示钢桶的重力；

GD：锚链每段链环的重力；GE：锚的重力；

0

B：第一节钢管对浮标拉力与水平方向的夹角；

B：第i节钢管下端所受拉力与水平方向的夹角(i1,2,3,4)；

iC(D)：钢桶所受锚链的拉力与水平方向的夹角；

0

Dj：锚链第j段(自上而下)下端所受拉力与水平方向的夹角(j1,2,L,N)；

A：浮标的水平倾斜角度；

B：第i节钢管的水平倾斜角度(i1,2,3,4)；

ic：钢桶的水平倾斜角度；

h1,h2：圆柱体浮标两端的浸水长度,其中h1h2；

h：浮标的吃水深度； vwind：水平方向的风速；

Swind：浮标在风向法平面的投影面积；

VoverflowA：浮标的排水体积；

Vwater：水流速度大小；

mA：浮标的质量；

mB：每节钢管的质量； mC：钢桶系统的质量；

m：重物球的质量；

：海水密度；

D：锚链的线密度；

CDN：法向阻力系数；

CDN：切向阻力系数；

rA：浮标的底面半径；

rB：每节钢管的底面半径；LB：每节钢管的长度；

rC：钢桶的底面半径； ：钢桶的长度；

LCLD：锚链每段链环的长度；

H：海水深度；

s：锚链某点至锚链上端顶点的弧长；

SA：浮标浸入水下部分在水流力法平面的投影面积；

SBi：第i节钢管在水流力法平面的投影面积（i1,2,3,4）；

SC：钢桶在水流力法平面的投影面积；

T：锚链某点处所受拉力的大小；

：锚链某点处所受拉力与水平方向的夹角；

XA：平衡状态下浮标的水平距离；

XB：平衡状态下钢管的水平距离； XC：平衡状态下钢桶的水平距离； XD：平衡状态下锚链的水平距离；

X：平衡状态下整体的水平距离； yD：锚链竖直方向上的高度；

max：锚不滑动时，锚链在锚点与海床的最大夹角； D：单位长度上锚链所受的法向水流力大小；

D：单位长度上锚链所受的切向水流力大小；

4.模型的建立与求解

4.1问题的分析

4.1.1系统平衡的两种状态——张紧状态和松弛状态及游动区域的理解

本文考虑的在风力及水流力作用下的系泊系统的设计问题，首先我们考虑一种简单情形，即海水静止且风向水平情况，此时整个系泊系统可能会因风力移动达到一个平衡的状态，由于风向恒定我们可以认为最终的平衡状态中浮标、钢管、钢桶、锚链和锚大致位于同一个平面内，如图1所示。

图1传输节点示意图

在平衡状态下，各系统所受力及力矩分别达到平衡，此时浮标会与锚形成一定的水平距离。此时的平衡状态可能分为两种：一种是张紧状态，即锚链的每一段都受到了前后端的拉力从而使整个系统平衡；另一种是松弛状态，即锚链上方的一部分存在拉力，而下方部分散落在水底呈无拉力状态达到平衡。

在张紧状态下达到平衡时，浮标到达某个定点，此时该点与锚得水平距离为半径以锚和浮标的水平距离为半径、锚所在位置为圆心做一个圆面，该区域即浮标的游动区域，如图2所示。但在松弛状态下由于松弛部分锚链的可伸缩性，浮标的位置会处于一定的区域，当松弛部分锚链长度最大时浮标位于最远位置，以此时浮标与锚的水平距离为半径、锚所在位置为圆心的圆形区域即其游动区域。

图2浮标游动区域俯视示意图

4.1.2浮标的倾斜问题及对吃水深度的理解

在对组成各个系统的零件受力分析的过程中，会出现无法将物体看作是均匀的情况，这样一来，物体的各力的受力点将不再是质点，为此应尽可能将受力点移动到质心上。这时，根据力的平移定理，在将各力平移到质心上去的过程中会产生力矩，且对于一个平衡的物体其力矩和应为零。

在力矩平衡的作用下，我们发现平衡状态下浮标会产生一定的倾斜角度。 假设在浮标的底面与海平面平行，即没有倾角时，其受力分析如图3所示。

图3浮标垂直海平面的受力分析图

由图可知除浮力FbuoyantA外，风力Fwind和钢管拉力FB都会产生非零力矩，且方向相

0同，导致力矩和不为零，从而会产生旋转，不会达到平衡状态，我们认为在平衡状态下，浮标底面相对海平面存在一个倾斜角度A。

吃水深度为浮标底部到海面的距离。由于浮标倾斜，导致浮标各点侵入海水的长度不同，如图4所示。假设浮标两边侵入水中的长度分别为h1,h2，其中h1h2，则吃水深度即两边浸水长度的平均值

h1h2。 2图4倾斜浮标的浸水示意图

4.2问题1模型的建立

根据对问题的分析发现，问题旨在讨论在不同情形达到平衡状态下，分析各钢管的倾斜角度、锚链形状、浮标的方程，通过已知的信息，求得系统各部分所受拉力大小及其角度情况和各部分本身的倾斜角度。 4.2.1浮标的力平衡方程

通过浮标的受力分析可得，浮标受到竖直向上的浮力FbuoyantA，水平向右的风力Fwind，重力GA以及钢管对其的拉力FB，其受力分析图如图5所示：

0图5浮标的受力分析图

根据平面力系平衡原理，可得浮标的力平衡方程如下： 在水平方向上，有

FB0cosB0Fwind,（1）

其中，B为钢管对浮标拉力与水平方向的夹角，Fwind为水平风力的大小，有如

0下的近似公式（S为物体在风法平面的投影面积）

Fwind0.625S(vwind)2;（2）

在竖直方向上，有

FB0sinB0mAgFbuoyantA,（3）

其中，有

FbuoyantAgVoverflowA,（4）

而

11VoverflowArA2h1rA2(h2h1)rA2(h1h2),（5）

22在这儿，mA为浮标质量，h1,h2分别表示浮标两边的浸水长度(见图4)，rA为浮标的底面半径，VoverflowA为浮标的排水体积，为海水密度，g为重力加速度。 4.2.2浮标在风向法平面的投影面积公式

在计算风力Fwind时，需用到浮标在风向法平面的投影面积S，由于浮标倾斜，该投影区域并非一标准的矩形，而是由两个部分构成，其中的下方部分为矩形，而上方部分为倾斜圆在水平面上投影的一半，即半椭圆的一半，最终得到浮标在风向法平面的投影面积如图6所示。

图6浮标在风向法平面的投影面积示意图

通过平面几何关系可得，第一部分矩形的面积S1满足：

S1=2rA[LA0.5(h1h2)]2rA(h1h2)(2rA)S222,（6）

第二部分半椭圆的面积S2满足：

(h2h1)rA4,（7）

其中，rA为浮标的底面半径，h1,h2分别表示浮标两边的浸水长度，LA表示浮标的长度。

4.2.3钢管的平衡方程

通过分析可得，每节钢管的受力分析图如图7所示。

图7第i节钢管的受力分析图

根据平面力系平衡原理及力的平移定理(力矩平衡公式)，有如下方程： 水平方向力平衡方程

FBicosBiFBi1cosBi1(i1,2,3,4),（8）

竖直方向力平衡方程

FBisinBimBgFBi1sinBi1FbuoyantB(i1,2,3,4),（9）

力矩平衡方程

11FBisin(BiBi)LBFBi1sin(Bi1Bi)LB(i1,2,3,4).（10） 22在这儿，FbuoyantB为每节钢管所受浮力大小，满足

FbuoyantBgrB2LB,（11）

而FB为第i节钢管下端所受拉力大小，mB为每节钢管的质量，g为重力加速度，B是

iiB是拉力FB与水平方向的夹角，rB第i节钢管与水平方向的夹角（钢管的倾斜角度），

ii为钢管的底面半径，LB为钢管的长度。 4.2.4钢桶系统的平衡方程 钢桶的受力分析如图8所示。

图8钢桶的受力分析图

根据平面力系平衡原理及力的平移定理(力矩平衡公式)，有如下方程： 水平方向力平衡方程

FCcosCFB4cosB4,（12）

竖直方向力平衡方程 力矩平衡方程

FCsinC(mCm)gFB4sinB4FbuoyantC,（13）

在这儿，FbuoyantC为钢桶系统的浮力，有

111FCsin(CC)LCFB4sin(B4C)LCmgcosCLC,（14）222

FbuoyantCgrC2LC,（15）

而FB表示钢管对钢桶拉力的大小，B表示钢管对钢桶拉力与水平方向的夹角，

44FC表示钢桶受下方锚链拉力的大小，C表示FC与水平方向的夹角，mC表示钢桶的质

量，m表示重物球的质量，rC表示钢桶的底面半径，LC表示钢桶的长度。 4.2.5锚链系统的平衡方程

由于锚链无档普通链环，我们可以将其看作无弹性悬垂线，这样可以将锚链视作柔性的。因此可用微元法的思想分析其受力平衡状态。

考虑锚链其中自上而下一小段弧长(s到sds)上的受力情况，如图9所示。

图9锚链微元的受力分析图

图中ds表示取的任意一段小弧长，T表示上段锚链对所取小弧长的拉力，而

TdT表示下段锚链对小弧长的拉力；是拉力T与水平方向的夹角；d是拉力TdT与水平方向的夹角。

根据这一弧长微元上的受力平衡，可得出以下方程： 水平方向上，有

Tcos(TdT)cos(d),（16）

竖直方向上，有

Tsin(TdT)sin(d)Ddsg（17）

其中D表示单位长度锚链的质量。 将其展开可得

TcosTcoscosddTcoscosdTsinsinddTsinsind（18） TsinTsincosddTsincosdTcossinddTcossinddsgD(d)3L可知，当角度变化d很小时，忽略高阶无穷小由泰勒公式sindd3!的影响，故sind近似等于d，cosd近似等于1，从而可得如下的近似方程：

0dTcosTsinddTsind（19） dsgdTsinTcosddTcosdD忽略高阶无穷小dTd项，（其中dTd都为无穷小，故它们的乘积为高阶的无

穷小）对上式化简可得如下一阶非线性微分方程组：

DgdcosdsT（20） dTgsinDds且有初值条件(0)C,T(0)FC 4.2.6锚的平衡方程及平衡状态条件

锚的受力直接关系着最终是否达到平衡状态，其受力分析如图10所示：

图10锚的受力分析图

水平方向上，有

FDNcosDNFf,（21）

在竖直方向上，有

FNFDNsinDNGEmEg,（22）

由于判断整个系统平衡的关键即对锚的水平拉力是否超过锚与海底之间的最大静摩擦力

Fmax，即得到平衡状态的判断条件为

FDNcosDNFmax.（23）

又由海水深度为18m，可知

4h1h2HsinALBsin(Bi)LCsinCyD（24）

2i1其中，H表示海水深度；

h1h2sinA表示浮标底面轴心到海平面的距离，h1,h2分24i1i别表示浮标两边浸入水中的长度，A浮标的倾斜角度，LBsin(B)表示4个钢管在水中的高度和，LB表示每节钢管的长度，B表示第i节钢管的倾斜角度；LC表示钢桶的长度，C表示钢桶的倾斜角度，yD表示锚链竖直方向的高度。

i4.2.7游动区域的分析

经过分析和求解可以知道在平衡状态下的各个系统水平距离分别为： 浮标的水平距离

(h2h1)2(2rA)2h2h1XA,（25） 2222(h2h1)(2rA)h2h1钢管的水平距离

XBLBicosBi,（26）

i14钢桶的水平距离

XCLCcosC,（27）

锚链的水平距离

XDCNDs0cos(d)dsd（28）

则整体的水平距离(游动半径)为

XXAXBXCXD（29）

4.3问题一模型的求解 4.3.1模型的求解思路

通过对系泊系统各部分进行单独受力分析和力矩分析，得到了系统各处拉力大小

及其倾角等未知量在平衡状态下的关系，由于方程数量较多，且相互影响，使得直接计算非常复杂。因此我们首先考虑各方程之间的关系，得到未知量满足一定规律的递推公式，对问题进行系列简化，然后再采用搜索算法进行求解。 4.3.2方程的简化及递推公式的推导

4

由于锚链是无档普通链环，实际中是环环相扣的离散系统，为了计算方便，我们首

先考虑对锚链系统对应的一阶非线性微分方程组(20式)进行转换[，由ds与水平方向

3]

及竖直方向上的位移之间的关系

dxcosds,dysinds,（30）

可得链环连接处的受力大小及角度的递推公式

FDj1cosDj1FDjcosDj(j1,2,L,N),（31） FDj1sinDj1FDjsinDjGDj其中，FD表示链环上方第一段所受钢桶带来拉力的大小，D表示链环上方第一段所

0

0受钢桶带来拉力与水平方向的角度，由于其与链环对钢桶拉力FC互为反作用力，有

FD0FC,D0C；FDj表示链环第j段(自上而下)下端所受拉力大小，Dj表示链环第j段

下端所受拉力与水平方向的夹角(j1,2,L,N)，N为链环的段数。 4.3.2.2水平力平衡公式的推导及简化

整个系泊系统的每个部分均在两个水平力的作用下平衡，有

FwindFBicosBiFCcosCFDjcosDj(i1,2,3,4;j1,2,L,N),（32）

其中风力Fwind由浮标两边的入水深度h1,h2决定。 4.3.2.3竖直方向力平衡公式的推导及简化

由竖直方向上受力平衡可得到如下递推公式

FB0sinB0FbuoyantAGAkkFBksinBkFB0sinB0FbuoyantAGAFbuoyantBiGBi(k1,2,3,4)i1i144（33）,其FCsinCFB0sinB0FbuoyantAGAFbuoyantBiGBiFbuoyangtCGCi1i144FDksinDkFB0sinB0FbuoyantAGAFbuoyantBiGBiFbuoyangtCGCkGD(k1,2,L,N)i1i1中浮标的浮力FbuoyantA也由浮标两边的入水深度h1,h2决定。 4.3.2.4系统各段所受拉力大小及其角度的计算公式

将各拉力竖直方向分力与水平分力相比，可得相应角度正切值tan的计算公式如

tanB0(FbuoyantAGA)/FwindkktanBk(FB0sinB0FbuoyantAGAFbuoyantBiGBi)/Fwind(k1,2,3,4)i1i144（34）,tanC(FB0sinB0FbuoyantAGAFbuoyantBiGBiFbuoyangtCGC)/Fwindi1i144tanDk(FB0sinB0FbuoyantAGAFbuoyantBiGBiFbuoyangtCGCkGD)/Fwind(k1,2,L,N)i1i1从而可得平衡状态下依赖于h1,h2的各段拉力角度及大小的计算公式。 4.3.2.5各个部件在竖直方向的高度的计算

（1）找出浮标的竖直方向高度，这里记为HA(HAQR)。如图所示：

图11.浮标的竖直方向高度示意图

图中，h1MN，h2OP，HAQR。 在直角RMQ和直角NKP中，因为

，所以这两个三角形另外两个角相等，即： RMQNPK（同位角相等）

RQMKNParctanKP（35） NK所以

HARQQMcosRQMh1h2PKcos(arctan())（36） 2KN即:

HArA(h1h2)(h1h2)(2rA)22（37）

其中，rA表示浮标的半径,为1m。

（2）钢管的竖直方向高度，这里记为HB。

i根据浮标的平衡方程、水平方向的方程（FBcosBFwind）和竖直方向的方程

00（FBsinBGAFbuoyantA）相除得：

00Fwind0.625S(vwind)2（38） cotB0FbuoyantAGAgVoverflowAGA即可以求出B，这样再带入浮标水平方向平衡方程则可求出FB。再把B和FB的

0000值带入钢管的平衡方程，我们可以求出钢管的FB拉力及其倾斜角度B和钢管的倾斜

ii角B，又可以根据钢管的几何关系（如图6所示）：

iHBHBilBisin(Bi)（39）

ii44其中HB表示四根钢管的总高度，HB表示第i节钢管的高度，B表示第i节钢管的倾斜角度。

（3）钢桶的竖直方向高度，这里记为HC。

ii把上个环节的计算出来的结果带入钢桶系统的平衡方程，同样也可以计算出钢桶倾斜角度C，锚链对钢桶的拉力FC及其方向（与X轴的夹角C）。则可以得出钢桶高度的公式：

HClCsinC（40）

其中lC为钢桶的长度，它是一个常数，大小为1m。 （4）钢桶的竖直方向高度，这里记为HD。

根据已求得的钢桶对锚链的拉力FC及其方向，把这对初始值带入锚链平衡微分方程，可解出函数关系。就可以表示出锚链的高度：

22.05HDdssin(d)（41） 0DNC其中，ds表示锚链一微段的长度，d表示那一微段的倾斜角度。

故整个系泊系统总高度为

HHAHBHCHD（42）

4.3.2.6吃水深度h1和h2的值范围的确定

根据题目可知浮标的高度为2m，则有：0h12m，0h22m，h2为较长的一条高。 只有当浮标在海水中时，浮标的重力与浮标的浮力二力平衡，即：

GAgVoverflowA（43）

可以计算出此时h1h20.31m。那么可以进一步确定范围：

0.31h1h22，（44）

4.3.3算法步骤

（1）给定一组浮标吃水深度值[h1,h2],，（得到0.31h1h22），可以依次计算得到H，D；

N（2）判断如果同时满足H18，（取较小值，比如0.1）和0D16两个条件，

N则输出结果，并停止计算；如果不满足，则取步长有h1h1h1,h2h2h2,并转到步骤（1）。 4.3.4计算结果及分析

（1）当风速为12ms时，得到的结果如下表所示：

表1.风速为12ms时的求解结果 h1 h2 海水深度 锚链水平倾角 钢桶的倾角 所用无档链环个数 0.68m 0.79m 17.9042m 3.4838° 1.0405° 144 钢管1的倾角 钢管2的倾角 钢管3的倾角 钢管4的倾角 游动区域半径 浮标的吃水深度 1.0086° 1.0146° 1.0206° 1.0267° 14.4429m 0.735m 得到锚链的形状如图12所示： 图12.风速为12ms时锚链的形状

以上求解得到的结果满足题目的要求，且所用的链环个数为144个，小于链的总个数，说明达到平衡时锚链并未完全绷紧。

（2）当风速为24ms时，得到的结果如下表所示：

表2.风速为24ms时的求解结果 h1 0.71m h2 0.79m 海水深度 17.9117m 锚链水平倾角 11.7855° 钢桶的倾角 3.9299° 所用无档链环个数 185 钢管1的倾角 钢管2的倾角 钢管3的倾角 钢管4的倾角 游动区域半径 3.814° 3.8356° 3.8575° 3.8797° 17.2252m 浮标的吃水深度 0.75m 得到锚链的形状如图13所示： 图13.风速为24ms时锚链的形状

以上求解得到的结果满足题目的要求，且所用的链环个数为185个，小于链的总个数，说明达到平衡时锚链并未完全绷紧。 4.4问题2的模型分析与建立

4.4.1风速为36ms时系泊状态及其分析

当风速达到36ms时，若不考虑角度的，可通过问题一的递推式及搜索算法求得如下结果：

表3.风速为36ms时的求解结果 h1 h2 海水深度 锚链水平倾角 钢桶的倾角 所用无档链环个数 0.759m 0.7810m 17.7796m 17.7796° 8.1268° 210 钢管1的倾角 钢管2的倾角 钢管3的倾角 钢管4的倾角 游动区域半径 浮标的吃水深度 7.99° 7.9424° 7.9853° 8.0287° 18.3712m 0.77m 得到锚链的形状如图14所示： 图14.风速为36ms时锚链的形状

此时，锚链末端与锚的链接处的切线方向与海床的夹角为17.7796度，超过了16度，钢桶的倾斜角度(钢桶与竖直线的夹角)为8.1268度，超过5度，已经不能够达到平衡状态，因此需考虑调节重物球质量，使其重新达到平衡状态。 4.4.2基于重物球质量的定常系统最优系泊模型

当风速达到36ms时，系泊系统已经不再满足角度的要求，从而不再处于平衡状态，因此我们考虑调整重物球质量。在此我们假设海水静止，考虑水平风速、水深、锚链型号、长度等均为定常数，基于重物球质量的选择，建立了单决策变量的多目标最优系泊模型。 4.4.2.1目标函数

（一）浮标的吃水深度尽可能小：

minh;（45）

m

其中吃水深度hh1h2，即浮标两边浸水长度的平均值； 2（二）游动区域尽可能小，即游动半径尽可能最小：

minmX（46）

其中游动半径XXAXBXCXD(见29式)； （三）钢桶的倾斜角度尽可能小；

minm2C（47）

其中C表示钢桶与水平面的夹角。 4.4.2.2约束条件

(一)海水的深度约束，即平衡状态下整个系泊系统的竖直高度之和应等于水深H，有

4h1h2sinALBsin(Bi)LCsinCyDH,（48） 2i其中

h1h2sinA表示浮标底面轴心到海平面的距离，h1,h2分别表示浮标两边浸入24i1i水中的长度，A浮标的倾斜角度，LBsin(B)表示4个钢管在水中的高度和，LB表示每节钢管的长度，B表示第i节钢管的倾斜角度；LC表示钢桶的长度，C表示钢

i桶的倾斜角度，yD表示锚链竖直方向的高度。

（二）锚链底端拉力的夹角约束，即锚链在锚点与海床的夹角不超过16度，有

(s)Dmax,（49）

N由问题1可知，T,,s满足如下微分方程组

dDgcosdsT,（50） dTgsinDds其中D表示锚链在锚点与海床的夹角，max表示使得锚不滑动时，锚链在锚点

与海床的最大夹角。

（三）平衡约束条件，即达到系统各部分力和力矩平衡的各等式约束条件，包括式？

综上：我们建立了基于重物球质量的多目标最优系泊模型如下

Nminh;minX;min（51）

2C;4h1h2H2sinALBsin(Bi)LCsinCyDi（52） maxDN力及力矩平衡约束条件(1)-(23)4.5问题2的模型求解

4.5.1多目标的转换

对于问题2，这是一个多目标的优化问题，不能够直接进行求解，需要把多目标转化为单目标进行求解。具体方法为：

先计算单目标

min2C（53）

4h1h2H2sinALBsin(Bi)LCsinCyDi（54） maxDN力及力矩平衡约束条件(1)-(23)得到倾角的最优值*，再将*5作为约束条件加入到优化问题里面，这样

2问题转化为双目标优化问题。

minminh;（55） X;4h1h2H2sinALBsin(Bi)LCsinCyDiDNmax(56)*52 力及力矩平衡约束条件(1)-(23)再考虑吃水深度的优化，利用类似的方法得到其最优值作为约束条件，即：

minX;(57)

4h1h2H2sinALBsin(Bi)LCsinCyDiDNmax(58)*5*2 hh2力及力矩平衡约束条件(1)-(23)4.5.2可行解的范围

因为该问所给的风速变大，在初始的条件下不能再得到可行解了。所以这里除了给定h1、h2的值，还需要给定重物球的质量。在问题一种重力球质量为1200kg已经不能满足题目要求，所以重物球的重力应该满足：GC112000N。为了找到重物球重力的上限，首先假设只有重物球的重力作用，且浮标完全淹没在水中。此时，对浮标进行受力分析，得到的受力分析图如图15所示：

图15.浮标的受力分析图

根据受力分析，当浮标处于平衡状态时有：

FbuyantAGAGC1(59)

其中，FbuoyantAgVoverflowA； 由于浮标完全被淹没，所以有： 则：

GC153082.6N。此时，因此设定重物球重力的取值范围是12000NGC153082.6N。

当重物球重力的值接近上限值时，浮标接近被完全淹没，此时浮标受风吹的面积小，

即风力小，整个系统在水平方向只受风力和锚上的摩擦力的作用，且这两个力是一对平衡力，当风力较小时，系统将不会在水平方向上移动，当浮标完全沉没时，系统将完全不在水平方向移动，这时钢桶倾斜角度达到最小为0度，锚链处于松弛状态，且此时为平衡状态。 4.5.3问题2的具体算法步骤

求解的基本思路同问题一是类似的，还是采用搜索算法。但是由于变量个数增加了一个，搜索的范围和计算量将大大增加，如果还采用问题一的定步长搜索的话，计算时间和结果可能难以保证。故采用变步长搜索算法，具体步骤为：

（1）在搜索参数的有效范围内，给定一个较大的步长G,h1,h2（比如[5,0.1,0.1]），与问题一的求解方法一样，得到对应的H,D；

N（2）判断如果同时满足条件H18，（取较小值，比如0.1）和0D16，

N； 则输出结果并结束计算；如果不满足，则记下最接近判断条件的一组G,h1,h2为新的初始点，（3）以G,h1,h2选定一个较小的步长G,h1,h2（比如，[1,0.01,0.01]）

按照上述方法再次计算H,D，如果满足判断条件则停止；否则可以类似方法继续缩

N小步长，总可以找到可行解。 4.5.4问题2的计算结果

根据上述的求解思路和方法，编写程序进行搜索求解，限定搜索的时间，得到以下结果：

第一种结果：

重物球的质量 海水深度 1968kg 1m 1m 17.9081m 锚链水平倾角 钢桶的倾角 所用链环数/总链环数 游动区域半径 13.7633° 4.3475° 207/210 18.2717m 表4.问题2的第一种求解结果 h1 h2 在得到以上结果的情况下，将参数代入可以得到锚链的形状为： 图16.第一种求解结果下锚链的形状

其他结果：

重物球的质量 海水深度 1969kg 1m 1m 17.9038m 锚链水平倾角 钢桶的倾角 所用链环数/总链环数 游动区域半径 12.7107° 4.3475° 210/210 18.3094m 表5.问题2的第二种求解结果 h1 h2 重物球的质量 海水深度 2102kg 1.04m 1.04m 17.9047m 锚链水平倾角 钢桶的倾角 所用链环数/总链环数 游动区域半径 12.0327° 3.9415° 210/210 18.1739m 表6.问题2的第三种求解结果 h1 h2 问题2的第二、第三种求解结果满足题目的要求，说明对于第二问的解题思路和方法是正确且可行的，如果按照此种方法就能找到最优的解，但由于算法的数据处理量过于庞大以及时间的，这里就不一一进行求解了。 4.6问题3的模型建立与求解 4.6.1最优系泊问题

问题一和问题二的模型都是在水深、风速恒定、水流静止的特殊情况下建立的，而对于此问，题目给定水深、风速、水流的范围，考虑实际系泊系统中水流速度的影响，在水速、风速、水深不定的情况下，基于锚链型号、长度及重物球质量的选择，建立了多决策变量的多目标最优系泊模型。

所以决策目标和第二题相同，考虑基于锚链型号、长度及重物球质量的变系数，那么模型的条件因不同情况而不同，下面以当水平风速和水流力同向的情况分析。

4.6.2问题3当水平风速和水流力同向的情况 1浮标的力平衡方程

根据理论力学的相关知识，对平衡时的浮标进行受力分析，得到浮标的受力分析图如图17所示：

图17.浮标的受力分析图

浮标受到竖直向上的浮力FbuoyantA；水平向右的风力Fwind；浮标所受水流力FwaterA；竖直向下的浮标重力GA；浮标受下方钢管的拉力FB。其中A表示浮标的倾斜角度，

0B表示拉力FB与X轴方向的夹角。根据平面力系平衡原理，x轴与y轴方向上受力

00达到平衡，可得出浮标的平衡方程：

水平方向：

FB0cosB0FwindFwaterA(60)

其中，有近海风载荷的近似公式：（S为物体在风法平面的投影面积）

Fwind0.625S(vwind)2(61)

S2rA[LA0.5(h1h2)]2rA(h1h2)(2rA)22+(h2h1)rA4(62)

近海水流力的近似公式：（SA为物体在水流法平面的投影面积）

由问题一中浮标在风向法平面投影面积的分析可知，投影面积可分为矩形和椭圆两

部分，可得：

SA2rA2(h1h2)(h1h2)(2rA)22rA(h2h1)4(63)

竖直方向：

FB0sinB0GAFbuoyantA()

其中，有

GAmAg，FbuoyantAgVoverflowA(65)

11VoverflowArA2h1rA2(h2h1)rA2(h1h2)(66)

22tanAh2h1(67) 2rA其中，h1、h2分别表示浮标左、右边进入水中桶壁的长度 2钢管的平衡方程

对第i节钢管对其进行分析，得到受力分析图如图18所示：

图18.第i节钢管的受力分析图

FbuoyantBi为第i节钢管的在水中的浮力；FBi1为第i节钢管受到上方物体的拉力；FBi：为

第i节钢管受到下方物体的拉力；钢管所受水流力FwaterB；GB为第i节钢管的重力；iii是第i节钢管与X轴的夹角（钢管的倾斜角度）；B是拉力FB与X轴的夹角；B是拉

iii1力FB与X轴的夹角。根据平面力系平衡原理和力的平移定理，x与y方向上受力平衡且力矩和为零，可得出钢管的平衡方程：

i1力平衡：(水平方向)

FBicosBiFBi1cosBi1FwaterBi(i1,2,3,4)(68)

其中，有

FwaterBi374SBi(vwater)2(69)

SB2rBLBsinBirB2cotBisinBi(70)

力平衡：(竖直方向)

FBisinBiGBiFBi1sinBi1FbuoyantBi(i1,2,3,4)(71)

其中，有

GBimBig，FbuoyantBigVBigrB2LB(i1,2,3,4)(72)

力矩平衡：

11FBisin(BiBi)LBFBi1sin(Bi1Bi)LB(i1,2,3,4)(73) 223钢桶系统的平衡方程

对于通信系统的钢桶，得到受力分析如图19所示：

图19.钢桶的受力分析图

节钢管对钢桶的拉力；B表示拉力FB与X轴的夹角；FC表示

钢桶受下方锚链的拉力；钢管所受水流力FwaterC；C表示FC力与X轴的夹角；GC1表示重球对钢桶的重力；GC表示钢桶的自身所受的重力；FbuoyantC表示钢桶在水中所受的浮

44FB4表示上面第4

力；C为钢桶的倾斜角度（钢桶与X轴的夹角）。根据平面力系平衡原理和力的平移定理，x与y方向上受力平衡且力矩和为零，可得出浮标的平衡方程：

力平衡：(水平方向)

FCcosCFB4cosB4FwaterC(74)

其中，有

FwaterC374SC(vwater)2(75)

SC2rCLCsinCrC2cotCsinC(76)

力平衡：(竖直方向)

FCsinCGCFB4sinB4FoverflowC(77)

其中，有

GC(mCm)g，FbuoyantCgVCgrC2LC(78)

力矩平衡：

111FCsin(CC)LCFB4sin(B4C)LCmgcosCLC(79) 2224锚链系统的平衡分析

假设锚链是无弹性悬垂线，这样可以将锚链简化看作是柔软的，利于模型的简化。再考虑其中一小段弧长(s到sds)上的受力情况，在有水流力的情况下，如图20所示：

图20.锚链的受力分析图

图中ds表示取的任意一段小弧长，T表示下段锚链对所取小弧长的拉力，而TdT表示上段锚链对小弧长的拉力；i是拉力T与水平方向的夹角；id是拉力TdT与水平方向的夹角。

根据这一微段的平衡，可得出以下力系平衡方程： 水平方向上，有

Tcos(TdT)cos(d)Dds(80)

竖直方向上，有

Tsinchaindsg(TdT)sin(d)Dds(81)

其中chain表示单位长度锚链的质量,D和D分别为单位长度上锚链所受的法向和切向水流力。

对上式化简可得如下一阶非线性微分方程组：

chaingdcosDdsT(82) dTchaingsinDds通过查找资料[4]，D和D按下列公式来计算：

D1CDNdVwater2cos2(83) 21DCDNdVwater2sin2(84)

2其中，CDN和CDN分别为法向和切向阻力系数，为水的密度，Vwater为水流速度 综上：我们建立了基于重物球质量的多目标变系数最优系泊模型如下

minh;minX;min(85)

2C;4h1h2H2sinALBsin(Bi)LCsinCyDi(86) maxDN力及力矩平衡约束条件(60)-(84)4.6.3模型求解

在水流速度为1.5m/s，风速为36m/s这种情况下，相比问题二多了两个决策变

量：锚链的型号及其长度，增加了三个变系数：水流速度、风速、水深。所以求解这个模型依然通过给定重物球重力GC1，锚链长度s，锚链型号，h1和h2。用问题二同样的思路去搜索解，最后搜索到结果为：

表7.问题三的第一组解 锚链型号 III 重物球的水流速 h1 h2 海水深度 质量 2800kg 1.5m/s 1.3 1.3 16.0057m 钢管1的倾钢管2的钢管3的钢管4的所用链环数/总角 倾角 倾角 倾角 链环数 6.1184° 6.172° 6.1874° 6.2029° 195/195 表8.问题三的第二组解 锚链型号 IV 重物球的水流速 h1 h2 海水深度 质量 2960kg 1.5m/s 1.4m 1.4m 18.7723m 钢管1的倾钢管2的钢管3的钢管4的所用链环数/总角 倾角 倾角 倾角 链环数 5.5703° 5.6178° 5.6306° 5.34° 131/168 锚链水平钢桶的倾倾角 角 5.284° 6.4744° 游动区域锚链长度 半径 23.4m 21.1156m 锚链水平倾角 0.0331° 锚链长度 25.2m 钢桶的倾角 5.887° 游动区域半径 26.5694m 得到的锚链形状如图所示：

图21.第一组求解结果下锚链的形状 图22.第二组求解结果下锚链的形状

5．模型的推广与改进方向

本文建立了单一决策变量的多目标以及多决策变量的多目标优化模型，推导出系泊系统各个部分的平衡方程，巧妙的运用搜索算法，搜索算法极大的减少了运算量，提高了运算效率。根据建立的模型可以求解并设计一个最优的系泊系统。

本文求解模型的方法可以推广应用到各种参数较多，方程数较多的方程组求解。以及本文建立的模型可以很好的运用于各种链状的物体受力分析的研究上去。但本文对更多优解的研究还不够深入，可以继续求解更有结进行研究。

6．模型的优缺点

通过对系泊系统各个部分严谨的受力分析，得出了比较完美的受力平衡关系式，进而建立了比较完善的模型，思路严晰，结构严谨。只是由于运算量太大，没有找到更好的求解方法。

如果时间足够充裕还可以对模型的求解进行更加深入的研究。

7.参考文献

[1]姜启源.数学模型（第三版）[M].北京：高等教育出版社，1999.

[2]韩中庚.数学建模方法及其应用（第二版）[M].北京：高等教育出版社，2009. [3]孙宁松，海上移动式平台锚泊定位系统锚索链受力分析[J].中国海洋平台，

2008,23(2)

[4]赫春玲，滕斌，不均匀可拉伸单锚链系统的静力分析[J].中国海洋平台，

2003,18（4）

[5]滕斌，赫春玲，韩凌，Chebyshev多项式在锚链分析中的应用[J].中国工程

科学，2005,7（1）

8．附录

1、第一问和第二问程序 %系泊系统的第一问和第二问 %搜索方法 clear clc closeall

g=9.8;%重力加速度 v=36;%风速 rou=1025;%海水密度 ma=1000;da=2;ha=2;%浮标

lb=1;mb=10;db=0.05;%钢管 mc=100;lc=1;dc=0.3;%钢桶 ball=1200;%重物球

ld=0.105;md=7;Ld=22.05;nd=Ld/ld;md=md*ld*g;%链环 me=600;%锚

Ffb=rou*g*pi*db^2/4*lb;%钢管的浮力 mb=mb*g;%钢管的重力 t=1;

forball=1200%1200:2:3000

forh1=0.759%0.9:0.01:1.1%0.3:0.05:1%0.31:0.01:1.5%0.68:0.0005:0.72 forh2=0.781%h1:0.01:1.5%0.75:0.0005:0.79 L=h1; R=0;

%------浮标-------------------------

S1=da*(ha-(h1+h2)/2)*da/sqrt((h2-h1)^2+da^2); S2=pi*da^2/8*(h2-h1)/da; S=S1+S2;%风载荷 Fw=0.625*S*v^2;%风力

V=pi*da^2/8*(h1+h2);%吃水体积 Ffa=rou*g*V;%浮标的浮力 Ta1=Ffa-ma*g;%竖直方向 Ta2=Fw;%水平方向

Ta=sqrt(Ta1^2+Ta2^2);%浮标的总拉力 thetaa=atand(Ta1/Ta2);%浮标拉力的角度

alphaa=atand((h2-h1)/da);%浮标的倾斜角 %------浮标------------------------- %------钢管------------------------- fork=1:4

Tb2(k)=Ta2;%水平方向 ifk==1

Tb1(k)=Ta1+Ffb-mb;%竖直方向 else

Tb1(k)=Tb1(k-1)+Ffb-mb;%竖直方向 end

thetab(k)=atand(Tb1(k)/Tb2(k));%拉力角度 Tb(k)=sqrt(Tb1(k)^2+Tb2(k)^2);%拉力大小 ifk==1%钢管的倾斜角

alphab(k)=atand((Ta*sind(thetaa)+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Ta*cosd(thetaa))); else

alphab(k)=atand((Tb(k-1)*sind(thetab(k-1))+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Tb(k-1)*cosd(thetab(k-1)))); end ifk==1

xb(k)=ld*cosd(alphaa);%钢管k的x长度 else

xb(k)=ld*cosd(alphab(k-1));%钢管k的x长度 end

R=R+xb(k);%目前游动总长度 L=L+sind(alphab(k))*lb; k=k+1; end

%------钢管------------------------- %------钢桶------------------------- Ffc=rou*g*pi*dc^2/4*lc;%钢桶的浮力 Tc2=Ta2;%水平方向

if(Tb1(4)+Ffc-(mc+ball)*g)<0 break; end

Tc1=Tb1(4)+Ffc-(mc+ball)*g;%竖直方向 Tc=sqrt(Tc1^2+Tc2^2);%钢桶的总拉力 thetac=atand(Tc1/Tc2);%钢桶拉力的角度

alphac=atand((Tc*sind(thetac)+Tb(4)*sind(thetab(4))+ball*g)... /(Tb(4)*cosd(thetab(4))+Tc*cosd(thetac)));%钢桶的倾斜角 %ifalphac<85 %continue; %end

R=R+lc*cosd(alphac);%目前游动总长度 L=L+lc*sind(alphac);

%------钢桶------------------------- %------链环------------------------- %fork=1:nd

%ifk==1

%Td(k)=Tc-md*sind(thetac);%拉力

%thetad(k)=(-md*cosd(thetac))/Tc+thetac;%拉力角度 %else

%Td(k)=Td(k-1)-md*sind(thetad(k-1));%拉力

%thetad(k)=-md*cosd(thetad(k-1))/Td(k-1)+thetad(k-1);%拉力角度 %end

%options=odeset('RelTol',1e-4,'AbsTol',[1e-81e-8]); %ifk==1

%[A,B]=ode45(@rigid,[0,0.105],[Tcthetac/180*pi],options); %else

%[A,B]=ode45(@rigid,[nd*(k-1)nd*k],[Td(k-1),thetad(k-1)/180*pi],options); %end

%Td(k)=B(end,1);

%thetad(k)=B(end,2)*180/pi; fork=1:nd

Td2(k)=Tc2;%水平方向拉力 ifk==1

Td1(k)=Tc1-md;%竖直方向拉力 else

Td1(k)=Td1(k-1)-md;%拉力 end

thetad(k)=atand(Td1(k)/Td2(k));%拉力角度

Td(k)=sqrt(Td1(k)^2+Td2(k)^2);%拉力大小 ifk==1

xd(k)=0-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=0-ld*sind(thetad(k));%锚链k的y坐标 else

xd(k)=xd(k-1)-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=yd(k-1)-ld*sind(thetad(k));%锚链k的y坐标 end

R=R+ld*cosd(thetad(k));%目前游动总长度 L=L+ld*sind(thetad(k));%

ifabs(L-18)<0.1&&thetad(k)>0&&thetad(k)<30%第一二问 disp('d'); numk(t)=k; BALL(t)=ball; THETA(t)=thetad(k); H1(t)=h1;H2(t)=h2; t=t+1; break; else continue; end

%ifthetad(k)<0

%L=L-ld*sind(thetad(k)); %ifabs(L-18)<0.1

%disp('d'); %H1=h1;H2=h2; %break; %else %continue; %end %end k=k+1; end

%------链环------------------------- end end end %plot(Td); %figure %plot(thetad)

[thetasort,s]=sort(THETA);

H1best=H1(s(1));H2best=H2(s(1));Ballbest=BALL(s(1)); numk=numk(s(1)) plot(xd,yd) 2、第三问程序： %系泊系统的第三问 %搜索方法 clear

clc closeall

g=9.8;%重力加速度 v=36;%风速 vs=1.5;%水流速度 rou=1025;%海水密度 ma=1000;da=2;ha=2;%浮标 lb=1;mb=10;db=0.05;%钢管 mc=100;lc=1;dc=0.3;%钢桶 %ball=1200;%重物球 lmd=[0.0783.2; 0.1057; 0.1212.5; 0.1519.5; %ld=0.105;md=7; %Ld=22.05;

%nd=Ld/ld;md=md*ld*g;%链环 me=600;%锚

Ffb=rou*g*pi*db^2/4*lb;%钢管的浮力 mb=mb*g;%钢管的重力 t=1; fori=3%1:5

ld=lmd(i,1);md=lmd(i,2)*ld*g; forLd=23.4%22:ld:30

nd=Ld/ld;

forball=2800%2600:40:3000%1200:2:3000

forh1=1.3%1.1:0.1:1.8%0.9:0.01:1.1%0.3:0.05:1%0.31:0.01:1.5%0.68:0.0005:0.72

forh2=h1%h1:0.01:1.5%0.75:0.0005:0.79 L=h1; R=0;

%------浮标-------------------------

S1=da*(ha-(h1+h2)/2)*da/sqrt((h2-h1)^2+da^2); S2=pi*da^2/8*(h2-h1)/da; S=S1+S2;%风载荷 Fw=0.625*S*v^2;%风力

Sa=da^2*(h1+h2)/2/sqrt((h1-h2)^2+da^2)+pi*da*(h2-h1)/8;%浮标的水受力面积 Fsa=374*Sa*vs^2;%水流力 V=pi*da^2/8*(h1+h2);%吃水体积 Ffa=rou*g*V;%浮标的浮力 Ta1=Ffa-ma*g;%竖直方向 Ta2=Fw+Fsa;%水平方向

Ta=sqrt(Ta1^2+Ta2^2);%浮标的总拉力 thetaa=atand(Ta1/Ta2);%浮标拉力的角度 alphaa=atand((h2-h1)/da);%浮标的倾斜角 %------浮标------------------------- %------钢管------------------------- fork=1:4

ifk==1

Sb=db*lb*sind(thetaa)+pi*db^2/4*cotd(thetaa)*sind(thetaa); Fsb=374*Sb*vs^2;%水流力 Tb2(k)=Ta2+Fsb;%水平方向 Tb1(k)=Ta1+Ffb-mb;%竖直方向 else

Sb=db*lb*sind(alphab(k-1))+pi*db^2/4*cotd(alphab(k-1))*sind(alphab(k-1));

Fsb=374*Sb*vs^2;%水流力 Tb2(k)=Ta2+Fsb;%水平方向 Tb1(k)=Tb1(k-1)+Ffb-mb;%竖直方向 end

thetab(k)=atand(Tb1(k)/Tb2(k));%拉力角度 Tb(k)=sqrt(Tb1(k)^2+Tb2(k)^2);%拉力大小 ifk==1%钢管的倾斜角

alphab(k)=atand((Ta*sind(thetaa)+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Ta*cosd(thetaa))); else

alphab(k)=atand((Tb(k-1)*sind(thetab(k-1))+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Tb(k-1)*cosd(thetab(k-1)))); end ifk==1

xb(k)=ld*cosd(alphaa);%钢管k的x长度 else

xb(k)=ld*cosd(alphab(k-1));%钢管k的x长度 end

R=R+xb(k);%目前游动总长度 L=L+sind(alphab(k))*lb; k=k+1; end

%------钢管------------------------- %------钢桶------------------------- Ffc=rou*g*pi*dc^2/4*lc;%钢桶的浮力

Sc=dc*lc*sind(alphab(4))+pi*dc^2/4*cotd(alphab(4))*sind(alphab(4)); Fsc=374*Sc*vs^2;%水流力 Tc2=Tb2(4)+Fsc;%水平方向 if(Tb1(4)+Ffc-(mc+ball)*g)<0 break; end

Tc1=Tb1(4)+Ffc-(mc+ball)*g;%竖直方向 Tc=sqrt(Tc1^2+Tc2^2);%钢桶的总拉力 thetac=atand(Tc1/Tc2);%钢桶拉力的角度

alphac=atand((Tc*sind(thetac)+Tb(4)*sind(thetab(4))+ball*g)... /(Tb(4)*cosd(thetab(4))+Tc*cosd(thetac)));%钢桶的倾斜角 %ifalphac<85 %continue; %end

R=R+lc*cosd(alphac);%目前游动总长度

L=L+lc*sind(alphac);

%------钢桶------------------------- %------链环------------------------- %fork=1:nd %ifk==1

%Td(k)=Tc-md*sind(thetac);%拉力

%thetad(k)=(-md*cosd(thetac))/Tc+thetac;%拉力角度 %else

%Td(k)=Td(k-1)-md*sind(thetad(k-1));%拉力

%thetad(k)=-md*cosd(thetad(k-1))/Td(k-1)+thetad(k-1);%拉力角度 %end

%options=odeset('RelTol',1e-4,'AbsTol',[1e-81e-8]); %ifk==1

%[A,B]=ode45(@rigid,[0,0.105],[Tcthetac/180*pi],options); %else

%[A,B]=ode45(@rigid,[nd*(k-1)nd*k],[Td(k-1),thetad(k-1)/180*pi],options); %end

%Td(k)=B(end,1);

%thetad(k)=B(end,2)*180/pi; fork=1:nd

Td2(k)=Tc2;%水平方向拉力 ifk==1

Td1(k)=Tc1-md;%竖直方向拉力

else

Td1(k)=Td1(k-1)-md;%拉力 end

thetad(k)=atand(Td1(k)/Td2(k));%拉力角度 Td(k)=sqrt(Td1(k)^2+Td2(k)^2);%拉力大小 ifk==1

xd(k)=0-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=0-ld*sind(thetad(k));%锚链k的y坐标 else

xd(k)=xd(k-1)-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=yd(k-1)-ld*sind(thetad(k));%锚链k的y坐标 end

R=R+ld*cosd(thetad(k));%目前游动总长度 L=L+ld*sind(thetad(k));%

%ifabs(L-18)<0.1&&thetad(k)>0&&thetad(k)<16%第一二问 ifL<=20&&L>=16&&thetad(k)>0&&thetad(k)<16%第三问 disp('d'); numi(t)=i; LD(t)=Ld; numk(t)=k; BALL(t)=ball; THETA(t)=thetad(k); Alphac(t)=alphac; H1(t)=h1;H2(t)=h2;

t=t+1; continue; end

%ifthetad(k)<0

%L=L-ld*sind(thetad(k)); %ifabs(L-18)<0.1 %disp('d'); %H1=h1;H2=h2; %break; %else %continue; %end %end k=k+1; end

%------链环------------------------- end end end end end %plot(Td); %figure %plot(thetad)

[thetasort,s]=sort(THETA); forj=1:length(s) ifAlphac(s(j))>85 H1best=H1(s(j)) H2best=H2(s(j)) Ballbest=BALL(s(j)) numk=numk(s(j)) num=numi(s(j)) ld=LD(s(j)) break; end end plot(xd,yd)