荐 2019全国大学生数学建模竞赛题目 A题 高压油管的压力控制

A题 高油管的压力控制

问题1

问题1的第一小问，如何设置单向阀每次开启的时长，使得高压油管内的压力尽可能稳定在$100MPa$左右。首先要明确高压油管的工作原理以及过程。燃油经过高压油泵从A处进入高压油管，再由喷口B喷出。燃油进入和喷出的间歇性工作过程会导致高压油管内压力的变化。

高压油管内部压力变化

进气状态

高压油管内的初始压力为$100MPa$，当高压油泵从A出进入高压油管时，此时高压油管的压力是增大的。质量守恒定律，高压油管内的燃油的体积不变，根据密度-质量-体积公式
$\mathrm{ \rho}=\frac{m}{V}$
质量增大，密度增大，此时高压油管内的压强是增大的。这是在起初一个短时间内的状态。在开启了一段时长后（这段时长的参数是自己设置的，在后面要用到遍历及优化的思想，本题按$0.01ms$遍历搜索最佳时长为$0.29ms$） t₁(0.287ms) 后，要关闭 10ms 。在进气阶段两个蓝色标注是非常重要的量，我们要记在心里，贯穿整个高压油管工作的始尾。下图取$500ms$为界，此时是高压油管只进气压强的变化图。

喷气状态

喷气状态题中已经给出状态图

每100ms除了喷油嘴不工作的时间，其他时间都是在这个状态下进行的，喷油嘴喷气，质量守恒定律，高压油管内的燃油的体积不变，根据密度-质量-体积公式
$\mathrm{ \rho}=\frac{m}{V}$
质量减小，密度减小，在这一时间段高压油管内的压强是减小的。这里比较重要的参数是$2.4ms$喷气开始的时间 t₀(52ms) ，下图取$500ms$为界，此时是高压油管只喷气压强的变化图。

进喷气状态结合

进气状态和喷气状态是一个相互的过程，在高压油管工作时，进喷气会存在同时进行的状态。这时只需更改喷气开始的时间 t₀ 与进气时长 t₁ ，便可得到想要的结果，这里取2000ms的时长。

压强与密度关系

首先要进行迭代计算，通过题目中所给压强与密度关系式递推出每一密度下的压强。这里实际用Excel表格就能操作附件3-弹性模量与压力，第三列第五行单元格公式：$=(B5*0.85)/(B5-A5+100)$，下拉即可。

通过python pandas库导入文件的两列数据，用polyfit函数进行数据拟合，这里比较简便，可以选择二次项或者三次项的系数，下面完整代码有体现到这一代码的实现。这里采用了二次拟合，效果不是特别的好，也可选择其他更良好的拟合方式，这样后续的数据更准确。这里只用到了密度转压强的函数，没有设计压强转密度的函数。

关于进喷气函数的实现

整个过程是在$t=0.01ms$的逐步递推进行的，首先要想到$t=0$的出状态，自己可以在纸上写一下，经过$t=0.01ms$高压油管的状态是什么，经过$t=0.01ms+0.01ms=0.02ms$高压油管的状态是什么，首先计算单位时间进油量，乘$160MPa$密度即为进入的质量，质量变化、体积不变、密度变化、压强变化，每走$0.01ms$更新一下高压油管内压强的状态，即可完成整个过程.

关于遍历及优化

这个考虑到高压油管内部压强比较稳定，所以进喷气的质量大致相同，取质量相同求得$t0=0.28ms$，在附近进行遍历即可。
优化这个可以自己设计，本题要求在$100MPa$较稳定的工作状态，我们可以计算整个函数与$P=100MPa$的面积，面积越小时优化结果越好。但是在实际代码设计优化状态下可能会出现问题，这个我们还要自己把握。我的想法每遍历一次输出一张图片到文件夹里，然后打开文件夹翻阅查看，再通过那段区间设置优化的限制。

稳定度分析（0.27ms,0.28ms,0.29ms）

第三张图片吻合的明显要前两张图片好，此时是$t0.29ms$的情况，$0.001ms$搜索会找到更稳定的工作状态。

平衡状态代码

代码比较灵活，里面参数可以自行修改，在使用时要导入python第三方库。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

data = '附件3-弹性模量与压力.xlsx'
data = pd.read_excel('附件3-弹性模量与压力.xlsx')
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

y = data[u'压力(MPa)']
x = data[u'燃油密度(mm3/ms)']

plt.plot(x, y, marker='.', markersize=0.1, label='real')
f = np.polyfit(x, y, 2)

k = f[0] * (x ** 2) + f[1] * x ** 1 + f[2]
plt.plot(x, k, label='fit')
print(f)

plt.ylabel('燃油密度(mm3/ms)')
plt.xlabel('压力(MPa)')
plt.title('压力与燃油密度的关系')
plt.legend()
plt.show()

C = 0.85  # 流量系数
A = math.pi * 0.7 ** 2  # 小孔的面积 1.539380374mm2
P2 = 160  # 高压油管端压力
P1 = 100  # 100MPa压力
rou160 = 0.868705986  # 高压油管端油气密度
rou = 0.85  # 100MPa油气密度

v0 = 500 * math.pi * 5 ** 2  # 高压油管体积
m = v0 * 0.85  # 高压油管油气质量
t = 0  # 开始时刻 t每增进1则进行0.01ms
P_list = []  # 压强时刻表


# 密度转压强
def P(n):
    return f[0] * (n ** 2) + f[1] * n ** 1 + f[2] - 1.6


# 进油函数 t1为进油时间段
def enter(t, t1):
    if 0 < t % (t1 + 10) < t1:
        global m
        global P1
        global rou
        Q1 = C * A * math.sqrt(2 * (P2 - P1) / rou160)  # 单位时间进油量
        det_m1 = Q1 * rou160 * 0.01  # 0.01为步长 质量改变量
        m = m + det_m1  # 更新质量
        rou = m / v0  # 更新密度
        P1 = P(rou)  # 更新压强
    else:
        pass


# 出油函数 t1为出油时间点
def out(t, t2):
    global rou
    global m
    global P1
    if t % 100 == t2:
        pass
    if t2 < (t % 100) < t2 + 0.2:
        Q2 = (t - t2) % 100 * 100
        det_m2 = Q2 * rou * 0.01  # 0.01为步长 质量改变量
        m = m - det_m2  # 更新质量
        rou = m / v0  # 更新密度
        P1 = P(rou)  # 更新压强
    elif t2 + 0.2 <= (t % 100) < t2 + 2.2:
        Q2 = 20
        det_m2 = Q2 * rou * 0.01  # 0.01为步长 质量改变量
        m = m - det_m2  # 更新质量
        rou = m / v0  # 更新密度
        P1 = P(rou)  # 更新压强
    elif t2 + 2.2 <= (t % 100) <= t2 + 2.4:
        Q2 = ((t - t2) * (-100) + 240) % 100
        det_m2 = Q2 * rou * 0.01  # 0.01为步长 质量改变量
        m = m - det_m2  # 更新质量
        rou = m / v0  # 更新密度
        P1 = P(rou)  # 更新压强
    else:
        pass


sum_min = 10000000000000  # 记录压力波动偏差最小值
P_real = []  # 记录压力最小时刻表
t_1 = 0  # 记录放气时间段
t_2 = 0  # 记录进气时间点

for i in np.arange(0.287, 0.288, 0.001):  # 遍历开阀门时间段 start->end
    for j in np.arange(52, 53, 1):  # 遍历出气时刻 start->end
        sum = 0  # 记录压力波动偏差
        t = 0
        P_list = []
        while (t <= 2000):  # 时间可修改，单位为ms
            t = t + 0.01
            enter(t, i)
            out(t, j)
            P_list.append(P1)
            sum = sum + abs(P1 - 100)
        print('i =', i, 'j =', j)
        print(sum)
        if sum_min > sum:
            sum_min = sum
            t_1 = i
            t_2 = j
            P_real = P_list

while (t <= 2000):  # 时间可修改，单位为ms
    t = t + 0.01
    enter(t, t_1)
    out(t, t_2)
    P_list.append(P1)

print(t_1, t_2)
x1 = np.arange(0, len(P_real) / 100, 0.01)
x2 = np.arange(0, len(P_real) / 100, 1000)
y_ = np.linspace(100,100,len(x2))
plt.xlabel('t(ms)')
plt.ylabel('P(KPa)')
plt.plot(x1, P_real, label='压强-时间')
plt.plot(x2, y_,marker='.',markersize=0.1)
plt.legend()
plt.show()

增压模型（0.9ms,0.75ms,0.73ms时对应的增压过程）

问题2

问题2比问题1复杂了许多，进气过程不再是简单的双段过程，由高压油泵供油提供进气，针阀的开启会使高压油管油气喷出。这里就要了解高压油泵的进气方式，凸轮的工作原理其实了解了这个工作原理后续工作就可以继续进行了。

进气状态

凸轮转动角速度会影响到凸轮供气周期。凸轮每运动一周会带动活塞的上下移动，向上移动首先会使柱塞腔内的压力增大，当达到$100MPa$时，下一时刻柱塞腔压力大于高压油管内压力，气体被压入。至于压入气体的周期由角速度决定。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
data = '附件1-凸轮边缘曲线.xlsx'
data = pd.read_excel(data)

y = data[u'极径（mm）']
x = data[u'极角（rad）']

f = np.polyfit(x, y, 6)  # 六次函数拟合


def R(x):  # 极径 #x为极角 弧度制
    return f[0] * x ** 6 + f[1] * x ** 5 + f[2] * x ** 4 + f[3] * x ** 3 + f[4] * x ** 2 + f[5] * x + f[6]


w = 0.027  # 弧度制下的角速度
T = 2 * math.pi / w
y_max = []

for j in np.arange(117, T + 117, 1):  # 时间遍历，以周期为遍历end结束点
    ymax = 0  # 更新最大值
    for i in np.arange(0, 2 * math.pi, 0.01):  # 弧度制遍历
        y_point = R(i) * math.cos(w * j + i)  # 每一时间点下不同弧度坐标位置 y_point y坐标
        if y_point > ymax:
            ymax = y_point  # 每次遍历下来求得每一时间点下 不同弧度纵坐标最大值
    y_max.append(ymax)  # 将此时刻最大值添加到列表中

x = np.arange(0, len(y_max), 1)
plt.plot(x, y_max, label='活塞底部在y轴上的位置')
plt.xlabel('t/ms')
plt.ylabel('ymax/mm')
plt.legend()
plt.show()

x__ = np.arange(0, len(y_max), 1)  # 时间变量
vy = -w * R(0) * np.sin(2 * math.pi - w * x__)  # 速度变化

plt.plot(x__, vy, label='活塞速度')
plt.xlabel('t/ms')
plt.ylabel('v/mm')
plt.legend()
plt.show()

print(min(y_max), max(y_max), max(y_max) - min(y_max))  # 2.4157081016847566 7.247512164299576 4.83180406261482

喷气状态

针阀开启后，首先针阀下底面会与密封座产生一个圆环口，油气从这个口喷出，起初喷出油气时，圆环口面积一直增大，是完全喷出状态。当圆环口面积增大到喷孔面积时，这时候喷气产生了阻塞，每一时刻状态下喷气量不再增加。

进喷气状态结合

依然同第一问类似，只不过模型更为复杂，喷油嘴在附件二中给出，利用喷油嘴的流量变动推导高压油泵引起的流量变动，清楚各部分之间左右的关系。

关于凸轮的拟合

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
data = '附件1-凸轮边缘曲线.xlsx'
data = pd.read_excel(data)

y = data[u'极径（mm）']
x = data[u'极角（rad）']


plt.plot(x,y,marker='.', markersize=0.1,label='real')

f= np.polyfit(x,y,6)
y_=f[0]*x**6+f[1]*x**5+f[2]*x**4+f[3]*x**3+f[4]*x**2+f[5]*x+f[6]
print(f)

plt.plot(x,y_ ,label='fit')
plt.legend()
plt.show()



# 1.圆半径
def R(x):
    return f[0]*x**6+f[1]*x**5+f[2]*x**4+f[3]*x**3+f[4]*x**2+f[5]*x+f[6]
# 2.圆心坐标
a, b = (0, 0)
# 参数方程
theta = np.arange(0, 2*np.pi, 0.01)
x = a + R(theta) * np.cos(theta)
y = b + R(theta) * np.sin(theta)
#绘图
plt.figure(figsize=(6,6))

plt.xlabel('x(mm)')
plt.ylabel('y(mm)')
plt.plot(x, y)
plt.show()

问题3

问题3要求给出高压油泵和减压阀的控制方案，大家好好想想吧QAQ。

本文地址：https://blog.csdn.net/qq_44864262/article/details/107351538