2022国赛数模A题思路以及解析（附源码 可供学习训练使用）

需要全部源码或者论文请点赞关注收藏后评论区留言

前言

发布已获得创作队伍的同意，论文最终斩获省一等奖，写的十分优秀，可供后面的数模比赛训练参考

摘要

                基于摇荡模型的波浪能装置最大输出功率设计问题研究

本文研究了波浪能转换装置的最大输出功率设计问题。通过建立摇荡运动模型， 利用拉氏变换、数轴法、质量投影法以及模拟运行计算实现了对波浪能装置最大输出 功率设计的求解。

针对问题一，考虑浮子中轴做垂荡运动，受到波浪激励力、兴波阻力及浮子位移 带来的浮力变化。通过垂荡运动的特性及牛顿第二定律，建立浮子位移与时间的关系 式，对目标关系式进行拉氏变换，求解出浮子的绝对位移：

再根据浮子与振子所构成的弹簧-质量-阻尼系统，利用数轴法对振子受力情况进行分 析，得到振子相对位移：

求得浮子与振子的对应函数后，在两种不同阻尼系数情况下求出浮子与振子的垂荡 位移与速度，将求得的详细数据在附件及论文中给出。

针对问题二，将问题一中所建立的运动模型进行参数的重新取值，得到问题二中 的浮子位移和振子位移，通过 P=FV 得到阻尼器功率：

 

设立起始步长 =20000,终止步长 =100000，模拟步长 20000 并代入函数 进行计算，求解出当阻尼器阻尼系数 C=40000 左右时 PTO 系统的平均输出功率最大。

针对问题三，在浮子所做纵摇运动完全由波浪激励力矩所提供情况下，引入质量 投影法求解得到浮子壳体对于圆锥形壳体顶点轴的转动惯量 J=25083.2。再根据力矩公 式求解出浮子所做纵摇运动的垂荡位移与速度以及纵摇角位移与角速度，将求得的详 细数据在附件及论文中给出。

针对问题四，可利用问题二的结论得到直线阻尼器的功率

2

由

转动功率为力矩和加速度的积得到旋转阻尼器的功率分别在区 间内对旋转、直线两种阻尼器的 c 取不同值代入功率函数绘图，比较功率大小得到近 似的最优阻尼系数分别为 500、40000。

关键词：垂荡运动 纵摇运动 弹簧系统 力矩

问题重述 

一、问题重述

1.1 问题背景

在经济和社会飞速发展的今天，人类正面临着能源需求和环境污染的双重挑战，发 展可再生能源产业已成为世界各国的共识。波浪能作为一种海洋中重要的可再生能源， 分布广泛，储量丰富，具有相当可观的应用前景，而波浪能装置的能量转换效率是波浪 能能否规模化利用的关键问题之一。

1.2 波浪能装置及参数

波浪能装置是由浮子、振子、中轴以及能量输出系统（PTO，包括弹簧和阻尼器） 构成，在波浪的作用下，浮子运动并带动振子运动，通过两者的相对运动驱动阻尼器做 功，并将所做的功作为能量输出。浮子在线性周期微幅波作用下会受到波浪激励力、附 加惯性力、兴波阻尼力和静水恢复力。在分析问题时，忽略中轴、底座、隔层及 PTO 的 质量和各种摩擦。

1.3 问题要求

问题一：初始时刻浮子和振子平衡于静水中，利用附件提供的参数值分别对以下两 种情况计算浮子和振子在波浪激励力 作用下前 40 个波浪周期内时间间隔为 0.2 s 的垂荡位移和速度。 1.直线阻尼器的阻尼系数为 10000 N·s/m 2. 直线阻尼器的阻尼系数与浮子和振子的相对速度的绝对值的幂成正比，其中比例 系数取 10000，幂指数取 0.5。

问题二：仍考虑浮子在波浪中只做垂荡运动，对以下两种情况建立确定直线阻尼器 的最优阻尼系数的数学模型，使得 PTO 系统的平均输出功率最大。 1. 阻尼系数为常量，阻尼系数在区间 [0,100000] 内取值 2. 阻尼系数与浮子和振子的相对速度的绝对值的幂成正比，比例系数在区间 [0,100000] 内取值，幂指数在区间 [0,1] 内取值。计算两种情况的最大输出功率及相应 的最优阻尼系数。

问题三：初始时刻浮子和振子平衡于静水中，考虑浮子只做垂荡和纵摇运动,利用附 件提供的参数值，假定直线阻尼器和旋转阻尼器的阻尼系数均为常量，分别为 10000 N·s/m 和 1000 N·m·s，计算浮子与振子在波浪激励力和波浪激励力矩 ，作 用下前 40 个波浪周期内时间间隔为 0.2 s 的垂荡位移与速度和纵摇角位移与角速度。 给出 10 s、20 s、40 s、60 s、100 s 时，浮子与振子的垂荡位移与速度和纵摇角位移与角 速度。

问题四：考虑浮子在波浪中只做垂荡和纵摇的情形，针对直线阻尼器和旋转阻尼器 的阻尼系数均为常量的情况，建立确定直线阻尼器和旋转阻尼器最优阻尼系数的数学模 型，直线阻尼器和旋转阻尼器的阻尼系数均在区间 [0,100000] 内取值。计算最大输出功 率及相应的最优阻尼系数。

问题分析 

2.1 问题一的分析

由于题设中浮子只做垂荡运动，于是对浮子在竖直 z 轴上的受力情况进行分析，根 据牛顿第二定律可以得到在 的激励下浮子的位移关于时间的函数。再在波浪能 装置内部建立弹簧-质量-阻尼系统数学模型，利用数轴法建立方程能够得到浮子位移和 振子位移的运动关系。根据两种情况对阻尼系数分别代入具体数值求解得出浮子和振 子前 40 个波浪周期内时间间隔为 0.2 s 的垂荡位移和速度。

2.2 问题二的分析

通过问题一所得垂荡运动模型，在问题二中代入新的参数值求出对应浮子垂荡位移 以及振子的相对位移。由对应浮子与振子的运动关系建立 PTO 系统输出功率模型，计 算出在不同阻尼系数的条件下 PTO 系统的输出功率图，并找出当 PTO 平均输出功率达

到最大值时，阻尼系数的详细值。

2.3 问题三的分析

为简化运算假设纵摇运动对垂荡运动的受力情况没有干扰,利用质量投影法算出浮 子外壳在顶点轴上的转动惯量，再对浮子进行力矩分析，得到浮子在纵摇运动下角位移 关于时间的函数。建立浮子与振子之间的力矩关系模型，对振子的角位移以及角速度进 行求解。而浮子在海浪中做垂荡运动时，直线阻尼器的阻尼系数变化对浮子本身的运动 不造成影响，所以在问题三中浮子的受力方程在形式上与问题一保持一致，对所得方程 代入新参数得到浮子与振子的垂荡位移及速度。

2.4 问题四的分析

对于问题四，直接由问题二的结论可得到直线阻尼器的功率函数，替换同一物理量 的参数数值大小即可求解。由物理学定律得旋转阻尼扭矩做功是扭矩和角加速度的积。 便可得到旋转阻尼器的功率函数。

三、模型假设

1.垂荡运动受到外部激励时，浮子与振子的加速度保持一致

 2.浮子做垂荡运动时，受到的外力规律稳定不变；振子的运动对浮子的运动没有影响

 3.纵摇运动对垂荡运动的受力情况没有干扰；垂荡运动对纵摇运动的影响视为对振子的 转动半径增加一个常量。

4.转轴架的高度忽略不计

                                                              符号与变量说明如下图

 

 模型建立与求解

5.1 问题一模型的建立与求解

5.1.1 浮子垂荡运动模型的建立与求解

针对浮子在海浪中做垂荡运动的运动情况，根据牛顿第二定律，列出浮子在受到波 浪激励力 作用下浮子的受力方程如下

其中 为浮子位移，N 为兴波阻尼系数， 为附加质量， 、 分别为浮子质量及 振子质量。 结合附件所给具体参数值，根据已知的初值：初速度为 0，初位移为 0，使用拉普拉 斯变换求解可以得到：  

 

化简后得浮子只做垂荡位移时，位移与时间的关系式以及图像： 

 

 

5.1.2 浮子与振子运动关系模型的建立与求解

浮子受到波浪激励力会导致浮子与振子之间存在相对运动，根据弹簧-质量-阻尼系 统的数轴法可对浮子及振子的受力情况进行分析：当弹簧或阻尼器两端均有位移时，弹 簧力的方向与所建立的数轴正方向相反，通过牛顿第二定律得到振子的受力图。

 

 

 

5.1.3 两种不同阻尼系数情况下浮子与振子的垂荡位移与速度

在直线阻尼器阻尼系数为 10000 的情况下，分别对浮子垂荡模型以及浮子与 振子的运动关系模型进行代入求解，得到浮子与振子的运动模型。将所得具体数值结果 通过 excel 导出至 result 1-1.xlsx。在结果中取出 t=10s、20s、40s、60s、100s 时浮子和振 子的具体位移与速度随附下表:

表 1 不同时间节点中浮子与振子的位移与速度值（阻尼系数为常数）

 

 

 

5.2 问题二模型的建立与求解 5.2.1 浮子和振子垂荡模型的建立和求解

由于浮子在海浪中做垂荡运动时，直线阻尼器的阻尼系数变化对浮子本身的运动不 造成影响，所以在问题二中浮子的受力方程在形式上与问题一保持一致。可以得到浮子 做垂荡运动的运动公式：

 

 

5.2.2 阻尼系数为未知常量时 PTO 系统输出功率模型的建立与求解

由 PTO 系统中阻尼器功率为

 

 

 

5.3.1 浮子纵摇运动模型的建立

对于假设我们有纵摇运动与垂荡运动相互独立，即浮子所做纵摇运动完全由波浪激 励力矩 产生。由附件运动简图我们得到浮子的纵摇运动是绕圆锥壳体顶点，各质 元到刚体转轴距离保持不变的转动。在计算浮子绕圆锥壳体顶点轴的转动惯量时引入质 量投影法，将刚体上各质元向垂直于转轴的平面投影，保持总质量不变，便得到新的刚 体，其转动惯量没有变化。这样，可将三维刚体转化为与其转动惯量相等的二维刚体，

质量投影后的转动惯量计算图如下：

 

 

 

 

5.3.2 振子纵摇运动模型的求解

对装置做纵摇运动时的振子进行分析。在纵摇时，振子不直接受到波浪激励力矩等 的作用，有保持静止的惯性力，浮子平面发生纵摇时，平面与中轴的夹角发生改变， 产生相对角位移，进一步旋转阻尼器和扭转弹簧产生扭矩????和????，受到的运动激励来 自这两者的扭矩。受激励后振子中轴沿一假想的固定轴进行旋转。振子有重力矩与非 匀速转动时产生的力矩????和????，由此得到力矩关系式

由题目信息可知，弹簧扭矩与角位移成正比，故有,阻尼器扭矩 与角速度成正比，角速度是角位移的时间微分，故有  .θ是振子沿轴 的角位移，??0是浮子沿轴进行的角位移，是振子与浮子平面之间的相对角位移， 角位移单位为 rad/s。

振子所受重力沿旋转切向的分量为由物理学知识的重力 矩为受力与力臂的积，计算重力矩时，由于振子为均匀分布的刚体，可等效为一个在 其几何中心的质点，故力臂长度为振子高度的一半与振子底部到转动轴长度的和。转 动轴长度可视为直线弹簧的原长。故计算重力矩时力臂长度为

?? ′ =0.5m+ 1 2 *0.5m=0.75m。得到重力矩?(57.3 ∗ ??1 )。

振子非匀速转动产生的力矩与其转动惯量和角加速度相关。均匀刚体的转动惯量 的计算公式为?? = ∫ ?? 2 0 ，由于旋转半径变化，且振子为均匀分布的刚体，故将其等 效为一质点，故?? = ∫ ?? 2???? ?? 0 = ?? 2 ∗ ??， r 是旋转半径，半径原长为弹簧原长和振子底 面中心到质心的距离，质心在其高度一半处，故半径原长为 0.5m+0.5*0.5m=0.75m。 但由于纵摇过程中伴随垂荡运动，振子到轴心的距离变化，转动半径发生变化具体过 程过于复杂此处无法分析，故半径的变化取一个常量。最后在垂荡伸长范围内取附加 半径大小，得到等效半径为 1.6m。振子旋转运动的角加速度为角速度的时间微分，

由于在实际求解时，式子中存在θ（t）的三角函数超出计算机处理能力，且重力矩相较 于其他力矩数值较小，故省略，得到简化后的式子

 

 

5.3.3 浮子及振子垂荡运动模型的求解

由于浮子在海浪中做垂荡运动时，直线阻尼器的阻尼系数变化对浮子本身的运动不 造成影响，所以在问题三中浮子的受力方程在形式上与问题一保持一致。可以得到浮子 做垂荡运动的运动公式：

 

 问题三后面部分省略

5.4 问题四模型的建立与求解

5.4.1 旋转阻尼器功率模型的建立及其求解

将问题三中得到的代入上式得到含未知常数参量 c 的函数。由于函数复杂， 无法直接计算得出功率最大值的 c，故对范围在[0,100000]内的 c 取数个值，代入功率 函数利用 MATLAB 绘图，比较不同 c 下的功率，得到近似最优阻尼系数。

5.4.2 直线阻尼器功率模型的建立及其求解

由问题三得到， 振子仍独立做垂荡运动。故直线阻尼器功率方程仍与问题二中的 一致，为

 

 

六、模型的优缺点与评价

6.1 模型优点 （1）利用数轴法和质量投影法准确的对方程中的参数进行具体求解，在一定程度 上保证了所求方程的合理性。

（2）将波浪能装置的运动方式分为垂荡和纵摇两个模型单独分析，对运动情况的 讨论由浅入深，从基础定理出发进行公式推导，使模型具有坚实的数学基础。尽管阻 尼器功率公式过于庞杂，但所建立模型抓住主要影响因素，最终对波浪能装置最大功 率的设计求解得到了假设条件下较为精确的结果。

6.2 模型缺点

（1）问题二的模型由于运动方程过于复杂，无法进行拉氏变换和利用计算机直接计 算得出最优解，可使用优化算法处理得到具体的解。

（2）在做垂荡和纵摇运动的分析时，没有计入除了垂荡引起振子运动半径变化以外 的互相影响，且垂荡引起的半径变化没有具体分析，仅做粗略估计，产生的误差较大。

附录代码

附录 1：matlab 绘制 Q1 浮子垂荡运动图像 
>> e=2.718; 
t=0:0.2:100; 
y=-0.4258.*(-cos(1.4.*t)+e.^(-0.038.*t).*cos(1.911.*t))+0.1404*2*sin(1.911.*t)+0.0134*2*sin(1.4.*t); 
plot(t,y); 
ans = zeros(501,3); 
ans(:,1) = t; 
 
 
 
附录 2：matlab 绘制 Q1 振子相对运动图像 
 
 
>> t=0:0.2:100; 
e=2.718; 
y1=0.1192*e.^(-0.038.*t).*sin(1.91.*t)-2*0.0286*e.^(-2.055.*t).*sin(5.35.*t)-0.4498*e.^(-0.038.* t).* 
cos(1.91.* t) + 0.1366 .*e.^(-2.055.* t).*cos(5.35.* t); 
y2=2*0.2164.*cos(1.4 .*t)-2*0.0545.*sin(1.4 .*t)-0.0396*2.*cos(1.91.*t)+2*0.1473.*sin(1.91.*t); 
>> plot(t,y1+y2); 
>> ans = zeros(501,2); 
a
ans(:,2) = y1+y2; 
 
 
 
附录 3：matlab 绘制 Q1 振子绝对运动图像 
>> t=0:0.2:100; 
e=2.718; 
y1=0.1192*e.^(-0.038.*t).*sin(1.91.*t)-2*0.0286*e.^(-2.055.*t).*sin(5.35.*t)-0.4498*e.^(-0.038.* t).* 
cos(1.91.* t) + 0.1366 .*e.^(-2.055.* t).*cos(5.35.* t); 
y2=2*0.2164.*cos(1.4 .*t)-2*0.0545.*sin(1.4 .*t)-0.0396*2.*cos(1.91.*t)+2*0.1473.*sin(1.91.*t); 
y=y1+y2; 
 
ans = zeros(501,3); 
ans(:,1) = t; 
ans(:,2) = y+g; 
ans(2:501,3) = diff(y+g); 
 
 
 
附录 4：matlab 求 Q1 变阻尼近似平均阻值 
>> e = 2.718; 
t = 0:0.2:100; 
y1=0.1192*e.^(-0.038.*t).*sin(1.91.*t)-2*0.0286*e.^(-2.055.*t).*sin(5.35.*t)-0.4498*e.^(-0.038.* t).* 
cos(1.91.* t) + 0.1366 .*e.^(-2.055.* t).*cos(5.35.* t); 
y2=21.4 .*t)-2*0.0545.*sin(1.4 .*t)-0.0396*2.*cos(1.91.*t)+2*0.1473.*sin(1.91.*t); 
y = y1+y2; 
g = diff(y); 
g = abs(g).^(1/2); 
sum = 0; 
 for i = 1:500 
sum = sum+g(i); 
end; 
sum = sum/500 
 
 
 
 
 
附录 5：matlab 绘制 Q1 变阻换平均振子绝对图像 
>> t=0:0.2:100; 
e=2.718; 
y0.038.*t).*sin(1.91.*t)-2*0.0517*e.^(-0.630.*t).*sin(5.7.*t)-0.2384*2*e.^(-0.038.* t).* 
cos(1.91.* t) + 2*0.0265 .*e.^(-0.630.* t).*cos(5.7.* t); 
y2=2*0.2249.*cos(1.4 .*t)-2*0.0271.*sin(1.4 .*t)-0.0129*2.*cos(1.91.*t)+2*0.1568.*sin(1.91.*t); 
g=-0.4258.*(-cos(1.4.*t)+e.^(-0.038.*t).*cos(1.911.*t))+0.1404*2*sin(1.911.*t)+0.0134*2*sin(1.4.*t); 
y = y1+y2+g; 
 plot(t,y); 
>> ans = zeros(501,2); 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
 
 
 
附录 6：matlab 绘制 Q2 浮子垂荡运动图像 
>> e=2.718; 
t=0:0.2:100; 
+e.^(-0.009.*t).*cos(1.93.*t))+0.009*2*sin(1.93.*t)+0.009*2*sin(2.21.*t); 
plot(t,y); 
ans = zeros(501,3); 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
 
 
 
附录 7：matlab 绘制 Q2 振子垂荡相对运动图像 
>> t=0:0.2:100; 
e=2.718; 
y1=0.1192*e.^(-0.038.*t).*sin(1.91.*t)-2*0.0286*e.^(-2.055.*t).*sin(5.35.*t)-0.4498*e.^(-0.038.* t).* 
e.^(-2.055.* t).*cos(5.35.* t); 
y2=2*0.2164.*cos(1.4 .*t)-2*0.0545.*sin(1.4 .*t)-0.0396*2.*cos(1.91.*t)+2*0.1473.*sin(1.91.*t); 
>> plot(t,y1+y2); 
>> ans = zeros(501,2); 
ans(:,1) = t; 
ans(:,2) = y1+y2; 
 
 
 
附录 8：matlab 绘制 Q2 变阻换平均振子绝对运动图像 
>> t=0:0.2:100; 
 
 
e=2.718;: 
((0.
0.009)^2 + 3.725))* 80000/(2433*(s^2) + x*s + 80000) 
y = y1+y2+g; 
 plot(t,y); 
>> ans = zeros(501,2); 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
 
 
 
附录 9：matlab 处理 Q2 功率-阻尼系数问题（1） 
:20000:100000 
syms s 
syms t 
F=[((0.498*s)/(s^2 + 4.8841) + 0.03474/(s^2 + 3.7249) + 0.03978/(s^2 + 4.8841) - (0.4982*(s + 
0.009))/((s + 0.009)^2 + 3.7249))* 80000/(2433*(s^2) + x*s + 80000)]; 
f(t)=ilap
p(t) = x*(diff(f(t))^2); 
ezplot(p,[0,30]); 
hold on; 
end; 
legend('阻尼系数为 20000','阻尼系数为 40000','阻尼系数为 60000','阻尼系数为 80000','阻尼系数为
100000'); 
 
 
 
附录 10：matlab 处理 Q2 功率-阻尼系数问题（2） 
for x=30000:10000:50000 
syms s 
syms t 
syms F 
syms f 
syms p 
syms w 
F=[(0.+s^2) + (0.426*s)/(1.96+(s^2)) + 0.537/(3.651+ (s^2)) - 
(0.426*(0.038+s))/(3.652+((0.038+s)^2))) * 80000/(2433*(s^2) + x*s + 80000)]; 
f(t)=ilaplace(F); 
p(t) = diff(f)*diff(f); 
 
hold on; 
ezplot(x*p(t),[0,30]); 
end; 
 
 
legend('30000','40000','50000'); 
 
 
 
附录 11：matlab 绘制 Q3 浮子垂荡运动图像 
>> e=2.718; 
t=0:0.2:100; 
y=0
0.04.*t).*sin(1.94.*t); 
plot(t,y); 
ans = zeros(501,3); 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
 
 
 
附录 12：matlab 绘制 Q3 浮子纵摇运动图像 
>> e=2.718; 
t=0:0.2:100; 
y=(0.176.*e.^(-0.01.*t).*sin(0.531.*t) - 0.1.* sin(1.7.* t) + 7.406.* e.^(-0.01.* t).* cos(0.531.* t) - 7.406.* 
cos(1.7.* t))/360; 
plot(t,y); 
ans = zero; 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
 
 
 
附录 13：matlab 绘制 Q3 振子垂荡绝对运动图像 
t = 0:0.2:100; 
e = 2.718; 
y =0.00 
t).*2 .*cos(1.94.* t) + 0.0288.* e.^(-0.04.* t)*2 .*cos(1.94.*t) + 0.256.* cos(1.7.* t) - 0.0192 .*sin(1.7 .*t); 
g=0.476.*cos(1.7.*t) - 0.476.*e.^(-0.04.*t).*cos(1.94.*t)- 0.074.*sin(1.7.*t) + 0.074* e.^(-
0.04.*t).*sin(1.94.*t); 
plot(t,y+g); 
ans = zeros(501,3); 
ans(:,1) = t; 
ans(:,2) = y+g; 
ans(2:501,3) = diff(y+g); 
 
 
 
 
 
附录 14：matlab 求解振子纵摇相对运动图像 
>> syms s 
syms z 
f =[((25000 + 1000*s)* (1/360 * (-(7.406 *s)/(s^2 + 2.89) - 0.17/(s^2 + 2.89) + (7.406* (s + 0.01))/((s 
+ 0.01)^2 + 0.281961) + 0.093456/((s + 0.01)^2 + 0.281961))))/(25000+ 1000 *s + 2433*1.6*1.6*s^2)]; 
y = ilaplace(f); 
e = 2.718; 
z=(0.176*e^(-0.01*t)*sin(0.531*t) - 0.1* sin(1.7* t) + 7.406* e^(-0.01* t)* cos(0.531* t) - 7.406* cos(1.7* 
t))/360; 
ezplot(y-z,[0,100]); 
 
 
 
附录 15：matlab 绘制 Q3 振子纵摇绝对运动图像 
e = 2.718; 
t = 0:
y = 0.014.*sin(1.7.*t) - 0.014.* e.^(-0.08.* t).* sin(2.* t) + 0.022.* e.^(-0.01.* t) .*cos(0.531.* t) - 0.07* 
cos(1.7.* t) + 0.048* e.^(-0.08.* t) .*cos(2.* t); 
ans =zeros(501,3); 
ans(:,1) = t; 
ans(:,2) = y; 
ans(2:501,3) = diff(y); 
>> plot(t,y); 
 
 
 
附录 16：matlab 求解 Q4 浮子垂荡运动图像 
>> syms s 
syms t )/((s^2 + 1.98^2)* (8390 *s^2 + 528.5* s + 31541.3))]; 
>> ezplot(y,[0,100]); 
>> p = ilaplace(y); 
>> ezplot(p,[0,100]); 
 
 
 
附录 17：matlab 求解 Q4 浮子纵摇运动图像 
>> syms s 
syms t 
y = [
>> ezplot(f,[0,100]); 
 
 
 
 
 
附录 18：matlab 求解 Q4 垂荡功率-阻尼问题 
>> syms s 
syms t 
for x = 20000:20000:100000 
y = [(17602 + 1.98^2)* (8390 *s^2 + 528.5* s + 31541.3)) * 80000/(2433*(s^2) + x*s + 80000)]; 
尼问题（1） 
>> syms s 
syms t 
syms z 
for x = 2000:2000:10000; 
f =[((25000 + x*s)* (2140 *s)/(s^2 + 1.98^2)*1/((7141.5 + 25082)* s^2 + 1655.9 *s + 8890))/(25000+ 
x*s + 2433*1.6*1.6*s^2)]; 
y = ilaplace(f); 
 
e = 2.718; 
z=(0.176*e^(-0.01*t)*sin(0.531*t) - 0.1* sin(1.7* t) + 7.406* e^(-0.01* t)* cos(0.531* t) - 7.406* 
t))/360; 
 
g= y-z; 
p = x*(diff(g)^2); 
ezplot(p,[0,100]); 
hold on; 
end; 
legen数为 2000','阻尼系数为 4000','阻尼系数为 6000','阻尼系数为 8000','阻尼系数为
10000'); 
 
 
 
附录 20：matlab 求解 Q4 纵摇功率-阻尼问题（2） 
>> syms s 
syms t 
 
syms z 
for x = 500:500:2000; 
f =[((25000 + x*s)* (2140 *s)/(s^2 + 1.98^2)*1/((7141.5 + 25082)* s^2 + 1655.9 *s + 8890))/(25000+ 
x*s 
z=(0.176*e^(-0.01*t)*sin(0.531*t) - 0.1* sin(1.7* t) + 7.406* e^(-0.01* t)* cos(0.531* t) - 7.406* cos(1.7* 
t))/360; 
g= y-z; 
p = x*(diff(g)^2); 
ezplot(p,[0,100]); 
hold on; 
end; 
legend('阻尼系数为 500','阻尼系数为 1000','阻尼系数为 1500','阻尼系数为 2000'); 
>> 

需要全部源码或者论文学习的请点赞关注收藏后评论区留言