©2020-2024 ioDraw All rights reserved
题目回顾与写在前面
* 首先，我们队在历经了千辛万苦之后，光荣得获得了 省三......
* 队伍构成 物理*2 + 计算机*1
* 队伍分工 计算机-->受力分析 物理-->数值计算
* 总评：图一乐，狠乐！物理系，计算机系嘛，不怎么看建模的啦！
* 如果只是考虑力学问题的话，我们分析得肯定还不太到位，但是这个题肯定很好分析
* 所以，这个题主要在数值计算上
思路
基本假设
1.海水是无粘及无旋的。
2.浮子在线性周期微幅波作用下会受到波浪激励力（矩）、附加惯性力（矩）、兴波阻尼力（矩）和静水恢复力（矩）。
3.忽略中轴、底座、隔层及 PTO的质量和各种摩擦。
4.初始浮子和振子平衡于静水中。
问题一的求解
问题一浮子振子整体的受力分析
问题一振子在浮子参考系下的受力分析
数值求解
4级显式Runge-Kutta方法
求解过程的代码：
import numpy as np import matplotlib.pyplot as plt from matplotlib.animation
import FuncAnimation import matplotlib matplotlib.rcParams["font.sans-serif"] =
["SimHei"] matplotlib.rcParams["axes.unicode_minus"] = False def
runge_kutta4(df, a, b, h, y0): num = len(y0) x = np.arange(a, b+h, h) w =
np.zeros( (x.size, num) ) w[0, :] = y0 for i in range(x.size - 1): s0 =
df(x[i], w[i, :],i*h) s1 = df(x[i] + h/2., w[i, :] + h * s0 / 2.,i*h) s2 =
df(x[i] + h/2., w[i, :] + h * s1 / 2.,i*h) s3 = df(x[i+1], w[i, :] + h *
s2,i*h) w[i+1,:] = w[i,:] + h * (s0 + 2*(s1+s2) + s3) / 6. return x, w def
df(x, variables,i): th1, th2, om1, om2 = variables A = np.zeros((2, 2)) b =
np.zeros(2) A[0, 0] = 2433+6201.535 A[0, 1] = 2433 A[1, 0] = 2433 A[1, 1] =
2433 b[0] =
6250*np.cos(1.4005*i)-10045*4*np.arctan(1)*th1*(th1<=0.999989731)-10045*4*np.arctan(1)*0.999989731*(th1>0.999989731)-656.3616*om1
b[1] = -80000*th2-10000*om2 dom1, dom2 = np.linalg.solve(A,b) return
np.array([om1,om2,dom1,dom2]) a, b = 0.0,180.0 h = 0.01 th10 = 0.1 th20 = 0.1
om10 = 0.1 om20 = 0.1 y0 = np.array([th10, th20, om10, om20]) # 计算求解 t, w =
runge_kutta4(df, a, b, h, y0) th1 = w[:, 0] th2 = w[:, 1] om1 = w[:, 2] om2 =
w[:, 3] plt.plot([h*i for i in range(len(th1))],th1,label="x1") plt.plot([h*i
for i in range(len(th2))],th2,label="x2") #plt.plot([i*h for i in
range(len(om1))],om1,label="x1'") #plt.plot([i*h for i in
range(len(om2))],om2,label="x2'") plt.xlabel("时间/s")
plt.ylabel("位移/m",rotation=True) #plt.ylabel("速度/m*s^-1",rotation=True)
plt.legend() plt.title("x1 x2 随时间的变化量；求解步长 -m1"+str(h)) plt.savefig("x1 x2
随时间的变化量 求解步长 -m1"+str(h)+".jpg") plt.show() #plt.title("x1' x2' 随时间的变化量；求解步长
"+str(h)) #plt.savefig("x1' x2' 随时间的变化量 求解步长 "+str(h)+".jpg") #plt.show()
import xlsxwriter as xls #th1 = np.array(th1[::20]) #th2 = np.array(th2[::20])
#om1 = np.array(om1[::20]) #om2 = np.array(om2[::20]) th1 = np.array(th1) th2 =
np.array(th2) om1 = np.array(om1) om2 = np.array(om2) workbook =
xls.Workbook("1-1-000-m1.xlsx") worksheet = workbook.add_worksheet("Sheet1")
headings = ["x1","x2","v1","v2"] worksheet.write_row("A1",headings)
worksheet.write_column("A2",th1) worksheet.write_column("B2",om1)
worksheet.write_column("C2",th1+th2) worksheet.write_column("D2",om1+om2)
workbook.close()
问题二
变步长搜索法
矩形法数值积分
问题二解法
60秒-80秒
80秒-100秒
100秒-120秒
最大平均功率/瓦特
290.141018
282.3656312
285.2766
289.2045
291.1627
291.2182
Beta值
38130
37850
37940
37970
37970
37970
Alpha值
0.098
0.1
0.1
0.1
0.1
0.1
问题三
数值解（类似一，二的解法）
物理量 \ 时间
10s
20s
40s
60s
100s
浮子垂荡位移
-0.685529938
-0.604544566
0.195435534
-0.156064651
0.142943636
浮子垂荡速度
0.5487404
-0.698097309
0.922311628
-0.867269801
-0.917710776
浮子纵摇角位移
0.0293679
0.001331413
-0.000590353
-0.002998229
-0.015015832
浮子纵摇角速度
-0.115251745
0.020961283
-0.03312735
0.042788285
0.051748796
振子垂荡位移
-0.764772476
-0.657517866
0.199876278
-0.159952226
0.168056577
振子垂荡速度
0.565213931
-0.788403868
1.017008287
-0.952503212
-0.997533033
振子纵摇角位移
0.030535835
0.001439153
-0.000618481
-0.003153401
-0.015074399
振子纵摇角速度
-0.117325818
0.019228743
-0.031151057
0.046614327
0.056193417
问题四
我们的优化模型：
小结
* 这道题受力分析的成分要更大一点哈
* 分析完了就数值计算
* 受力分析不完整，我们就不放出来了，以后再完善吧。