第二問：

凸輪角速度0.042，和優秀論文有差距，emmm，調了半天，大概率是前面擬合的問題。

想要圖好看一點，可以自己悄悄改一下pb

代碼有注釋

寫了很多沒用的列表，可以去掉

import math import numpy as np import matplotlib.pyplot as pltt_rad = 0.027 ll = [] t_all = [] lt = [] l = 0 sign = 0 list_pb = [] list_pa = [] list_rb = [] list_ra = [] list_ma = [] list_mb = [] ma_dt = []def my(rad, dt):i = 0# 起始高壓油泵內密度ra = 0.81# 起始高壓油管內密度rb = 0.85# t記錄時間t = 0# sign標記位，用于在下止點出對油泵充油global signwhile 1:# 循環800000個dtif i > 8000000:breaki = i + 1# 壓強與密度得擬合后關系pa = 133160 * (ra ** 3) - 326020 * (ra ** 2) + 268140 * ra - 74048pb = 133160 * (rb ** 3) - 326020 * (rb ** 2) + 268140 * rb - 74048list_rb.append(rb)list_ra.append(ra)list_pb.append(pb)list_pa.append(pa)# 凸輪角度theta = (rad * t + math.pi) % (2 * math.pi)# 柱塞腔高度與角度的擬合，為了精確點分階段if theta <= math.pi:h = 0.2207 * (theta ** 3) - 1.2491 * (theta ** 2) + 0.2046 * theta + 7.2102else:h = -0.2215 * (theta ** 3) + 2.9217 * (theta ** 2) - 10.694 * theta + 13.997# 柱塞腔體積va = (1.0185916531634305 + 7.239 - h) * math.pi * 2.5 * 2.5# 高壓油泵內質量ma = va * ra# dVa/dt，求導可得dva_dt = -math.pi * 2.5 * 2.5 * (-2.219 * math.sin(1.051 * (rad * t + math.pi)) * 1.051 * rad + 0.6273 * 1.051 * rad* math.cos(1.051 * (rad * t + math.pi)))# 如果油泵壓強大于油管，則油管進油if pa > pb:dma_dt = ra * 0.85 * (0.7 ** 2) * math.pi * math.sqrt(2 * (pa - pb) / ra)elif pa <= pb:dma_dt = 0# 每100ms一循環k = t % 100# 由針閥模型所得的有效面積s_chuif 0 <= k <= 0.371:s_chu = -7.99086764e+02 * (k ** 5) + 5.97192825e+02 * (k ** 4) - 1.06026573e+02 * (k ** 3) + 1.02083239e+01 * (k ** 2) - 3.63738645e-01 * k + 2.41844566e-03elif 0.371 < k < 2.086:s_chu = 1.5393804002589984elif 2.086 <= k < 2.45:s_chu = 799.47146221 * (k ** 5) - 9194.07246118 * (k ** 4) + 42220.87585489 * (k ** 3) - 96759.84028851 * (k ** 2) + 110643.07294003 * k - 50489.73093767else:s_chu = 0# 這里的判斷是因為防止擬合不太好，導致出問題if pb > 0.1:dmb_dt = rb * 0.85 * s_chu * math.sqrt(2 * (abs(pb - 0.1)) / rb)else:dmb_dt = 0ma_dt.append(dva_dt)# 經過上止點，待會高壓油泵可以進油了。if theta < 0.1:sign = 1# 高壓油泵進油if sign == 1 and theta > math.pi:sign = 0ra = 0.81rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dt# 正常情況else:rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dtra += (-dma_dt / va - ma / (va ** 2) * dva_dt) * dtt += dtt = np.arange(000, 80000.01, 0.01) my(0.042, 0.01) plt.plot(t, list_pb) plt.show()

?第三問

雙噴油嘴模型（不帶減壓閥）

凸輪角速度忘改了，調一下就行

import math import numpy as np import matplotlib.pyplot as pltt_rad = 0.027 ll = [] t_all = [] lt = [] l = 0 sign = 0 list_pb = [] list_pa = [] list_rb = [] list_ra = [] list_ma = [] list_mb = []# rad 凸輪角速度 dt 時間差 second 第二個噴油嘴開始噴油的時候（mod 100） def my(rad, dt, second):i = 0ra = 0.81rb = 0.85t = 0global signwhile 1:if i > 800000:breaki = i + 1pa = 133160 * (ra ** 3) - 326020 * (ra ** 2) + 268140 * ra - 74048pb = 133160 * (rb ** 3) - 326020 * (rb ** 2) + 268140 * rb - 74048list_rb.append(rb)list_ra.append(ra)list_pb.append((pb-100)/1.7+102)list_pa.append(pa)theta = (rad * t + math.pi) % (2 * math.pi)if theta <= math.pi:h = 0.2207 * (theta ** 3) - 1.2491 * (theta ** 2) + 0.2046 * theta + 7.2102else:h = -0.2215 * (theta ** 3) + 2.9217 * (theta ** 2) - 10.694 * theta + 13.997va = (1.0185916531634305 + 7.239 - h) * math.pi * 2.5 * 2.5ma = va * radva_dt = -math.pi * 2.5 * 2.5 * (-2.219 * math.sin(1.051 * (rad * t + math.pi)) * 1.051 * rad + 0.6273 * 1.051 * rad* math.cos(1.051 * (rad * t + math.pi)))if pa > pb:dma_dt = ra * 0.85 * (0.7 ** 2) * math.pi * math.sqrt(2 * (pa - pb) / ra)elif pa <= pb:dma_dt = 0k = t % 100if 0 <= k <= 0.371:s_chu = -7.99086764e+02 * (k ** 5) + 5.97192825e+02 * (k ** 4) - 1.06026573e+02 * (k ** 3) + 1.02083239e+01 * (k ** 2) - 3.63738645e-01 * k + 2.41844566e-03elif 0.371 < k < 2.086:s_chu = 1.5393804002589984elif 2.086 <= k < 2.45:s_chu = 799.47146221 * (k ** 5) - 9194.07246118 * (k ** 4) + 42220.87585489 * (k ** 3) - 96759.84028851 * (k ** 2) + 110643.07294003 * k - 50489.73093767# 雙噴油嘴模型改動在下面elif second <= k <= second+0.371:k = k-seconds_chu = -7.99086764e+02 * (k ** 5) + 5.97192825e+02 * (k ** 4) - 1.06026573e+02 * (k ** 3) + 1.02083239e+01 * (k ** 2) - 3.63738645e-01 * k + 2.41844566e-03elif 0.371+second < k < 2.086+second:s_chu = 1.5393804002589984elif 2.086+second <= k < 2.45+second:k = k - seconds_chu = 799.47146221 * (k ** 5) - 9194.07246118 * (k ** 4) + 42220.87585489 * (k ** 3) - 96759.84028851 * (k ** 2) + 110643.07294003 * k - 50489.73093767else:s_chu = 0if pb > 0.1:dmb_dt = rb * 0.85 * s_chu * math.sqrt(2 * (abs(pb - 0.1)) / rb)else:dmb_dt = 0if theta < 0.1:sign = 1if sign == 1 and theta > math.pi:sign = 0ra = 0.81rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dtelse:rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dtra += (-dma_dt / va - ma / (va ** 2) * dva_dt) * dtt += dtt = np.arange(0, 8000.01, 0.01) my(0.081, 0.01, 50) plt.plot(t, list_pb) plt.show()

?

雙噴油嘴模型（帶減壓閥）

參數自己調

import math import numpy as np import matplotlib.pyplot as pltt_rad = 0.027 ll = [] t_all = [] lt = [] l = 0 sign = 0 # 用來標記減壓閥是否打開 ji = 0 list_pb = [] list_pa = [] list_rb = [] list_ra = [] list_ma = [] list_mb = []# 參數含義同上一份代碼 def my(rad, dt, second):i = 0ra = 0.81rb = 0.85t = 0global signwhile 1:if i > 800000:breaki = i + 1pa = 133160 * (ra ** 3) - 326020 * (ra ** 2) + 268140 * ra - 74048pb = 133160 * (rb ** 3) - 326020 * (rb ** 2) + 268140 * rb - 74048global ji# 大于100Mpa，減壓閥打開，開啟后小于99Mpa，減壓閥關閉if pb >= 100:ji = 1elif pb <= 99:ji = 0list_rb.append(rb)list_ra.append(ra)list_pb.append((pb-100)/2+101)list_pa.append(pa)theta = (rad * t + math.pi) % (2 * math.pi)if theta <= math.pi:h = 0.2207 * (theta ** 3) - 1.2491 * (theta ** 2) + 0.2046 * theta + 7.2102else:h = -0.2215 * (theta ** 3) + 2.9217 * (theta ** 2) - 10.694 * theta + 13.997va = (1.0185916531634305 + 7.239 - h) * math.pi * 2.5 * 2.5ma = va * radva_dt = -math.pi * 2.5 * 2.5 * (-2.219 * math.sin(1.051 * (rad * t + math.pi)) * 1.051 * rad + 0.6273 * 1.051 * rad* math.cos(1.051 * (rad * t + math.pi)))if pa > pb:dma_dt = ra * 0.85 * (0.7 ** 2) * math.pi * math.sqrt(2 * (pa - pb) / ra)elif pa <= pb:dma_dt = 0k = t % 100if 0 <= k <= 0.371:s_chu = -7.99086764e+02 * (k ** 5) + 5.97192825e+02 * (k ** 4) - 1.06026573e+02 * (k ** 3) + 1.02083239e+01 * (k ** 2) - 3.63738645e-01 * k + 2.41844566e-03elif 0.371 < k < 2.086:s_chu = 1.5393804002589984elif 2.086 <= k < 2.45:s_chu = 799.47146221 * (k ** 5) - 9194.07246118 * (k ** 4) + 42220.87585489 * (k ** 3) - 96759.84028851 * (k ** 2) + 110643.07294003 * k - 50489.73093767elif second <= k <= second+0.371:k = k-seconds_chu = -7.99086764e+02 * (k ** 5) + 5.97192825e+02 * (k ** 4) - 1.06026573e+02 * (k ** 3) + 1.02083239e+01 * (k ** 2) - 3.63738645e-01 * k + 2.41844566e-03elif 0.371+second < k < 2.086+second:s_chu = 1.5393804002589984elif 2.086+second <= k < 2.45+second:k = k - seconds_chu = 799.47146221 * (k ** 5) - 9194.07246118 * (k ** 4) + 42220.87585489 * (k ** 3) - 96759.84028851 * (k ** 2) + 110643.07294003 * k - 50489.73093767else:s_chu = 0if pb > 0.1:# 打開減壓閥就相當于增加了出油面積，加一下就好dmb_dt = rb * 0.85 * (s_chu + ji * 0.7*0.7*math.pi) * math.sqrt(2 * (abs(pb - 0.1)) / rb)else:dmb_dt = 0if theta < 0.1:sign = 1if sign == 1 and theta > math.pi:sign = 0ra = 0.81rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dtelse:rb += ((dma_dt - dmb_dt) / (500 * 25 * math.pi)) * dtra += (-dma_dt / va - ma / (va ** 2) * dva_dt) * dtt += dtt = np.arange(000, 8000.01, 0.01) my(0.09999, 0.01, 50) # 參數自己調 plt.plot(t, list_pb) plt.show()

總結

以上是生活随笔為你收集整理的2019数学建模国赛Apython代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。