2022年高教社杯全国大学生数学建模竞赛A 题代码和思路

基于解析和数值方法的波浪能发电装置最大功率求解

摘要：

本文研究了波浪能发电装置的运动状态与最大平均输出功率，根据牛顿第二定律和刚体转动定律，分别对浮子与振子建立常微分方程模型，根据阻尼系数是否恒定采用不同方法求解(角)位移与(角)速度：以平均输出功率为目标函数.根据阻尼系数是否恒定求解最优化问题，得到最大输出功率及相应参数。

问题一，首先对平衡状态下的装置进行受力分析，得到浮子与振子初始位置，并将其作为参照点，然后对垂荡运动下的浮子与振子分别进行受力分析，根据牛顿第二定律，可以得到波浪能装置垂荡运动状态下的二元二阶常系数非齐次微分方程组。

在直线阻尼器的阻尼系数恒定的情况下，采用拉普拉斯变换法，对微分方程组求拉普拉斯变换，求解出浮子振子的位移与速度的象函数，再作拉氏反变换回到时域得到微分方程组的解析解，求得第100s时浮子位移为-0.0836m，速度为-0.643m/s,振子位移为-0.084m,速度为-0.6042m/s,其他时刻数据见正文。

在阻尼系数非恒定的情况下，将微分方程改写成四元一阶非常系数微分方程组，采用四阶 Runge- Kutta法求微分方程组的数值解，得到100s时浮子位移为-0.0894m,速度为-0.6066m/s,振子位移为-0.0952m, 速度为-0.6467m/s,其他时刻数据见正文。

问题二，基于问题一中给出的微分方程组，在阻尼系数恒定的情况下，由于电学方程和力学方程具有对偶性质，采用力电比拟的方法，将力学量等价为电学量，力学方程比拟成电学方程，画出电路图，使用电路理论和阻抗匹配原理，将阻尼器的输出功率转化为电阻的最大功率进行求解，求得当阻尼系数为37193.81N·s·m⁻¹时, PTO系统的最大平均功率为229.334W。

在阻尼系数非恒定的情况下，建立优化模型计算阻尼器的功率的极值，以比例系数和幂指数为决策变量，取速度稳定后的十个周期的平均功率作为目标函数，浮子和振子的运动方程作为控制函数，浮子和振子的运动限制为约束条件，建立最优化模型。采用梯形法计算目标函数的数值积分，求解最优化模型得到比例系数为100000, 幂指数为0.3975, PTO的最大平均输出功率为228.89W。

问题三，采用从整体到局部的方法分析纵摇运动，首先根据浮子质量均匀的条件，计算出浮子的质心位于圆锥底部上方1.564m处，再由垂荡运动下振子的坐标得到系统整体质心的坐标。由此表示出质心与转轴的距离，并计算出相对于转轴的转动惯量。考虑力矩作用点不同，根据假设对力矩进行了不同参考系下的转化。参考垂荡运动的微分方程，根据刚体转动定理列出角位移的微分方程，并与垂荡运动的方程联立，采用四阶 Runge- Kutta法求解，得到第100s时相对于转轴,浮子垂荡位移为-0.0501m,垂荡速度为-0.9466m/s,纵摇角位移为-0.0022, 纵摇角速度为0.0401s⁻¹;           振子垂荡位移为-0.0426m，垂荡速度为-1.0365m/s,纵摇角位移为-0.00024, 纵摇角速度为0.0527s⁻¹。

问题四，根据问题三列出的关于浮子与振子平动与转动的方程，参考问题二列出的最优化模型，目标函数变为直线阻尼器和旋转阻尼器的功率之和，控制函数为浮子和振子的转动方程，再次使用梯形法求得当直线阻尼器的阻尼系数为62140N·s/m,旋转阻尼器的阻尼系数为95541N·s²时,PTO的最大平均输出功率为327.96W。

关键字：波浪发电装置 力电比拟法 拉普拉斯变换 最优化 龙格库塔法 梯形法

一、问题重述

1.1 问题背景

随着人类科技和社会经济的迅速发展，传统化石能源已经不能满足发展的需要，能源需求和环境污染将是人类面临的双重挑战，因而，发展可再生能源产业受到了世界各国的重视。波浪能作为一种储量大、无污染的可再生能源，具有很大的开发潜力。本题给出了一种振荡浮子式波浪发电装置，其由浮子、振子、中轴及能量输出系统(PTO) 组成，在波浪的作用下，浮子在运动的同时带动振子运动，两者的相对运动能够驱动阻尼器做功，并将所做的功作为能量输出。

1.2 问题要求

问题1仅考虑垂荡，在波浪激励力fcosωt的作用下，要求分别计算两种情况下浮子和振子在前40个波浪周期内，时间间隔为0.2s的垂荡位移和速度。情况1为直线阻尼器的阻尼系数恒为10000N·s/m，情况2 为阻尼系数与浮子和振子的相对速度绝对值的0.5次幂成正比，且比例系数取10000。

问题2在问题1的基础上，给出阻尼范围为[0，100000]，幂指数区间为[0，1].要求在阻尼恒定与变化两种情况下，计算最大输出功率及相应的最优阻尼系数。

问题3要求考虑垂荡与纵摇，在波浪激励力fcosωt和波浪激励力矩Lcosωt的作用下，要求计算直线阻尼器和旋转阻尼器的阻尼系数分别恒定为 10000N.s/m和1000N·m·s时，浮子和振子在前40个波浪周期内，时间间隔为0.2s的垂荡位移与速度和纵摇角位移与角速度。

问题4在问题3的基础上，给出直线阻尼器和旋转阻尼器的阻尼系数的取值范围均为[0，100000]，要求计算最大输出功率及相应的最优阻尼系数。

二、问题分析

2.1 对问题一的分析

首先对平衡状态下的浮子与振子进行受力分析，可以求得平衡状态下浮力、弹簧弹力、浮子与振子重力之间的关系，求出平衡时浮子的吃水深度与弹簧压缩量。

然后对垂荡运动下的浮子与振子进行受力分析，浮子在线性周期微幅波作用下会受到重力、浮力、波浪激励力、兴波阻尼力、弹簧弹力以及阻尼器产生的阻尼力，振子受到重力、弹簧弹力与阻尼器产生的阻尼力。

由于在平衡状态时，浮子的重力，浮力和弹簧拉力已经处于平衡状态；在垂荡状态时，我们又只关心浮子相对平衡位置的位移，故可将浮子的重力，浮力和一部分弹簧拉力合成为静水恢复力，对振子也可做相似简化。

根据牛顿第二定律，可以分别对浮子与振子列动力学方程。根据各个力的定义，列出的方程组为二阶非齐次线性微分方程组。在阻尼系数恒定时，方程组是常系数的：在阻尼系数与浮子振子相对速度有关时，即阻尼系数非恒定时，方程组是非常系数的。问题一需要解这两个微分方程组。

对于二阶常系数非齐次线性微分方程组，其激励是正弦函数，解析解容易求得，可以用拉普拉斯变换，将时域的微分方程组变换到复频域中得到象函数，解象函数的方程，再做拉普拉斯反变换恢复回时域。

对于二阶非常系数非齐次微分方程组，其解析解难以求得，因此可以采用微分方程的数值解法，譬如龙格-库塔法。数值解法是通用的解法，因此还可以应用在问题二的优化模型中。但是数值解法是近似的，可能与解析解存在误差，因此需要对其进行误差分析，证明其可靠性。

2.2对问题二的分析

问题二是一个优化问题，要求求阻尼器的最大功率。阻尼器的瞬时功率可以表示为阻尼系数乘以相对速度的平方，平均功率为瞬时功率在一个周期内的积分除以周期长度。和问题一一样，问题二分为阻尼系数恒定与阻尼系数非恒定两种情况。

阻尼系数恒定时，理论上可以再次使用拉普拉斯变换求解含参的微分方程组，但此方法略繁琐。考虑到力学系统和电学系统间存在对偶性，我们可以将遵循牛顿第二定律的二阶常系数非齐次线性微分方程组中的质量换成电感，阻尼系数换成电阻，弹簧刚度换成电感的倒数，位移换成电荷，便可得到对偶的遵循回路电流方程的二阶常系数非齐次线性微分方程组。这种方法被称为力电比拟法，力学问题化为了电学问题。

通过方程组，可以画出对应的电路图，特别是由于阻尼被比拟成电阻，阻尼的功率就是电阻的功率，因此只需要计算电阻阻值取值为何值时其功率最大，可以利用阻抗匹配的知识求解。

阻尼系数非恒定时，无法通过力电比拟法求解，因此只能列出优化模型，决策变量是阻尼系数前的比例系数与幂指数，用阻尼系数和相对速度表示出目标函数。因为目标函数是一个定积分，而被积函数是通过微分方程的数值解求得，只有离散点的取值，因此需要采用数值积分法，譬如梯形法才能求得目标函数的取值。

2.3对问题三的分析

问题三要求装置除了垂荡运动外，还作纵摇运动，因此需要一组新的变量来描述其运动。将系统视为一个质量不变，但是质心位置变化的刚体：根据刚体运动的规律，可将系统的运动分为平动和转动，分别可以对应垂荡和纵摇运动。

为了便于分析，作出假设，即浮子和振子进行纵摇运动的角位移很小，因而在垂荡运动时浮子和振子近似处在同一条直线上，纵摇运动不会影响到竖直方向上的受力，在此种假设下垂荡运动仍可以沿用第一问的建立的模型求解。而纵摇运动，可以根据刚体的转动规律来求解。只要计算出装置的各个部分相对于转轴的转动惯量，根据已知的力矩条件就可以得到关于角位移的微分方程。初始位置，浮子、振子、转轴、质心均处于同一直线上。

2.4对问题四的分析

问题四可以参照问题二建立最优化模型。问题四浮子与振子的运动和问题三一致，因此决策变量是直线阻尼和旋转阻尼的阻尼系数；从问题三可以得到浮子

三、   模型假设

1. 不考虑中轴、底座、隔层及 PTO 的质量和各种摩擦。

2. 假设海水是无粘、无旋且不可压缩的理想流体，可采用微幅波理论进行计算。

3. 假设运动过程中圆锥部分始终全在水下。

4. 假设直线阻尼器和旋转阻尼器的阻尼力(矩)转化为电能的效率是100%。

5. 除了弹簧和阻尼器，装置上的所有零件均不发生形变，可以视作刚体。

6. 假设浮子和振子进行纵摇运动的角位移较小，在垂荡运动时浮子和振子近似处在同一条直线上，且纵摇运动不会影响分析垂荡运动的受力。

四、   符号说明

五、     模型的建立与求解

5.1问题一模型的建立

由于浮子和振子之间通过阻尼器所产生的阻尼力与弹簧弹力相衔接，故将浮子与振子分别进行受力分析。记竖直向上为正方向。

5.1.1平衡状态下的受力分析

首先，对平衡状态下的浮子进行受力分析，令浮子质量为m₁，振子质量为m₂，可以得到：

 Ffloot-m1g+FN=0                  (1-1)

 Ffloot=ρgV0       &n