基于解析和数值方法的波浪能发电装置最大功率求解
摘要：
本文研究了波浪能发电装置的运动状态与最大平均输出功率，根据牛顿第二定律和刚体转动定律，分别对浮子与振子建立常微分方程模型，根据阻尼系数是否恒定采用不同方法求解(角)位移与(角)速度：以平均输出功率为目标函数.根据阻尼系数是否恒定求解最优化问题，得到最大输出功率及相应参数。
问题一，首先对平衡状态下的装置进行受力分析，得到浮子与振子初始位置，并将其作为参照点，然后对垂荡运动下的浮子与振子分别进行受力分析，根据牛顿第二定律，可以得到波浪能装置垂荡运动状态下的二元二阶常系数非齐次微分方程组。
在直线阻尼器的阻尼系数恒定的情况下，采用拉普拉斯变换法，对微分方程组求拉普拉斯变换，求解出浮子振子的位移与速度的象函数，再作拉氏反变换回到时域得到微分方程组的解析解，求得第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. 假设浮子和振子进行纵摇运动的角位移较小，在垂荡运动时浮子和振子近似处在同一条直线上，且纵摇运动不会影响分析垂荡运动的受力。
四、 符号说明
| 符号 | 说明 | 单位 |
| m₁ | 浮子的质量 | kg |
| m₂ | 振子的质量 | kg |
| x₁ | 浮子位移 | m |
| x₂ | 振子位移 | m |
| F₁ | 静水恢复力 | N |
| F₂ | 修正后的弹簧弹力 | N |
| F₃ | 阻尼器产生的阻尼力 | N |
| F₄ | 垂荡兴波阻尼力 | N |
| F₅ | 波浪激励力 | N |
| k₁ | 垂荡静水恢复力系数 | N/m |
| k₂ | 弹簧刚度 | N/m |
| k₃ | 阻尼器阻尼力系数 | N·s/m |
| k₄ | 垂荡兴波阻尼系数 | N·s/m |
| mc | 垂荡附加质量 | kg |
| x₀ | 平衡时弹簧长度 | m |
| J₁ | 浮子相对于转轴的转动惯量 | N·m² |
| 1₂ | 振子相对于转轴的转动惯量 | N·m² |
| θ₁ | 浮子角位移 | rad |
| θ₂ | 振子角位移 | rad |
| M₁ M₂ | 静水恢复力矩 弹簧扭矩 | N·m N·m |
| M₃ | 阻尼器阻尼力矩 | N·m |
| M₄ | 纵摇兴波阻尼力矩 | N·m |
| M₅ | 波浪激励力矩 | N·m |
| l₁ | 静水恢复力矩系数 | N/m |
| l₂ | 扭转弹簧刚度 | N·m |
| l₃ | 阻尼器阻尼力矩系数 | N·s² |
| l₄ | 纵摇兴波阻尼系数 | N·m·s |
| Jc | 纵摇附加转动惯量 | kg·m² |
五、 模型的建立与求解
5.1问题一模型的建立
由于浮子和振子之间通过阻尼器所产生的阻尼力与弹簧弹力相衔接，故将浮子与振子分别进行受力分析。记竖直向上为正方向。
5.1.1平衡状态下的受力分析
首先，对平衡状态下的浮子进行受力分析，令浮子质量为m₁，振子质量为m₂，可以得到：
Ffloot-m1g+FN=0
(1-1)
Ffloot=ρgV0
(1-2)
其中，Ffloat 为浮力，FN为弹簧对浮子的弹力，V₀为排开水的体积，FN可表示为：
FN=-k2x0
(1-3)
x₀为弹簧此时相对原长的缩短量。
其次，对平衡状态下的振子进行受力分析：
-m2g-FN=0
(1-4)
将(1-1)(1-2)(1-3)式合并有:
ρgV₀-m₁g-k₂x₀=0
(1-5)
相当于将浮子和振子当作一个整体进行受力分析。由(3)(4)可以算出：
x₀=0.2980m
假设浮子的圆锥部分全在水下，令浮子在平衡状态下的吃水深度为Δh，浮子圆柱部分高度和圆锥部分高度分别为H₁和H₂，浮子低半径为R，则可以求得：
V0=13πR2H2+πR2Δh-H2
(1-6)
进而由(1-5)(1-6)可以求得此时浮子吃水深度为：
Δh=2.8m
因为 ah=2.8m>H₂,
故此时圆锥部分的确全在水下，且水下的高度为0.2m
5.1.2垂荡运动状态下的受力分析
设浮子偏离平衡位置的距离为x₁，振子偏离平衡位置的距离为x₂，如图1所示。设竖直向上为正方向。
针对垂荡运动状态下的浮子进行受力分析，如图2所示。首先，浮子受到波浪激励力，记作F₅。F₅对浮子的初始作用方向为竖直向上，其表达式为：
F₅=fcosωt
(1-7)
其中，f为垂荡激励力振幅，ω为角频率。
浮子在海水中做垂荡运动时，会受到浮力 Ffloat，比平衡状态下所受浮力大△Ffloat：
Ffloat=ρgV0+Δffioot
(1-8)
△Frtoat可以用浮子位移为x₁时，浮标在海平面以下的体积V(x₁)计算得到，因此△Ffoat仅为x₁的函数:
ΔFμoat=0gVx1-ρgV0
(1-9)
考虑到浮标的几何形状是圆柱下拼圆锥，且几何参数已知， △Frioat可以计算如下:
ΔFfloot=-ρgπR2x1≤-1-ρgπR2x1-1≤x1<22≤x2
(1-10)
由式(1-10)，假设浮子可以处于全部浸入海面、海面位于圆柱部分两种状态，不考虑其他状态。在求解后面的优化问题中，会将这个条件作为优化问题的约束条件之一。
PTO 系统中包含一个弹簧和一个阻尼器。此时弹簧对浮子产生的弹簧弹力仍然记作FN，比平衡状态下所受弹力大△FN，可得：
FN=-k2x1-x2+x0=-k2x0+ΔFA
(1-11)
将阻尼器对浮子产生的阻尼力记作F₃，因为阻尼器的阻尼力与浮子和振子的相对速度成正比，比例系数为阻尼器的阻尼系数，记作k₃，故：
F3=-k3x1-x2
(1-12)
浮子在海水中做垂荡运动时，会兴起波浪，产生对垂荡运动的阻力，称为兴波阻尼力，记作F₄。F₄与垂荡运动的速度成正比，方向相反，故：
F₄=-k₄x₁
(1-13)
此外，浮子运动时会推动周围海水运动，产生附加惯性力，其对应产生一个附加质量，记作 mc。
根据以上分析，可以得到浮子的受力分析关系式：
Fftout+FN+F3+F4+F5-m1g=m1+mcx1
(1-14)
代入(1-8)与(1-11)式:
pg V0+ΔFfloat-k2x0+ΔFN+F3+F4+F5-m1]g=m1+mcx11-1
5)再代入(1-5)式，消去常数项：
ΔFfloac+ΔFN+F3+F4+F5=m1+mcx1
(1-16)
ΔFμoot
就是由浮体在垂荡运动时所受到的浮力变化引起的静水恢复力，记做 F₁:
FN
为修正后的弹簧对浮子的拉力，记做 F₂,
因此最终有：
F1+F2+F3+F4+F5=n1+mcx1
(1-17)
对垂荡运动状态下的振子进行受力分析，如图3所示。
振子受到弹簧与阻尼器对它的力及其本身的重力，关系式如下：
[--k₂x₂-x₁+x₀+F₃+c₂=m₂x₂
(1-18)
代入(1-3)(1-4)两式, 得到
-F₂+F₃=m₂x₂
(1-19)
综上，可以得到如下微分方程组模型：
|
代入数据得： |
|
|
|
c)x₁ |
|
(1-20) |
m1+mcx1=-k1x1-k2x1-x2-k3x1-x2-k3x1-x2-kix1+fcosωt}1-21
7
5.2问题一模型的求解
上述模型是一个二元二次常系数非齐次微分方程组，问题一给出了两种情况，在第一种情况下，直线阻尼器的阻尼系数为定值，这使得方程组的形式比较简单，可以通过求拉普拉斯变换的方法求得其解析解：在第二种情况下，直线阻尼器的阻尼系数与浮子和振子的相对速度的绝对值的幂成正比，表达式相对复杂，使用拉普拉斯变换计算难度较大，所以使用四阶 Runge- Kutta 法求其数值解。
5.2.1阻尼系数为定值情况下的求解
运用拉普拉斯变换的线性性质和微分性质可将复杂的常微分方程运算过程简单化。微分方程的拉普拉斯变换解法的方法是：
1、先取根据拉氏变换把微分方程化为象函数的代数方程：
2、根据代数方程求出象函数：
3、再取逆拉氏变换得到原微分方程的解。
在阻尼系数为定值的情况下，(1-20)为二阶微分非齐次线性方程组，其象函数容易求出，因此我们采用拉普拉斯变换求解。
先假设在整个运动过程中永远只有一部分圆柱和所有圆锥在水下，即-1≤x₁≤2, 此时:
F₁=-ρgπR²x₁=-k₁x₁
(2-1)
我们会在计算结束后验证这个假设。
将(1-21)式移项并合并同类项：
{m1+mcX1+k¯3+k4x1-k3x2+k1+k2x1-k2x2-fcosωt=0}(2-2
…2 2 -3^1 ' “3 ^2 -2^1 ' “2^2 - ”
根据拉普拉斯变换的定义， 一个定义在区间[0，∞)的函数f(t)，它的单边拉普拉斯变换式F(s)为：
Fs=0- fte-5tdt
(2-3)
其单边拉普拉斯逆变换为：
ft=12πiβ=-jβ+j Fsexds
(2-4)
对式(2-2)求单边拉普拉斯变换得象函数：
代入第一问的相关参数和初值：
解得：
分别对X₁(s),X₂(s),sX₁(s),sX₂(s)做单边拉普拉斯逆变换, 并将结果中的复指数化为三角函数的形式，可得：
xi00-7578×10>98<10γ2-83γ-10≤30x210>10γ2-39γ+10≤53,8≤0<60γ2-62γ2-18γ-14≥8080<61γ2+63γ-16≥69kt=19<52γ2+16γ-69≥798c=(27<10γ2-97γ2+60γ-9w
因为四个变量都只在t≥0的时间不为0，所以每一个变量都乘以阶跃函数：
注意到每一个变量都由形如 Asinωt+ϕeᵇᵗ
的三项组成，其中( ω₁=6.247
与 ω₂=1.881
为此系统的特征频率. ω₃=1.4005
为外力的频率。
取时间步长为0.2s, x₁(t)与x₂(t)-x₁(t). x₁(t)与x₂(t)如图4所示:
可见，在整个运动时间内始终有 -1≤x₁≤2,
说明永远只有一部分圆柱和所有
06
圆锥在水下的假设是正确的。
取t= 10s,20s,40s,60s,100s 计算这些时刻的x₁(t)与x₂(t),x₁(t)与x₂(t)如下表:
| 时间(s) | 浮子 | 振子 | ||
| 位移(m) | 速度(m/s) | 位移(m) | 速度(m/s) | |
| 10 | -0.190713821 | -0.693953476 | -0.211681303 | -0.64100737 |
| 20 | -0.590687694 | -0.272774738 | -0.634251241 | -0.240950128 |
| 40 | 0.285373218 | 0.332912377 | 0.296497416 | 0.312970676 |
| 60 100 | -0.314507575 -0.083615586 | -0.51572832 -0.643002735 | -0.331437723 -0.084068954 | -0.479454391 -0.604211053 |
5.2.2阻尼系数为变值情况下的求解
5.2.2.1 四阶 Runge- Kutta法简介
Runge- Kutta法用于求解形如 y'=fxy
的微分方程组，其中x，y可以是长度一样的向量。基本思想是利用f(x，y)在某些特殊点上的函数值的线性组合，来估算高阶单步法的平均斜率。四阶 Runge-Kutta法使用如下方法进行计算：
|
|
|
(2-9) |
其中
y₁为迭代下标为i 的因变量，x₁为迭代下标为i 的自变量，步长为 h=xᵢ₊₁-
x₁, d₁~d₄为 Runge- Kutta 法中间变量。
5.2.1.3 模型的求解
第二种情况下，直线阻尼器的阻尼系数与浮子和振子的相对速度的绝对值的0.5次幂成正比，其中比例系数取10000，即：
k3=1000|x1-x2|05
(2-10)
代入式(1-20)中，并将二元二阶常系数非齐次微分方程组改写为四元一阶常系数非齐次微分方程组, 以符合 Runge- Kutta的要求:
在这个方程中， x₁(t),x₂(t),x₃(t),x₄(t)分别为浮子的位移、振子的位移、浮子的速度、振子的速度。
代入第一问的相关参数与初值(2-6), 使用 Runge- Kutta 法迭代 40 个周期，取时间间隔为 tᵢ₊₁-tᵢ=h=0.2s,
计算得到的x₁(t)与 x₂t-x₁t,x₃t
与
x₄(t)如图5所示:
| 时间(s) | 浮子 | 振子 | ||
| 位移(m) | 速度(m/s) | 位移(m) | 速度(m/s) | |
| 10 | -0.205920238 | -0.645809997 | -0.234964788 | -0.691513475 |
| 20 | -0.61064358 | -0.246366756 | -0.660770395 | -0.267815187 |
| 40 | 0.266694073 | 0.30147726 | 0.277778099 | 0.320606667 |
| 60 100 | -0.328078998 -0.089379868 | -0.485631323 -0.606618652 | -0.350656448 -0.095246409 | -0.517789228 -0.646739304 |
5.3问题二模型的建立
问题二分为两种情况，第一种情况阻尼系数恒定，微分方程组是线性的，因此可以使用力电比拟的方法，不去求解含参的微分方程，而是将整个振荡浮子式波浪发电系统等效为一个电感-电容-电阻的二阶电路，即将问题转化为求解可变线性电路元件的最大功率的问题，通过函数极值的方法求解析解。
第二种情况，阻尼系数与浮子和振子的相对速度的绝对值的幂成正比，无法采用力电比拟法，故建立最优化模型进行求数值解。
5.3.1力电比拟法模型的建立与求解
力学系统和电学系统在一些情况下是类似的，如果两个系统都是二阶齐次常系数微分方程，那么两个系统具有对偶关系[2]：
mx+ px+ kx=0 (3-1)
Lā+ Rq+1/ Cq=0 (3-2)
其中力学方程遵循牛顿第二定律，电学方程遵循无源RLC 回路压降等于0的方程。质量对应电感，阻尼系数对应电阻，弹簧刚度对应电感的倒数，位移对应电荷，不难发现两个方程是等价的。因此，可以用电学系统的语言描述力学系统。
对于情况一，考虑阻尼系数恒定时的方程组：
11
{m1+mcx1+k1x1+k2x1-x2+k3x1-x2+k4x1-fcosωt=0}3-3
…………………………………
两个力学方程，对应两条电学回路。电路中各变量与机械系统中各物理量的对应关系如表1所示。
表1 振荡浮子式波浪发电系统机械量与电量等效对照图
| 机械物理量 | 电路变量 | ||
| 名称 | 符号 | 名称 | 符号 |
| 浮子质量与垂荡附加质量之和 | m₁+m。 | 电感1 | L₂ |
| 振子质量 | m₂ | 电感2 | L₂ |
| 浮子位移 | x₁ | 回路1电荷量 | q₁ |
| 振子位移 | x₂ | 回路2电荷量 | 92 |
| 浮子速度 | x₁ | 回路1电流 | l₂ |
| 振子速度 | x₂ | 回路2电流 | l₂ |
| 垂荡静水恢复力系数 | k₁ | 电容1⁻² | C₁¹ |
| 弹簧刚度 | k₂ | 电容2⁻² | C₂¹ |
| 垂荡兴波阻尼系数 | k₄ | 电阻1 | R₁ |
| 阻尼器阻尼力系数 | k₃ | 电阻2 | R₂ |
| 垂荡激励力振幅 入射波浪频率 | f ω | 余弦交流电源振幅 余弦交流电源频率 | A ω |
由表1转化为电路方程组可得：
(L1dl1dt+1C1∫I1dt+1C2∫I1-Izdt+R2I1-I2+R1I1=Acosωt
(3-4)
这是根据回路电流法列出的电路方程组，两个方程对应两条支路，可根据方程组画出电路图，如图6所示。
这是一个二阶 RLC 电路，应用相量法求解[2]。将L₁，C₁，R₁三个器件看作一个整体，其阻抗记作Z₁，根据三者串并联关系，用Z₁₈和Z₁ₓ表示Z₁的实部和虚部，可得：
Z1=jωL1+1jωC1+R1=Z18+jZ1x
(3-5)
Z1R=R1
Z₁ₓ=ωL₁-1/ωC₁
将 L₂,c₂,R₂
三个器件看作另一个整体，其阻抗记作Z₂，根据三者串并联关系，用Z₂₈和Z₂x表示Z₂的实部和虚部，可得：
Z2=jωL21ac2+Rzjωb2+1jωC2+R2=Z2R+jZ2x
(3-6)
Z2R=C2l2R2ω2-C2L2R2ω21-C2l2ω2c22Rz2ω2+1-C2l2ω22
Z2=C22L2R22ω3+12ω1-c2l2ω2C22R2ω2+1-c2l2ω22
在力电比拟模型中，阻尼器阻尼力系数k₃被比拟成电阻R₂，振子与浮子的弹簧被比拟成C₂，故PTO的输出功率为C₂与R₂的功率[2]。在L₂，C₂，R₂三个器件组成的负载中，L₂与C₂在一个周期内充电和放电的功率相等，平均功率为0[1]，而R₂在一个周期中始终对外做功，故 L₂,C₂,R₂三个器件的平均功率即为R₂的平均功率[1]:
P2=I2R2=A22Z28Z18+Z2R2+Z1x+Z2x2
(3-7)
其中Z₂₈与Z₂x都是R₂的函数，因此P₂也是R₂的函数，求其极值点即可求得R₂的最大功率和此时k₃的值，进而求得PTO的最大输出功率。
由问题二的力学参数(3-8)可以得到电路参数(3-9)：
曼● (3-8)
l1=02i-24≥93c1=31c2+168cC21=168c+10>5R1=16π2+395≥7a=484n≤481
q10=0q20=0710=0I20=0
(3-9)
13
将(3-9)相关参数代入(3-5)(3-6)(3-7)中, 可以算得此时:
Z₁₈=16784
Z₁ₓ=-89495
|
|
|
(3-10) |
|
Z2x=4584+4127×10-6R220.724+7661×10-10R22 |
P2=1.717×107×9450×108R2+R231.307×1018R2+2328×109R22+4820R23+R24(3-12
令 dP2dR2=0,
求得P₂的极值点为 R₂=37193812,
此时 P₂ₘₐₓ=2293341V.
为了验证结果合理性，绘制出P₂关于R₂的函数图像，如图8所示。
从图7中可以看出，P₂的极值点为( R₂=37193812,
此时 P₂ₘₐₓ=229334W,
结果合理。故当 k3=R2=3719381N⋯⋅m-1
时，阻尼器平均功率最大为229.334W。力电比拟法求最大功率不需要用到数值方法，求出的结果是精确解。
5.3.2最优阻尼系数模型的建立与求解
设发电系统一个周期内的平均功率为P。已知位移x₁和x₂是时间的函数，即：
x1=x1tx2=x2t
(3-13)
因此在一个周期内阻尼器的功率可以表示为：
P=c0t0+7 F3dx1-x27
(3-14)
其中F₃为阻尼器的力的大小。将F₃代入，可以得到：
P=t0t0+7 k3x1-x22dtT
(3-15)
由第一问的数据可知，在时间较小的地方浮子与振子的速度幅度尚不稳定，为了
求得稳定的平均功率，减小误差，故这里计算100个周期的浮子与振子的速度，并取最后十个周期的平均功率作为目标函数：
P=∞νr10or k3x1-x22dt107
(3-16)
取α与β作为决策变量，且 α∈010000,β∈01,
k3=α|x1-x2|β
(3-17)
第一个约束条件为振子不能碰到浮子顶端，振子高度为H₃，可得：
x₀+x₂-x₁+M₃≤H₁
(3-18)
第二个约束条件为浮子不能飞出水面，而问题一中求得浮子在平衡位置时的吃水深度为2.8m, 可得:
x₁≤2,8
(3-19)
第三个约束条件为弹簧长度大于零，可得：
x₁-x₂+x₀≤1₀
(3-20)
最终建立如下最优阻尼系数模型：
αβ=argarαβmaxP=00r10cosr k1x1-x22dt10T
(3-21)
|
(m₁+m₀)x₁=-k₁x₁-k₂(x₁-x₂)-k₃(x₁-x₂)-k₄x₁+fcosωt m₂x₂=k₂(x₁-x₂)+k₃(x₁-x₂) |
|
k₃=α|x₁-x₂|β |
|
s. t.( x₀+x₂-x₁+H₃≤H₁ |
|
x1≤2.8 |
L₀
α∈[0,100000]
β∈[0,1]
由于目标函数为定积分，但是被积函数只有离散点处的值，因此只能使用数值积分计算该定积分。梯形法通过将一个区域分为包含多个更容易计算的梯形区域，对区间内的积分计算近似值.对于具有N+1个均匀分布的点的积分，近似值为：
ab fxdx=b-a2N∑n=1Nfxn+fxn+1
=b-a2Nfx1+2fx2+⋯+2fxN+fxN+1
(3-22)
由于只有两个决策变量，求解此优化问题较简单，使用 MATLAB 的数值积分函数计算目标函数，使用优化工具箱编程求解，得到最优解与决策变量的值：
Pₘₐₓ=228.89W
α=0.3975α=100000
为了验证结果的可靠性，以0.001作为α的步长，1000作为β的步长，且β∈
15
[0,200000], 对P₂进行遍历求值, 结果如图8所示:
可以看出，图像存在一个明显的“山脊”，高度在228W左右，并且山脊的方向沿着β的正方向，并向α的正方向偏：在β∈[0，100000]时，极值点的确在β的边界处取得：如果β的上界能进一步提高，极值点的β还会进一步增加，对应的 Pmax也能跟着提高。
5.4问题三模型的建立
由题目可知，该波浪能装置的系统除了垂荡运动外，还作纵摇运动，因此需要一组新的变量来描述其运动。将系统视为一个质量不变，但是质心位置变化的刚体：根据刚体运动的规律，可将系统的运动分为平动和转动，分别可以对应垂荡和纵摇运动。
假设浮子和振子进行纵摇运动的角位移很小，因而在垂荡运动时浮子和振子近似处在同一条直线上，纵摇运动不会影响到竖直方向上的受力，垂荡运动仍可以沿用第一问的建立的模型求解。
5.4.1质心位置和各点相对于转轴转动惯量的计算
浮子由一个圆锥壳体和一个壳体组成，它们的总面积为：
S=12⋅2πRR2+H22+3⋅πR2+2πR⋅H1=9.64m2
根据题设，浮子与振子质量均匀，则浮子处处具有相等的面密度：
σ=m15
(4-1)
则在第一问的坐标系下，浮子的质心：
xo1=1m1-28-2 σ⋅2π⋅54⋅xR2+H22H2dx+-21 σ⋅πR2⋅dx=-1.236m;
振子的质心
代码部分
支撑材料列表：
result1-1. xlsx
result1-2. xlsx
result3. xlsx
module1 1. m 用于给出第一问第一种情况的微分方程
modulel 2. m 用于给出第一问第二种情况的微分方程
module2 1. m 用于给出第二间第一种情况的微分方程
module2 2. m 用于给出第二问第二种情况的微分方程
module3 1. m 用于给出第三间的微分方程
module4 1. m 用于给出第四问的微分方程
v1 jiexi. m用于给出第一问第一种情况浮子速度的解析方程
v2 jiexi. m用于给出第一问第一种情况振子速度的解析方程
x1 jiexi. m用于给出第一问第一种情况浮子位移的解析方程
x2 jiexi. m用于给出第一问第一种情况振子位移的解析方程
laplace. nb 计算拉普拉斯逆变换
circuit2 1. nb计算力电比拟电路的数值解
circuit2 2. nb计算力电比拟电路的符号解
CreateProblem2 1. m 给第二问第一种情况建立优化问题模型
CreateProblem2 2. m 给第二问第二种情况建立优化问题模型
CreateProblem4 1. m 给第四问建立优化模型
task1. ml x 第 一问主函数
task2. ml x 第二问主函数
task3. ml x 第三问主函数
task4. ml x 第四问主函数
error estimate. ml x 误差检验主函数
modulel 1. m 用于给出第一问第一种情况的微分方程
function dxdt =module1 1(t,x)
m1=4866;
m2=2433;
mc=1335.535;
rho=1025;
g=9.81;
R=1;
H 1=3;
k1=31557.298;
k2=80000;
k3=10000;
k4=656.3616;
x0=m2*g/k2;
fmax=6250;
omega=1.4005;
x1=x(1);
x2=x(2);
x3=x(3);%dx1dt
x4=x(4);%dx2dt
dxdt=[x3;
x4;
k1*×1*x1=-1)-k1*-1*(x1<-1)-k2*x1-x2+0-k3*x3-x4
(- -
k4⁺x3+fox⁺cosonega⁺)/m1+nc;
-g+k2*x1-x2+x0+k3*x3-x4/n2;
];
end
modulel 2. m 用于给出第一问第二种情况的微分方程
function dxdt =module1 2(t,x)
m1=4866;
m2=2433;
mc=1335.535;
rho=1025;
g=9.81;
R=1;
H 1=3;
x0=0.2980425;
k1= pi*g*R*R* rho;
k2=80000;
k4=656.3616;
fmax=6250;
omega=1.4005;
x1=x(1);
x2=x(2);
x3=x(3);% dxidt
x4=x(4);%dx2dt
k3=10000* sqrt( abs(x3-x4));
dxdt=[x3;
x4;
(-k1*x1*(x1>=-1)-k1*(-1)*(x1<-1)-k2*(x1-x2)-k3*(x3-x4)-
k4*x3+ fmax* cos( omega*t))/(m1+ mc);
-g+k2*x1-x2+x0+k3*x3-x4/m2;
];
end
module2 1. m 用于给出第二问第一种情况的微分方程
function dxdt =module2 1(t,x,k3)
m1=4866;
m2=2433;
mc=1165.992;
rho=1025;
g=9.81;
R=1;
H 1=3;
x0=0.2980425;
k1= pi*g*R*R* rho;
k2=80000;
k4=167.8395;
fmax=4890;
omega=2.2143;
x1=x(1);
x2=x(2);
x3=x(3);%dx1dt
x4=x(4);%dx2dt
dxdt=[x3;
x4;
(-k1*×1-k2*x1-x2-k3*x3-x4.
k4*x3+fmax*cosonega+c/m1+nc;
-g+k2*x1-x2+x0+k3*x3-x4/m2;
];
end
module2 2. m 用于给出第二问第二种情况的微分方程
function dxdt =module2 2(t,x, beta, alpha)
m1=4866;
m2=2433;
mc=1165.992;
rho=1025;
g=9.81;
R=1;
H 1=3;
x0=0.2980425;
k1= pi*g*R*R* rho;
k2=80000;
k4=167.8395;
fmax=4890;
omega=2.2143;
x1=x(1);
x2=x(2);
x3=x(3);% dxidt
x4=x(4);%dx2dt
k3= beta*( abs(x3-x4))^( alpha);
dxdt=[x3;
x4;
(-k1⁴+x1-k2*x1-x2-k3.°x3-x4.
k4*x3+ fmax* cos( omega*t))/(m1+ mc);
-g+(k2*(x1-x2+x0)+k3.*(x3-x4))/m2;
];
end
module3 1. m 用于给出第三问的微分方程
function dxdt =module3 1(t,x)
omega=1.7152;
m1=4866;
m2=2433;
mc=1028.876;
rho=1025;
g=9.81;
R=1;
H 1=3;
k1=31557.298;
k2=80000;
k3=10000;
k4=683.4558;
fmax=4890;
x0=m2*g/k2;
%J1=1040.6332;
%J1=1420.6332;
Jc=7001.914;
L1=8890.7;
L2=250000;
L3=1000;
L4=654.3383;
Mmax=1690;
x1=x(1);
x2=x(2);
x3=x(3);% dxidt
x4=x(4);%dx2dt
theta1=x(5);
theta2=x(6);
theta3=x(7);%dtheta1dt
theta4=x(8);%dtheta2dt
J2=2433*0.5/3*(3*(x2-x1+x0)^2+3*0.5*(x2-x1+x0)+0.5*0.5);
J1=0.5*4866*(1.236+x2-x1+x0+0.25).^2/9;
001=1/3*(0.25+x2-x1+x0+1.236);
扫一扫
2006
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
