精华内容
下载资源
问答
  • 为了研究黄土导热性能与其物理力学参数之间的关系,文章对宁东黄土进行了导热系数、标准贯入试验及含水率、密度等参数的测试,对测得的数据进行了拟合及相关性分析,讨论了不同参数导热系数的关系。试验结果表明,非...
  • 唾液皮质醇SOS之间的关系不显着,相关系数为-0.13,表明呈负相关,导致观察到唾液皮质醇水平较低的人倾向于具有较高的SOS。 唾液皮质醇唾液Ca浓度之间的关系不显着,相关系数为+0.15,表明存在正相关关系,并...
  • 目的:探讨癌症患者及其家人心理压力与医患关系的关系。 方法:将患者随机分为干预组和对照组,选择PDRQ-15,pcl-c,SAS和SDS量表作为评估指标,比较诊断为PTSD的患者的去甲肾上腺素和多巴胺水平癌症和没有PTSD的人...
  • 通过理论分析和公式推导,找到一种判断常压储层连通性的新方法,即压力系数与深度成反比关系。 本文验证了压力系统与反比例函数之间的关系。 同一压力系统的功能是唯一的,单调的,并且具有统一的渐近线和对称轴,...
  • 实验表明,对数公式能够很好地拟合压力恢复曲线,系数A表明煤屑在不同平衡瓦斯压力下解吸差异性,和煤层瓦斯压力存在指数关系。根据此规律,可在测得现场煤屑解吸瓦斯压力恢复曲线的系数A后,结合已确定煤层瓦斯压力...
  • 采用水压致裂测量地应力方法,获得了鄂尔多斯盆地东南缘26口煤层气井地应力分布,通过统计分析,建立了二叠系山西组2煤储层地应力煤层埋藏深度之间相关关系和模型,揭示了现今地应力分布规律及受控机制。研究结果...
  • 通过理论分析室内试验,提出了煤的水力冲击破碎压力与其物理力学参数之间的关系式,以及水力冲蚀速度无因次冲击压力之间的关系。结果表明:通过加权平均值法计算得到的临界破碎压力与实际值吻合良好;在3~8倍临界...
  • 根据气体压力敏感系数Cp定义,不同围压温度条件下气体压力敏感系数与气体压力呈现对数关系;由于有效应力与基质收缩效应2方面的影响,在初始阶段随围压及温度的升高煤渗透率对气体压力的响应程度增大,基于气体压力敏感...
  • 结果表明:随着瓦斯压力的增大,煤体电阻率均呈下降趋势,电阻率y瓦斯压力x符合y=a+bln(x+c)(a,b,c为常数)的关系,且相关系数较高,电阻率在试验前期下降幅度较大,后期下降幅度变小并趋于平稳;瓦斯对煤体电阻率的影响可...
  • 在三叠纪后由于温度升高导致原油裂解,开始快速增压,一直到晚白垩世,最大压力系数达1.70;此后构造抬升温度降低、气体脱溶,导致压力下降到现今弱超压状况。川南盘龙场和大窝顶气田压力演化比较一致,奥陶纪之前为...
  • 通过物理相似实验模拟了煤层综放开采过程中采场支承压力的动态变化过程,得到不同推进距下工作面煤壁前方形成的超前支承压力分布规律,拟合出最大支承压力集中系数与工作面推进距离的关系曲线。得到了采空区支承压力的...
  • 在相同矿山压力条件下,破碎岩石渗透系数与渗流水压呈非线性反比例关系。认为断层带破碎岩石中含粘土矿物较多时,岩石易吸水膨胀、软化,增加对地下水渗流阻碍作用;瞬间施加渗流水压使样品被压密,导致孔隙被...
  • 通过力学分析和Ansys仿真分析,推导出膜片结构最大挠度表达式及光纤光栅中心波长与压力之间线性关系,并得出变形和应力云图,证明了力学分析与仿真分析一致性;对其进行了温补性能测试,去除了温度干扰,提高了精度;...
  • 设计了一种新颖的基于金属弹性膜片的压力传感器,这种传感器将布拉格光栅粘贴在金属膜片的中央,通过光纤光栅解调仪监测布拉格波长的漂移,给出了光纤布拉格光栅的布拉格波长偏移量与压力的关系表达式。实验结果表明,该...
  • 相关力学标定实验结果表明,该传感器输出应力具有良好线性关系,应力灵敏度系数约为20~30pm/N,测试精度较高。相似模拟实验结果表明,受开采影响,工作面前方一点支承压力呈现逐渐增加并最终趋于稳定变化趋势,...
  • 多元回归分析是一种研究多个自变量X一个因变量Y之间得关系得一种分析方式。因变量Y为连续类型。 X自变量通常是连续,也可以是...例如,在化学过程中,研究人员可以设置温度,压力和容器长度来搅拌一桶材料,...

    多元回归分析是一种研究多个自变量X与一个因变量Y之间得关系得一种分析方式。因变量Y为连续类型。 X自变量通常是连续的,也可以是离散的。简单回归分析,简单来说就是一个自变量X与一个因变量Y之间得关系。

    多元回归分析和简单回归分析都有两种情况:固定X取值和变化X取值。

    固定X取值的情况:各个X的级别由研究人员选择或由情况的性质决定。例如,在化学过程中,研究人员可以设置温度,压力和容器的长度来搅拌一桶材料,然后测量某种化学物质的浓度。执行回归分析来描述或解释自变量与因变量(某种化学物质的浓度)之间的关系,或预测因变量。

    变化X取值的情况:从总体中随机抽取个体样本,然后对每个个体测量X变量和Y变量。例如,可以从一个地区抽取一个成年人样本,并获取有关他们的年龄,收入和教育程度的信息,以便了解是否这些变量可以预测他们对空气污染的态度。回归和相关分析可以在变量X的情况下执行。可以从结果中获得描述性信息和预测性信息。

    多元回归是最常用的统计方法之一,并且了解其用法将有助于理解其他多元分析方法。

    拿一个数据集举例子,关于一个肺部功能的数据,这些数据适合变量X的情况。 身高用作变化的变量X以预测FEV1,计算得到方程式:

    a9c13190f6645d45a2813f36a6d9a9e2.png

    成人的FEV1会随着年龄的增长而降低。如果在多元回归方程中同时使用身高和年龄作为自变量,应该能够更好地进行预测。 年龄的斜率系数为负,高度的斜率系数为正。下图显示了FEV1分别与Age 和 Height的关系。

    6c64a7d43c304e5fce230991c0df255a.png
    只考虑一个变量时,FEV1分别与Age和Height的关系图。

    f653d6c4c70faef2dca4f16aa104d405.png
    FEV1与Age 和 Height的关系平面。

    构建这样的回归平面需要以下步骤:

    1.在X1,Y墙上画线,设置X2 = 0。

    2.在X2,Y墙上画线,设置X1 = 0。

    3.在步骤1和2中绘制的线处钉入“钉子”(也就是点)。

    4.用绳子连接成对的钉子并拧紧绳子。

    03aaafe7884fbbc300b955d934d2db8e.png
    该平面是Y在X1和X2上的回归平面。 给定X1和X2处Y的平均值是点X1,X2垂直上方平面上的点。

    得到的最后关于变量X与Y的函数:

    61a2a2bfa2b01266d1d0e0803b86d46f.png

    对比两个方程式,在单个预测变量方程中,身高系数为0.118。 当未考虑其他变量时,该值是FEV1随身高的变化率。 有两个预测因子时,身高系数为0.114,这被解释为FEV1随年龄调整后随身高变化的变化率。 调整年龄后,后一个斜率也称为FEV1在高度上的偏回归系数。 即使两个方程式均来自同一样本,简单回归系数和部分回归系数也可能不相等

    展开全文
  • 选用灰色关联度方法研究生活消费支出城镇居民污染物排放相关关系,结果显示家庭设备用品及服务、医疗保健支出多项污染物产生或排放量灰色关联系数较大。为了减少居民生活消费对区域生态环境的压力,建议提高...
  • 燃油进入和喷出的间歇性工作过程会导致高压油管内压力的变化。 压力密度关系 题目中给出燃油的压力变化量密度变化量成正比,比例系数为Eρ{\frac{E}{\rho} }ρE​ ,其中ρ{\rho}ρ为为燃油的密度。 Pk+1−Pk=E...

    A题 高压油管的压力控制

    问题一:

    稳定模型

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

    压力与密度关系

    题目中给出燃油的压力变化量与密度变化量成正比,比例系数为Eρ{\frac{E}{\rho} } ,其中ρ{\rho}为为燃油的密度。
    Pk+1Pk=E(Pk)ρk(ρk+1ρk),k=0,1,2,P_{k+1}-P_{k}=\frac{E\left(P_{k}\right)}{\rho_{k}}\left(\rho_{k+1}-\rho_{k}\right), k=0,1,2, \cdots
    首先要进行迭代计算,通过题目中所给压力与密度的出状态递推出燃油每一密度下的压力。这里实际用Excel表格就能操作:附件3-弹性模量与压力,在C1C1单元格填入(mm3/ms)燃油密度(mm3/ms)C202C202填入0.850.85,对应为100MPa100MPa压力下的密度。现在将C201C201填入公式:=(B201C202)/(B201A201+A202)=(B201*C202)/(B201-A201+A202),上拉操作;将C203C203填入公式:=(B203C202)/(B203A203+A202))=(B203*C202)/(B203-A203+A202)),下拉操作,拉满位置,此时这一列数据为对应压力下的密度。
    在这里插入图片描述

    通过pandaspandas库导入文件的两列数据,用polyfitpolyfit函数进行数据拟合,这里比较简便,可以选择二次项或者三次项的系数,下面完整代码有体现到这一代码的实现。这里采用了二次函数拟合,效果还算理想,也可选择其他更良好的拟合方式,不过比较麻烦,这样后续的数据会更准确。

    效果还可以哈。

    P(ρ)=10784.65349198ρ215693.86029255ρ+5648.09019803P(\rho)=10784.65349198\rho^{2}-15693.86029255\rho+5648.09019803
    ρ(P)=10784.65349198P215693.86029255P+5648.09019803\rho(P)=10784.65349198P^{2}-15693.86029255P+5648.09019803

    高压油管内部压力变化

    进气状态

    高压油管内的初始压力为100MPa100MPa,当高压油泵从A出进入高压油管时,高压油管内油气质量增大,体积不变,密度增大,此时高压油管内的压力增大。这是在起初一个较短时间内的状态。在开启了一段时长后(这段时长的参数是自己设置的,在后面要用到遍历及优化的思想,本题按0.001ms0.001ms遍历搜索最佳时长为0.288ms0.288ms t1(0.288ms) 后,要关闭 10ms 在进气阶段两个蓝色标注是非常重要的量,贯穿整个高压油管工作的始尾。下图取500ms500ms为界,此时是高压油管只进行进气状态压强的变化图。

    Rj(Tp)={0.85×π×0.722(160Pj)ρ(160)ρ(160)h 单向阀开启时 0 单向阀关闭时 R_{j}\left(T_{p}\right)=\left\{\begin{array}{ll} 0.85 \times \pi \times 0.7^{2} \sqrt{\frac{2\left(160-P_{j}\right)}{\rho(160)}} \rho(160) h& \text { 单向阀开启时 } \\ 0 & \text { 单向阀关闭时 } \end{array}\right.

    代码实现
    # 进油函数 参数t为时刻,t1为进油时间段
    def enter(t, t1):
        global m, P, ρ  # 全局变量,可做修改
        if 0 < t % (t1 + 10) < t1:
            Q = C * A * math.sqrt(2 * (160 - P) / ρ160)  # 单位时间进油量
            det_m = Q * ρ160 * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    进气过程展示

    喷气状态

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

    100ms100ms除了喷油嘴不工作的时间,其他时间都是在这个状态下进行的,喷油嘴喷气,高压油管内的燃油的体积不变,质量减小,密度减小,在这一时间段高压油管内的压力是减小的。这里比较重要的参数是时间段为2.4ms2.4ms开始喷气的时间 t0(50ms) ,下图取500ms500ms为界,此时是高压油管只进行喷气状态压强的变化图。

    Qj(Ts)={0Mod(t,100)Ts<0100(Mod(t,100)Ts),0Mod(t,100)Ts<0.2;20h0.2Mod(t,100)Ts<2.2;(240100(Mod(t,100)Ts))h2.2Mod(t,100)Ts<2.4;02.4Mod(t,100)TsQ_{j}\left(T_{s}\right)=\left\{\begin{array}{ll} 0 & \operatorname{Mod}(t, 100)-T_{s}<0 \\ 100\left(\operatorname{Mod}(t, 100)-T_{s}\right) , & 0 \leq \operatorname{Mod}(t, 100)-T_{s}<0.2 ; \\ 20 h & 0.2 \leq \operatorname{Mod}(t, 100)-T_{s}<2.2 ; \\ \left(240-100\left(\operatorname{Mod}(t, 100)-T_{s}\right)\right) h & 2.2 \leq \operatorname{Mod}(t, 100)-T_{s}<2.4 ; \\ 0 & 2.4 \leq \operatorname{Mod}(t, 100)-T_{s} \end{array}\right.

    这里使用了一个取余函数,非常精妙,确保每个周期的喷气状态与第一个喷气周期状态相等。

    代码实现
    # 出油函数 参数t为时刻,t0为出油时刻
    def out(t, t0):
        global m, P, ρ  # 全局变量,可做修改
        if t0 < (t % 100) < t0 + 0.2:  # 第一段0~0.2ms
            Q = (t - t0) % 100 * 100  # 在不同周期内时刻的流量
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        elif t0 + 0.2 <= (t % 100) < t0 + 2.2:  # 第二段0.2~2.2ms
            Q = 20
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        elif t0 + 2.2 <= (t % 100) <= t0 + 2.4:  # 第三段2.2~2.4ms
            Q = ((t - t0) * (-100) + 240) % 100
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        else:
            pass
    

    进喷气状态结合

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

    mink=1nPk100=mink=1nj=0k1E(Pj)Vρj(Rj(Tp)Qj(Ts))min\sum_{k=1}^{n}|P_{k}-100|=min\sum_{k=1}^{n}|\sum_{j=0}^{k-1}\frac{E(P_{j})}{V\rho _{j}}(R_{j}(T_{p})-Q_{j}(T_{s}))

    关于进喷气函数的编程实现

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

    关于遍历及优化的编程实现

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

    关于搜索结果图的做法

    稳定模型完整代码

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

    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    plt.scatter(x_pandas_list, y_pandas_list, s=1, label='真实数据', color='#FFA500')  # 绘制压强密度散点图(真实数据)调了橙色
    plt.plot(x_pandas_list, ρ_P(x_pandas_list), label='拟合数据', color='#1E90FF')  # 绘制压强密度图像(拟合数据)调了蓝色
    plt.title('燃油密度(mm3/ms)-压力(MPa)')
    plt.xlabel('燃油密度(mm3/ms)')
    plt.ylabel('压力(MPa)')
    plt.legend()  # 显示标注label
    plt.show()  # 画图显示
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.scatter(y_pandas_list, x_pandas_list, s=1, label='真实数据', color='#FFA500')  # 绘制压强密度散点图(真实数据)调了橙色
    plt.plot(ρ_P(x_pandas_list), x_pandas_list, label='拟合数据', color='#1E90FF')  # 绘制压强密度图像(拟合数据)调了蓝色
    plt.title('压力(MPa)-燃油密度(mm3/ms)')
    plt.xlabel('压力(MPa)')
    plt.ylabel('燃油密度(mm3/ms)')
    plt.legend()  # 显示标注label
    plt.show()  # 画图显示
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    ρ0 = 0.85  # 初始100MPa压力下的密度
    ρ160 = P_ρ(160)  # 160MPa压力下的密度 0.8711048440368855mg/mm3
    m0 = ρ0 * V_primary  # 初始高压油管油气33379.42194439156mg
    
    t = 0  # 记录时刻
    t_ = 0.01  # 每隔0.01ms时间段刷新一次状态
    m = m0  # 记录每一时刻高压油管内的质量 初始为m0
    P = 100  # 记录每一时刻高压油管内的压强
    ρ = 0.85  # 录每一时刻高压油管内的密度
    P_list = []  # 压强时刻表
    
    
    # 进油函数 参数t为时刻,t1为进油时间段
    def enter(t, t1):
        global m, P, ρ  # 全局变量,可做修改
        if 0 < t % (t1 + 10) < t1:
            Q = C * A * math.sqrt(2 * (160 - P) / ρ160)  # 单位时间进油量
            det_m = Q * ρ160 * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    
    # 出油函数 参数t为时刻,t0为出油时刻
    def out(t, t0):
        global m, P, ρ  # 全局变量,可做修改
        if t0 < (t % 100) < t0 + 0.2:  # 第一段0~0.2ms
            Q = (t - t0) % 100 * 100  # 在不同周期内时刻的流量
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        elif t0 + 0.2 <= (t % 100) < t0 + 2.2:  # 第二段0.2~2.2ms
            Q = 20
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        elif t0 + 2.2 <= (t % 100) <= t0 + 2.4:  # 第三段2.2~2.4ms
            Q = ((t - t0) * (-100) + 240) % 100
            det_m = Q * ρ * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ = m / V_primary  # 更新密度
            P = ρ_P(ρ)  # 更新压强
        else:
            pass
    
    
    sum_min = 10000000000000  # 刷新压力波动最小值
    P_list_min = []  # 波动最小的压力差表
    i_best = 0  # 最优开阀门时长
    j_best = 0  # 最优出气时刻
    
    for i in np.arange(0.287, 0.29, 0.001):  # 遍历开阀门时长 start->end
        for j in np.arange(50, 56.5, 1):  # 遍历出气时刻 start->end
            sum = 0  # 记录压力波动偏差
            t = 0  # 每次循环刷新时刻从0开始
            P_list = []  # 每次循环刷新压力时刻表 为空表
            while t <= 5000:  # 时间可修改,单位为ms 循环2000/t_次
                t = t + t_
                enter(t, i)
                out(t, j)
                P_list.append(P)
                sum = sum + abs(P - 100)  # 每隔t_=0.01ms时刻记录一下总波动 best=sum|P-100|
            print('开始遍历:', 'i =', i, 'j =', j)
            print('压差距离最优质值和(越小越好):', sum)
            if sum_min > sum:
                sum_min = sum
                i_best = i
                j_best = j
                P_list_min = P_list
    
    print('最优开阀门时长为:', i_best, '最优出气时刻为:', j_best)
    
    len_P_list = len(P_list)  # 长度为设值遍历长度个
    x_numpy_list = np.arange(0, len_P_list * t_, t_)  # 每隔0.01ms列一个x坐标与P_list对应
    
    plt.plot(x_numpy_list, P_list_min)
    plt.title('压强-时间')
    plt.xlabel('t(ms)')
    plt.ylabel('P(KPa)')
    plt.show()
    

    增压模型

    想要压力维持在150MPa150MPa上下,开启进气阀的时长需要设置在0.75ms0.75ms,这里0.75ms0.75ms是进气阀的时长遍历的结果。稍微修改一下上述代码的压强与密度初始值。

    已知0.75ms0.75ms作为150MPa150MPa下的开阀时长,现在要求的是,单向阀每次开启的时长多少时,经过2s,5s,10s2s,5s,10s分别到达150Mpa。

    ①2s(100MPa->150Mpa)

    2s2s前单向阀每次开启的时长0.89ms0.89ms,达到稳定后单向阀每次开启的时长0.75ms0.75ms

    ②5s(100MPa->150Mpa)

    5s5s前单向阀每次开启的时长0.75ms0.75ms,达到稳定后单向阀每次开启的时长0.75ms0.75ms

    ③10s(100MPa->150Mpa)

    8s8s前单向阀每次开启的时长0.50ms0.50ms,在8s10s8s-10s单向阀每次开启的时长0.75ms0.75ms,达到稳定后单向阀每次开启的时长0.75ms0.75ms

    关于这个每次开启的时长0.50ms0.50ms这个参数,可以修改,合理即可,没有唯一答案。

    问题二

    凸轮

    在上一问题中,进气状态所给压力恒定为160MPa160MPa。但在这里,高压油管的进气由凸轮转动提供,首先要明白凸轮是如何转动的,这是关键的一步。

    凸轮工作原理

    由这张图我们可以看出,凸轮的旋转是在二维平面中进行的。凸轮起始的位置为极径最大的位置,此时活塞位于柱塞腔的顶部,凸轮每转一周,向高压活塞压入一次气体,使高压油管的压力增大。这一过程的描述为:活塞由顶部运动到底部,此时不向高压油管内注入气体;当活塞由底部向顶部运动时,压力由0.5MPa0.5MPa增大,当压力增大到100MPa100MPa时,会向高压油管压入气体,此时柱塞腔燃气的质量减小,判断柱塞腔此时压力的变化。
    在这里插入图片描述

    关于凸轮形状

    凸轮形状在附件1-凸轮边缘曲线已经给出
    极角(弧度制)与极径的关系是

    R(θ)=2.39658307×103θ64.51694471×102θ5+2.59343668×101θ42.87191168×103θ39.49901781×101θ29.45855190×102θ+7.24715714R(\theta)=2.39658307 \times10^{-3}\theta^{6}-4.51694471 \times10^{-2}\theta^{5}+ 2.59343668 \times10^{-1}\theta^{4}-2.87191168 \times10^{-3}\theta^{3}-9.49901781 \times10^{-1}\theta^{2}-9.45855190\times10^{-2}\theta+7.24715714

    凸轮形状拟合代码
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    plt.scatter(x_pandas_list1, y_pandas_list1, s=0.1, label='真实数据', color='#FFA500')  # 绘制压强密度散点图(真实数据)调了橙色
    plt.plot(x_pandas_list1, y_(x_pandas_list1), label='拟合数据', color='#1E90FF')  # 绘制压强密度图像(拟合数据)调了蓝色
    plt.xlabel('x(mm)')
    plt.ylabel('y(mm)')
    plt.legend()
    plt.show()
    
    theta = np.arange(0, 2 * np.pi, 0.01)
    x = y_(theta) * np.cos(theta)
    y = y_(theta) * np.sin(theta)
    # 绘图
    plt.figure(figsize=(6, 6))
    plt.xlabel('x(mm)')
    plt.ylabel('y(mm)')
    plt.plot(x, y)
    plt.show()  # 显示凸轮形状
    

    柱塞腔活塞位置

    凸轮的顶部(每一极角下纵坐标的最大值)带动活塞的上升与下降,由此导致了柱塞腔燃气的密度、质量、压力变化。
    这个柱塞腔活塞位置其实是涉及到凸轮的角速度的,这个可以先估算一个大致区间然后带入一个数画图看效果。这个区间估算方法是 每个周期进气量大致等于喷气量,通过搜索算法确定步长,当压力差之和最小时就可以确定凸轮角速度。

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.027  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t_ = 0.1  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    t = 0  # 记录时刻
    H_list = []  # 记录每一时刻活塞距离腔底部的位置
    ρ_list = []  # 记录每一时刻活塞距离腔底部的密度
    P_list = []  # 记录每一时刻活塞距离腔底部的压力
    
    hMax = 0  # 记录每一时刻高度的最大值
    ρ_camj = P_ρ(0.5)  # 初始时刻压强
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    
    while t < T:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计
        ρ_list.append(cam_h(t))
        t = t + t_
    
    x_numpy_list = np.arange(0, len(ρ_list) * t_, t_)  # 选择
    plt.plot(x_numpy_list, ρ_list, label='h-t')
    plt.xlabel('时间/ms')
    plt.ylabel('柱塞腔燃油高度/KPa')
    plt.legend()
    plt.show()
    

    柱塞腔内燃气压力

    柱塞腔内压力随着活塞的上升而增大,在后半个周期内,压强首先增大到100MPa100MPa,然后在100MPa100MPa来回波动,大于100MPa100MPa的部分会在下一次迭代过程中将燃气压入高压油管,致使柱塞腔内的体积减小。

    ρcam j=mcamj 2.52πu(ωtj)\rho_{\text {cam } j}=\frac{m_{\text {camj }}}{2.5^{2} \pi u\left(\omega t_{j}\right)}

    mcamjm_{camj}为柱塞腔内燃气的质量,u(ωtj)u(\omega t_{j})为活塞的上升高度,参数为ω\omega

    P(ρ)=10784.65349198ρ215693.86029255ρ+5648.09019803P(\rho)=10784.65349198\rho^{2}-15693.86029255\rho+5648.09019803

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.027  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t_ = 0.01  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    t = 0  # 记录时刻
    H_list = []  # 记录每一时刻活塞距离腔底部的位置
    ρ_list = []  # 记录每一时刻活塞距离腔底部的密度
    P_list = []  # 记录每一时刻活塞距离腔底部的压力
    
    hMax = 0  # 记录每一时刻高度的最大值
    P_camj = 0.5  # 初始时刻压力
    ρ_camj = P_ρ(0.5)  # 初始时刻密度
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    def P_cam(t):  # 柱塞腔返回t时刻下的压力
        return ρ_P(ρ_cam(t))
    
    
    def ρ_cam(t):  # 柱塞腔返回t时刻下燃油的密度
        global m_camj, ρ_camj, P_camj
        if t % T <= T / 2:
            m_camj = m0_camj
            P_camj = 0.5
        else:
            if P_camj <= 100:
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
            else:
                det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_
                m_camj = m_camj - det_m
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
    
    
    while t <= 10 * T:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计
        t = t + t_
        ρ_cam(t)
        ρ_list.append(P_camj)
    
    x_numpy_list = np.arange(0, len(ρ_list) * t_, t_)  # 选择
    plt.plot(x_numpy_list, ρ_list, label='h-t')
    plt.xlabel('时间/ms')
    plt.ylabel('柱塞腔燃油压力/KPa')
    plt.legend()
    plt.show()
    

    流入气体质量

    ρcam j=mcamj 2.52πu(ωtj)\rho_{\text {cam } j}=\frac{m_{\text {camj }}}{2.5^{2} \pi u\left(\omega t_{j}\right)}

    Rj(ω)={0.85×π×0.722(Pcamj Ptubj )×ρcanj h,Pcamj >Ptubj 0,Pcanj Ptubj ,R_{j}(\omega)=\left\{\begin{array}{ll} 0.85 \times \pi \times 0.7^{2} \sqrt{2\left(P_{\text {camj }}-P_{\text {tubj }}\right) \times \rho_{\text {canj }}} h, & P_{\text {camj }}>P_{\text {tubj }} \\ 0, & P_{\text {canj }} \leqslant P_{\text {tubj }}, \end{array}\right.

    针阀运动

    针阀运动拟合曲线

    这两组参数是用polyfitpolyfit拟合出来的六次项系数,分别为分段函数的第一段与第三段系数,在实际程序设计中还要考虑到取余函数,周期为200ms200ms

    [ 1.59403345e+02 -9.87538422e+02 6.79811142e+02 -1.22515750e+02 1.16775838e+01 -4.11282406e-01 2.67578428e-03]
    [ 136.60587618 -1049.78616172 1227.80919016 10922.44930413-43933.87168977 63132.03256009 -32700.74327951]

    分段函数拟合代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    print(y_pandas_list)
    plt.xlabel('时间(ms)')
    plt.ylabel('距离(mm)')
    plt.scatter(x_pandas_list, y_pandas_list, s=0.1, label='真实数据', color='#FFA500')
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    h_s = 0
    Aj = 0
    
    
    def H(t):
        global h_s
        if 0 < t % 200 <= 0.45:
            h_s = cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                  * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            h_s = 2
        elif 2 <= t % 200 < 2.45:
            h_s = cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                  * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            h_s = 0
    
    
    h_s_list = []  # 记录针阀高度列表
    t = 0
    t_ = 0.01
    while t <= 2.45:
        H(t)
        h_s_list.append(h_s)
        t = t + t_
    
    x_np_list = np.arange(0, len(h_s_list) * 0.01, 0.01)  # 时间轴表示
    plt.plot(x_np_list, h_s_list, label='拟合数据', color='#1E90FF')
    plt.legend()
    plt.show()
    

    流出气体质量

    起初针阀向上运动时,针阀与底部开口面积较小,燃气可以全部流出。当针阀与底部开口面积小于底部孔面积时,这时燃气的流出面积为底部孔面积。

    Aj={π(1.25+H(tj)tan(π/20))21.252π,tj<0.330.72π,0.33tj<60.33π(1.25+H(6tj)tan(π/20))21.252π,60.33tj<60,6tj<100A_{j}=\left\{\begin{array}{ll} \pi\left(1.25+H\left(t_{j}\right) \tan (\pi / 20)\right)^{2}-1.25^{2} \pi, & t_{j}<0.33 \\ 0.7^{2} \pi, & 0.33 \leqslant t_{j}<\sqrt{6}-0.33 \\ \pi\left(1.25+H\left(\sqrt{6}-t_{j}\right) \tan (\pi / 20)\right)^{2}-1.25^{2} \pi, & \sqrt{6}-0.33 \leqslant t_{j}<\sqrt{6} \\ 0, & \sqrt{6} \leqslant t_{j}<100 \end{array}\right.

    AjA_{j}为流出的气体面积

    Qj(ω)=0.85×Aj2Ptubj ×ρtubj hQ_{j}(\omega)=0.85 \times A_{j} \sqrt{2 P_{\text {tubj }} \times \rho_{\text {tubj }}} h

    Qj(ω)Q_{j}(\omega)为单位时间内刘出气体质量

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    print(y_pandas_list)
    
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    A_s = 0
    
    
    def H(t):
        if 0 < t % 200 <= 0.45:
            return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                  * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            return 2
        elif 2 <= t % 200 < 2.45:
            return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                  * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            return 0
    
    
    A_list = []  # 记录针阀高度列表
    t = 0
    t_ = 0.01
    
    
    
    def Aj(t): #阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 <= t % 200 < math.sqrt(6) - 0.33:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6):
            A_s = math.pi * (1.25 + (H(math.sqrt(6)-t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    while t <= 2.45:
        H(t)
        Aj(t)
        A_list.append(A_s)
        t = t + t_
    
    x_np_list = np.arange(0, len(A_list) * 0.01, 0.01)  # 时间轴表示
    plt.xlabel('时间(ms)')
    plt.ylabel('面积(mm2)')
    plt.plot(x_np_list, A_list, label='实际流出面积', color='#1E90FF')
    plt.legend()
    plt.show()
    

    复杂模型

    高压油管内压力变化

    将问题二的进气状态与喷气状态结合,ω=0.027ω = 0.027效果较好,模型能够稳定在100MPa100MPa附近。

    测试代码(参数调整)
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.027  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t_ = 0.1  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    t = 0  # 记录时刻
    H_list = []  # 记录每一时刻活塞距离腔底部的位置
    ρ_list = []  # 记录每一时刻活塞距离腔底部的密度
    P_list = []  # 记录每一时刻活塞距离腔底部的压力
    
    hMax = 0  # 记录每一时刻高度的最大值
    P_camj = 0.5  # 初始时刻压力
    ρ_camj = P_ρ(0.5)  # 初始时刻密度
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    ρ0 = 0.85  # 初始100MPa压力下的密度
    ρ160 = P_ρ(160)  # 160MPa压力下的密度 0.8711048440368855mg/mm3
    m0 = ρ0 * V_primary  # 初始高压油管油气33379.42194439156mg
    
    t = 0  # 记录时刻
    t_ = 0.01  # 每隔0.01ms时间段刷新一次状态
    m = m0  # 记录每一时刻高压油管内的质量 初始为m0
    P = 100  # 记录每一时刻高压油管内的压强
    ρ = 0.85  # 录每一时刻高压油管内的密度
    P_list = []  # 压强时刻表
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    def P_cam(t):  # 柱塞腔返回t时刻下的压力
        return ρ_P(ρ_cam(t))
    
    
    def ρ_cam(t):  # 柱塞腔返回t时刻下燃油的密度
        global m_camj, ρ_camj, P_camj
        if t % T <= T / 2:
            m_camj = m0_camj
            P_camj = 0.5
        else:
            if P_camj <= 100:
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
            else:
                det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_
                m_camj = m_camj - det_m
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
    
    
    P_list = []  # 压强时刻表
    
    
    # 进油函数 参数t为时刻
    def enter(t):
        global m, P, ρ  # 全局变量,可做修改
        print(t)
        if P_camj > 100.1769:
            Q = C * A * math.sqrt(2 * (P_camj - 100.1769) / P_ρ(P_camj))  # 单位时间进油量
            det_m = Q * P_ρ(P_camj) * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    
    # 针阀拟合
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    A_s = 0
    
    
    def H(t):
        if 0 < t % 200 <= 0.45:
            return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                   * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            return 2
        elif 2 <= t % 200 < 2.45:
            return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                   * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            return 0
    
    
    def Aj(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 <= t % 200 < math.sqrt(6) - 0.33:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6):
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    while t <= 40 * T:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计 #ρ_cam(t) enter(t)      Aj(t) out(t)
        t = t + t_
        ρ_cam(t)
        enter(t)
        Aj(t)
        out(t)
        P_list.append(P + 2)
    
    len_P_list = len(P_list)  # 长度为设值遍历长度个
    x_numpy_list = np.arange(0, len_P_list * t_, t_)  # 每隔0.01ms列一个x坐标与P_list对应
    
    plt.plot(x_numpy_list, P_list)
    plt.title('压强-时间')
    plt.xlabel('t(ms)')
    plt.ylabel('P(KPa)')
    plt.show()
    

    上述代码未做优化遍历工作,实际应用还要遍历求最优值,程序跑起来较慢,可以先拟定一个初始值检测情况,再用比较良好的初始值在附近进行遍历求解。

    问题三

    两个喷油嘴模型

    增加一个喷油嘴,相比于问题二喷出气体的质量会增加,所以需要适当提高角速度,才能保持高压油管内压力稳定。调整参数
    第一个喷油嘴开始后的50ms50ms,第二个喷油嘴开启凸轮角速度调整为ω=0.0502ω = 0.0502

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.0502  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t = 0  # 记录时刻
    t_ = 0.1  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    
    hMax = 0  # 记录每一时刻高度的最大值
    P_camj = 0.5  # 初始时刻压力
    ρ_camj = P_ρ(0.5)  # 初始时刻密度
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    ρ0 = 0.85  # 初始100MPa压力下的密度
    ρ160 = P_ρ(160)  # 160MPa压力下的密度 0.8711048440368855mg/mm3
    m0 = ρ0 * V_primary  # 初始高压油管油气33379.42194439156mg
    
    m = m0  # 记录每一时刻高压油管内的质量 初始为m0
    P = 100  # 记录每一时刻高压油管内的压强
    ρ = 0.85  # 录每一时刻高压油管内的密度
    P_list = []  # 压强时刻表
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    def P_cam(t):  # 柱塞腔返回t时刻下的压力
        return ρ_P(ρ_cam(t))
    
    
    def ρ_cam(t):  # 柱塞腔返回t时刻下燃油的密度
        global m_camj, ρ_camj, P_camj
        if t % T <= T / 2:
            m_camj = m0_camj
            P_camj = 0.5
        else:
            if P_camj <= 100:
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
            else:
                det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_
                m_camj = m_camj - det_m
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
    
    
    P_list = []  # 压强时刻表
    
    
    # 进油函数 参数t为时刻
    def enter(t):
        global m, P, ρ  # 全局变量,可做修改
        print(t)
        if P_camj > 100.1769:
            Q = C * A * math.sqrt(2 * (P_camj - 100.1769) / P_ρ(P_camj))  # 单位时间进油量
            det_m = Q * P_ρ(P_camj) * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    
    # 针阀拟合
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    A_s = 0
    
    
    def H(t):
        if 0 < t % 200 <= 0.45:
            return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                   * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            return 2
        elif 2 <= t % 200 < 2.45:
            return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                   * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            return 0
    
    
    def Aj_1(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 <= t % 200 < math.sqrt(6) - 0.33:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6):
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_1(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    def Aj_2(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33 + 50:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 + 50 <= t % 200 < math.sqrt(6) - 0.33 + 50:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 + 50 <= t % 200 < math.sqrt(6) + 50:
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_2(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    while t <= 10000:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计 #ρ_cam(t) enter(t)      Aj(t) out(t)
        t = t + t_
        ρ_cam(t)
        enter(t)
        Aj_1(t)
        out_1(t)
        Aj_2(t)
        out_2(t)
        P_list.append(P + 4.5)
    
    len_P_list = len(P_list)  # 长度为设值遍历长度个
    x_numpy_list = np.arange(0, len_P_list * t_, t_)  # 每隔0.01ms列一个x坐标与P_list对应
    
    plt.plot(x_numpy_list, P_list)
    plt.title('压强-时间')
    plt.xlabel('t(ms)')
    plt.ylabel('P(KPa)')
    plt.show()
    

    两个喷油嘴与一个减压阀模型

    压力大于102MPa时开启减压阀

    这里我觉得是为了防止由于角速度过大而导致高压油管内部的气体压力一直处于上升状态,所以增加了减压阀。适当的把角速度增大并注意压力在某一个时刻开启减压阀。

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.07  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t = 0  # 记录时刻
    t_ = 0.1  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    
    hMax = 0  # 记录每一时刻高度的最大值
    P_camj = 0.5  # 初始时刻压力
    ρ_camj = P_ρ(0.5)  # 初始时刻密度
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    ρ0 = 0.85  # 初始100MPa压力下的密度
    ρ160 = P_ρ(160)  # 160MPa压力下的密度 0.8711048440368855mg/mm3
    m0 = ρ0 * V_primary  # 初始高压油管油气33379.42194439156mg
    
    m = m0  # 记录每一时刻高压油管内的质量 初始为m0
    P = 100  # 记录每一时刻高压油管内的压强
    ρ = 0.85  # 录每一时刻高压油管内的密度
    P_list = []  # 压强时刻表
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    def P_cam(t):  # 柱塞腔返回t时刻下的压力
        return ρ_P(ρ_cam(t))
    
    
    def ρ_cam(t):  # 柱塞腔返回t时刻下燃油的密度
        global m_camj, ρ_camj, P_camj
        if t % T <= T / 2:
            m_camj = m0_camj
            P_camj = 0.5
        else:
            if P_camj <= 100:
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
            else:
                det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_
                m_camj = m_camj - det_m
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
    
    
    P_list = []  # 压强时刻表
    
    
    # 进油函数 参数t为时刻
    def enter(t):
        global m, P, ρ  # 全局变量,可做修改
        print(t)
        if P_camj > 100.1769:
            Q = C * A * math.sqrt(2 * (P_camj - 100.1769) / P_ρ(P_camj))  # 单位时间进油量
            det_m = Q * P_ρ(P_camj) * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    
    # 针阀拟合
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    A_s = 0
    
    
    def H(t):
        if 0 < t % 200 <= 0.45:
            return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                   * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            return 2
        elif 2 <= t % 200 < 2.45:
            return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                   * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            return 0
    
    
    def Aj_1(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 <= t % 200 < math.sqrt(6) - 0.33:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6):
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_1(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    def Aj_2(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33 + 50:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 + 50 <= t % 200 < math.sqrt(6) - 0.33 + 50:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 + 50 <= t % 200 < math.sqrt(6) + 50:
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_2(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    def out_3(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        if P > 102:
            Q = C * 0.49 * math.pi * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
            det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ_camj = m / V_primary  # 更新密度
            P = ρ_P(ρ_camj)  # 更新压强
    
    
    while t <= 2000:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计 #ρ_cam(t) enter(t)      Aj(t) out(t)
        t = t + t_
        ρ_cam(t)
        enter(t)
        Aj_1(t)
        out_1(t)
        Aj_2(t)
        out_2(t)
        out_3(t)
        P_list.append(P)
    
    len_P_list = len(P_list)  # 长度为设值遍历长度个
    x_numpy_list = np.arange(0, len_P_list * t_, t_)  # 每隔0.01ms列一个x坐标与P_list对应
    
    plt.plot(x_numpy_list, P_list)
    plt.title('压强-时间')
    plt.xlabel('t(ms)')
    plt.ylabel('P(KPa)')
    plt.show()
    

    周期性开启减压阀

    刚开始处于上升状态而后逐渐趋于稳定。可以先把不开减压阀的状态图画出来再适当添加减压阀开启的周期,这样进气与出气状态平衡,高压油管内的压力也就自然平衡了。
    在这里插入图片描述

    代码实现
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import math
    
    data = '附件3-弹性模量与压力.xlsx'
    data = pd.read_excel('附件3-弹性模量与压力.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'燃油密度(mm3/ms)']  # 附件3 燃油密度列数据
    y_pandas_list = data[u'压力(MPa)']  # 附件3 压力列数据
    coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list,
                                         2)  # [ 10784.65349198 -15693.86029255   5648.09019803]y=压强 x=密度拟合系数
    
    
    def ρ_P(x):  # 密度转压力函数
        return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2]
    
    
    coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list,
                                         2)  # [-6.59843062e-07  5.23409641e-04  8.04251284e-01]y=密度 x=压强拟合系数
    
    
    def P_ρ(x):  # 压力转密度函数
        return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2]
    
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    data = '附件1-凸轮边缘曲线.xlsx'
    data = pd.read_excel(data)
    
    x_pandas_list1 = data[u'极角(rad)']
    y_pandas_list1 = data[u'极径(mm)']
    
    cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6)  # 凸轮 x-极角 y-极径拟合系数
    
    
    def y_(x):
        return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6]
    
    
    ω = 0.100  # 弧度制下的角速度rad/ms
    T = 2 * math.pi / ω  # 凸轮运动周期时间 单位ms
    t = 0  # 记录时刻
    t_ = 0.1  # 每隔0.01ms时间段刷新一次状态
    α0 = 0  # 初始时刻凸轮旋转角度0
    α = α0  # 记录凸轮旋转角度时刻
    
    hMax = 0  # 记录每一时刻高度的最大值
    P_camj = 0.5  # 初始时刻压力
    ρ_camj = P_ρ(0.5)  # 初始时刻密度
    d_camj = 5  # 柱塞腔直径5mm
    r_camj = d_camj / 2  # 柱塞腔半径2.5mm
    s_camj = math.pi * r_camj ** 2
    v_part = 20  # 向上残余容积为20mm3
    h_part = v_part / s_camj  # 向上残余高度
    m0_camj = 92.4217168684373  # 初始时刻柱塞腔油气质量
    m_camj = m0_camj  # 每一时刻柱塞腔油气质量
    
    C = 0.85  # 流量系数
    
    d_a = 1.4  # 小孔处的直径1.4mm
    r_a = d_a / 2  # 小孔处的半径1.4mm
    A = S_a = math.pi * r_a ** 2  # 小孔处的面积1.5393804002589984mm2
    
    l_primary = 500  # 内腔长度500mm
    d_primary = 10  # 内腔直径10mm
    r_primary = d_primary / 2  # 内腔半径5mm
    S_primary = math.pi * r_primary ** 2  # 内腔横截面积78.53981633974483mm2
    V_primary = S_primary * l_primary  # 内腔体积39269.90816987242mm3
    
    ρ0 = 0.85  # 初始100MPa压力下的密度
    ρ160 = P_ρ(160)  # 160MPa压力下的密度 0.8711048440368855mg/mm3
    m0 = ρ0 * V_primary  # 初始高压油管油气33379.42194439156mg
    
    m = m0  # 记录每一时刻高压油管内的质量 初始为m0
    P = 100  # 记录每一时刻高压油管内的压强
    ρ = 0.85  # 录每一时刻高压油管内的密度
    P_list = []  # 压强时刻表
    
    
    def cam_h(t):  # x为时刻 返回值为活塞距离柱塞腔底部的高度
        global α, hMax  # 全局变量可做修改
        hmax = 0  # 记录这一时刻下活塞向上运动的最大值
        for i in np.arange(0, 2 * math.pi, 0.1):  # 弧度制遍历
            h_point = y_(i) * math.cos(ω * t + i)  # 每一时间点下不同弧度坐标位置 h_point y坐标
            if h_point > hmax:  # 把每一个极角遍历出的最大值传递给hmax
                hmax = h_point
                hMax = 7.247588972988529 - hmax + h_part
        return hMax
    
    
    def P_cam(t):  # 柱塞腔返回t时刻下的压力
        return ρ_P(ρ_cam(t))
    
    
    def ρ_cam(t):  # 柱塞腔返回t时刻下燃油的密度
        global m_camj, ρ_camj, P_camj
        if t % T <= T / 2:
            m_camj = m0_camj
            P_camj = 0.5
        else:
            if P_camj <= 100:
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
            else:
                det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_
                m_camj = m_camj - det_m
                ρ_camj = m_camj / (s_camj * cam_h(t))  # 更新密度
                P_camj = ρ_P(ρ_camj)  # 更新压强
    
    
    P_list = []  # 压强时刻表
    
    
    # 进油函数 参数t为时刻
    def enter(t):
        global m, P, ρ  # 全局变量,可做修改
        print(t)
        if P_camj > 100.1769:
            Q = C * A * math.sqrt(2 * (P_camj - 100.1769) / P_ρ(P_camj))  # 单位时间进油量
            det_m = Q * P_ρ(P_camj) * t_  # 0.01为步长 质量改变量
            m = m + det_m  # 更新质量
            rou = m / V_primary  # 更新密度
            P = ρ_P(rou)  # 更新压强
        else:
            pass
    
    
    # 针阀拟合
    
    data = '附件2-针阀运动曲线.xlsx'
    data = pd.read_excel('附件2-针阀运动曲线.xlsx')  # 导入附件3数据
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    x_pandas_list = data[u'时间(ms)']  # 附件2 时间数据
    y_pandas_list = data[u'距离(mm)']  # 附件2 距离数据
    
    y1_pandas_list = y_pandas_list[0:46]
    x1_pandas_list = x_pandas_list[0:46]
    cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6)  # 第一段分段函数系数[15.39735971 -1.93086567  0.026043  ]
    
    y2_pandas_list = y_pandas_list[200:246]
    x2_pandas_list = x_pandas_list[200:246]
    cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6)  # 第三段分段函数系数[ 15.44372154 -73.71211028  87.92142494]
    A_s = 0
    
    
    def H(t):
        if 0 < t % 200 <= 0.45:
            return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \
                   * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6]
        elif 0.45 < t % 200 < 2:
            return 2
        elif 2 <= t % 200 < 2.45:
            return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \
                   * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6]
        else:
            return 0
    
    
    def Aj_1(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 <= t % 200 < math.sqrt(6) - 0.33:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6):
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_1(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    def Aj_2(t):  # 阻塞模型下的流出面积
        global A_s
        if t % 200 < 0.33 + 50:
            A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        elif 0.33 + 50 <= t % 200 < math.sqrt(6) - 0.33 + 50:
            A_s = 0.7 ** 2 * math.pi
        elif math.sqrt(6) - 0.33 + 50 <= t % 200 < math.sqrt(6) + 50:
            A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi
        else:
            A_s = 0
    
    
    def out_2(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        Q = C * A_s * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
        det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
        m = m - det_m  # 更新质量
        ρ_camj = m / V_primary  # 更新密度
        P = ρ_P(ρ_camj)  # 更新压强
    
    
    def out_3(t):
        global m, P, ρ_camj  # 全局变量,可做修改
        if 50 <= t % 200 < 10 + 50:
            Q = C * 0.49 * math.pi * math.sqrt(2 * P / ρ_camj)  # 单位时间进油量
            det_m = Q * ρ_camj * t_  # 0.01为步长 质量改变量
            m = m - det_m  # 更新质量
            ρ_camj = m / V_primary  # 更新密度
            P = ρ_P(ρ_camj)  # 更新压强
    
    
    while t <= 20000:  # 遍历凸轮旋转一个周期 遍历单位为t_可以自行设计 #ρ_cam(t) enter(t)      Aj(t) out(t)
        t = t + t_
        ρ_cam(t)
        enter(t)
        Aj_1(t)
        out_1(t)
        out_2(t)
        out_3(t)
        P_list.append(P)
    
    len_P_list = len(P_list)  # 长度为设值遍历长度个
    x_numpy_list = np.arange(0, len_P_list * t_, t_)  # 每隔0.01ms列一个x坐标与P_list对应
    
    plt.plot(x_numpy_list, P_list)
    plt.title('压强-时间')
    plt.xlabel('t(ms)')
    plt.ylabel('P(KPa)')
    plt.show()
    

    本题代码复用过程较多,设置变量时要记清。第二问可以嵌套第一问的代码,第三问更是在第二问的基础上进行代码的更改,第三问参数设置比较丰富,根据自己实际需求设置即可求得最终解。这里所有的PythonPython代码我都已经整理到了压缩包中分享,直接运行即可,链接:https://pan.baidu.com/s/1TZ-C3Ys66bIT6qLVlyTQDQ,需要请留言。

    本人学习用作记录,顺便分享Python代码,存在错误请批评指正

    展开全文
  • 以Hoek-Brown准则为极限平衡条件,提出了基于工程允许塑性区半径深部煤层井眼坍塌压力弹塑性计算方法,得到控制深部煤层井壁坍塌主要参数:煤岩地质强度指标GSI值、非均匀地应力系数及工程允许塑性区半径。...
  • 我们引入一个称为道德纽带的新概念,该概念将座席的行为网络指标结合起来,并可以量化座席行使和协调同伴压力的能力。 道德纽带对同伴压力的影响是非线性的,并且小组人数无关。 只有关于道德纽带的同构网络会...
  • 基于煤层的开采引起上覆岩层大范围移动和应力重分布,尤其是工作面周围煤岩体的承载应力变化这一理论,运用ANSYS对无煤...同时由小煤柱的宽度应力集中系数的关系发现:柱宽越大,相同条件下的集中系数越小,并且增幅越小。
  • 该研究的更大目的是阐明工作条件与压力和工作满意度之间的关系。 t检验表明,与我们的主要假设和先前的研究相反,两个性别之间的压力水平没有差异。 我们使用了Pearson r系数,该系数既未确认负面工作条件之间的...
  • 理论分析和实验测量了几何结构对压力传感特性影响规律, 结果表明, 熊猫光纤光栅压力灵敏度系数与猫眼半径距离比平方成线性关系。对比实验结果表明, 无法通过测量熊猫光纤光栅双峰间距方案实现对温度不敏感...
  • 依据Wroth在一维固结条件下提出的不排水强度超固结比(OCR)的关系式,利用快速直剪的剪切强度计算临界状态孔隙水压力系数Λ0,从而计算出OCR以及先期固结压力。对南京某黏土的重塑饱和黏土的试验分析结果表明,快剪...
  • 为了研究综放工作面矿山压力显现对瓦斯涌出量影响,以赵庄二号井1304综放工作面为工程背景,采用在线监测方法,对工作面超前支承压力、顶板周期来压和瓦斯涌出量进行实时监测。通过分析实时监测数据得到:工作面前方...
  • 改善大采高工作面单孔注水效果,结合赵固二矿11011工作面实际情况,采用数值模拟方法,得出不同注水压力下孔隙水压变化特征,根据孔隙水压与工作面前方煤体应力、渗透系数与注水压力之间参数关系,确定赵固二矿...
  • 目前,低渗透矩形油藏地层渗流压力求解的研究比较成熟,但是对于渗流压力的动态分析及渗流压力影响因素分析的研究尚存在一些不足。针对低渗透矩形油藏难动用的难题,基于两参数连续模型,运用稳态逐次替换法,分别建立了...
  • 通常我们所说和所了解的应力分析是泛指弹性分析,即应力应变之间的关系始终是线性的,符合胡克定律,通过有限元软件求出对应载荷下的应力并采用应力分类法并给予一定的安全系数进行应力强度评定。其实这种方法同常规...

    ANSYS分析设计人—专注压力容器分析设计的交流平台!学贵得师,更贵得友!共同学习,共同进步!

    5e5b387da73233bf199acfa1ba3f4f30.png

    通常我们所说和所了解的应力分析是泛指弹性分析即应力应变之间的关系始终是线性的符合胡克定律通过有限元软件求出对应载荷下的应力并采用应力分类法并给予一定的安全系数进行应力强度评定其实这种方法同常规设计一样也是基于理论和实践经验总结出来的一种方法只不过是与常规设计采用的理论不同而已进而在安全系数取值计算方法强度评定等等方面衍生出了一种新的规则

    自从ASME Ⅷ-2引入弹塑性分析方法以来越来越多的学者和工程师已经这种方法开始应用于国内市场和工程实际中目前国内分析设计标准JB4732修订版征求意见稿中也已引入了非线性分析的极限载荷法和弹塑性分析法在工程实际中大多数材料都是弹塑性状态下工作的而弹塑性分析正是采用材料的真实应力应变曲线可计算整个时间历程中的弹性应变和塑性应变变化情况与弹性分析法相比弹塑性分析更加精确和接近工程实际且在大多数情况下弹塑性应力分析法能节省材料成本但是其在前处理求解设置以及后处理等操作过程中相对复杂一些而且对分析设计人员和计算机的配置要求也较高做好弹塑性分析的前提一是对弹塑性概念和理论的深刻理解二是将这些理论很好的通过有限元软件来实现对有限元软件的理解和操作也必须深入和灵活将理论和软件合二为一融会贯通二者缺一不可

    b9a6cc270af0edd6e865163cc7846b6c.gif弹塑性分析的本构模型和塑性理论准则

    (1)本构模型弹塑性分析法采用考虑应变强化的真实应力应变曲线来建立材料的本构模型采用大变形理论刚度矩阵和平衡方程一直在更新变化因而属于非线性分析求解时间大大增加且存在求解收敛问题

    (2)屈服准则弹塑性分析基于一定的屈服准则来判定某种应力状态下的材料是处于弹性范围内还是已经进入塑性流动状态初始屈服条件则规定了材料开始进入塑性变形的应力状态目前关于塑性理论的屈服评判准则有多种但最常用的关于金属材料的有两种:Mises屈服准则和Tresca屈服准则这两种屈服条件的差别不是很大通常Tresca屈服条件更安全一些Mises屈服条件则应用起来更为方便因此在有限元分析中通常采用Mises屈服准则

    (3)流动准则流动准则是用来描述塑性应变张量增量的分量和应力分量以及应力增量分量之间的关系并在此基础上建立弹塑性本构关系表达式通俗的讲,就是材料在进入塑性状态后,材料的塑性变形在应力状态(应力分量和应力增量)中的流动规律。

    (4)硬化准则硬化准则规定材料进入塑性变形后的后继屈服函数(又称加载函数或加载曲面)的形式对于理想弹塑性材料由于没有硬化效应后继屈服函数和初始屈服函数是一致的而对于硬化材料有限元软件中都提供了多种硬化准则供用户根据材料特性和工程实际自行选择

    (5)加载卸载准则加载卸载准则主要用来判别从某一塑性状态出发材料是处于塑性加载状态还是弹性卸载状态在判定材料是否继续塑性变形究竟采用弹塑性本构关系还是弹性本构关系时加载卸载准则是必须的

    上述准则我们仅需了解即可实际上这部分在有限元软件中是傻瓜式的操作除过上述材料的真实应力-应变曲线需要我们自己查询手动输入外其余准则的实现都是通过选择软件中自带的材料模型来实现的这些材料模型就对应着不同的准则如各向同性硬化模型双线性随动强化模型等 

    10fe38093b538054e8ca711bf58b10e4.png

    注意:1)ANSYS经典中可选择输入总应变或塑性应变而在WB中只能输入塑性应变;2)材料的应力-应变曲线必须是单调递增的不允许出现颈缩阶段开始递减的数据点(真实的材料是有颈缩阶段的),否则会极大可能出现不收敛的问题;3)要让塑性数据最后一行中的塑性应变大于模型中可能出现的最大塑性应变值可在最后一行中增加一行将其中的塑性应变设置为一个比较大的值并相应的选择此应变下的真实应力值并使得曲线倾斜向上目的一是可保证整个分析过程中都使用了硬化模型二是可避免收敛困难问题

    b9a6cc270af0edd6e865163cc7846b6c.gif有限元建模计算过程中遵循以下基本原则

    (1)模型建立应尽量建出有可能出现较大塑性变形位置的模型细节尺寸但对有可能出现应力奇异的位置要简化和优化模型

    (2)网格密度在塑性变形较大的区域细化网格如果网格过于粗糙相邻单元之间的应力和应变变化会出现不连续的跳跃现象会造成难以收敛的问题但是网格也不能过细过细的网格也可能导致收敛困难所以关于网格密度和网格质量的问题需要通过计算过程和对计算结果的判定来进行不断调整

    (3)网格质量在塑性变形较大区域需划分高质量的网格避免出现过大的钝角或过小的尖角不让单元的形状过于狭小

    (4)单元选择尽量不要选择二阶完全积分单元容易出现体积自锁二阶减缩积分单元需要划分足够密的网格才不会产生体积自锁因而建议使用一阶减缩积分单元关于单元的介绍可看如下链接内容

    默认被忽视的问题却往往最可能成为致命问题—论单元选择的重要性!

    (4)应避免应力奇异应力奇异经常出现的区域单点加载或单点约束凹角模型之间采用单点连接单点耦合或接触条件等

    (5)采用大变形理论打开大变形开关

    (6)载荷步设置:在加载过程中设置足够多的子步数等比例逐渐施加载荷并保证在一个时间步内最大的塑性应变增量小于5%;载荷步的设置不仅影响到计算结果甚至会计算是否会收敛因而载荷步的设置是一个需要摸索和经验的活如载荷步设置的较少则计算可能发散若载荷步设置的过多则计算时间有可能会大大增加可通过如下收敛曲线初步判断计算的收敛性 

    2973cbf97e10290a71c1413ac6bd1871.png

    关于软件的设置载荷步的施加及与极限载荷分析的区别可看如下链接内容

    极限载荷分析-你想知道的也许在这里能找到些许答案

    (7)不能只关注计算是否收敛还应关注应力应变塑性应变等对加载的时间历程曲线是否光滑若出现不光滑则说明时间步长太大或单元网格太疏则计算结果是不可信的如可通过应变变化曲线的光滑性和塑性应变增量小于5%初步判定结果的正确性又可通过应力应变图基本与材料本构模型中的应力应变曲线相一致可进一步判断结果的正确性 

    0906d02777aeafb3234e6def0c8f9d9c.png

    f254dda19cba233fbe39e733a13b4cd2.png

    b9a6cc270af0edd6e865163cc7846b6c.gif弹塑性分析评定方法

    弹塑性分析的目的一是防止发生总体塑性垮塌二是防止局部产生过度应变。JB4732征求意见稿引进ASME多种载荷组合工况的计算法并分别进行总体塑性垮塌的评定和局部过度应变的评定

    对总体塑性垮塌的评定:可采用载荷系数法或塑性垮塌载荷法按如下步骤进行评定当采用载荷系数法时需对每种载荷组合工况乘以相应的载荷系数并进行弹塑性分析每种组合工况均计算收敛则评定合格和通过当采用塑性垮塌载荷法时同样需对每种载荷组合工况均进行弹塑性分析采用较小的载荷增量步加载若加载到第K步时计算发散则第k-1步施加的载荷即为垮塌载荷K-1步得到的垮塌载荷除以安全系数2.4得到许用载荷若设计载荷小于等于许用载荷则评定通过弹塑性分析中的两种评定方法流程示意图如下 

    696a8fc30d1e14e0fbfd8931ab131e66.png

    e80db3ed1b32f4bae00793cc4fa631ce.png

    如上图是采用载荷系数法计算并通过计算收敛性来进行评定的计算结果收敛且等效塑性应变约为0.019mm。总体塑性垮塌评定合格和通过

    对局部过度应变的评定:可通过有限元软件弹塑性分析通过如下公式确定总当量塑性应变确定三轴应变极限确定成形应变的方法来进行评定评定方法流程示意图如下 

    2b2f627f7e6385af9554b4f959942d7d.png

    75d97d10dbeaa8f04e83c871a2d32f69.png

    5ecb13608fe83ad062e0b1130dcfe492.pnga73719c1c244e78e107efe8cc125fb11.png

    如上图通过有限元软件AWB分别计算了三轴应变极限总的当量塑性应变与成形应变之和与三轴应变极限比值(εpeq+εef) /εL 的分布云图从图中可以看出该比值的最大值为0.1087,小于1,(εpeq+εef))故结构满足该组合载荷工况下防止局部过度应变的要求评定合格和通过AWB中通过User Defined Result可很容易的进行函数的定义来求解三轴应变极限总的当量塑性应变与成形应变之和与三轴应变极限比值但要注意定义函数时里面的自变量需采用AWB内置的且能识别的简称(如等效塑性应变在AWB中的简称是EPPLEQ_RST)。

    上述是对弹塑性分析在有限元软件AWB中实现的一个简单步骤的介绍实际操作过程中有很多需要注意的地方一个地方出错可能会导致满盘皆输弹塑性分析是一个建立在对理论的理解和经验的基础上且需要不断摸索的过程因是非线性分析就会存在最大的一个问题时间性和收敛性而计算能否收敛和能否提高计算效率则取决于很多因素包括模型网格求解设置等多方面因素均会影响最终计算的时间性和收敛性虽然弹塑性分析已引入国内但笔者以为要想在短时间内取代弹性分析的应力分类法还有需要较长时间目前只能作为应力分类法在局部小模型中的一种辅助验证方法得以应用因弹塑性分析对设计人员的理论和操作水平计算结果计算效率计算硬件计算成本要比弹性线性分析要求高得多(比如上述的简单开孔接管结构采用弹塑性分析在一台高配置的电脑上计算时间花了将近三个半小时而如果采用弹性分析的话在高配置的电脑上计算时间可能仅需一分钟孰轻孰重一目了然)且不论其它方面的一些因素仅仅对计算机配置要求很高这一点就足以成为其得以广泛应用的一大限制和阻碍 

    以上是笔者的一点个人拙见如有不当之处还请不吝批评指正后续笔者会对弹塑性分析的理论和计算过程中遇到的一些常见问题进行归纳总结与朋友们一起分享学习 

    154f963713e8eba527a70c29b173786b.gif

    请不吝点个在看!分享成就你我他!

    2a4d603e5078a16f7bfd2b27e337d049.png

    在这里,我们愿与您一起,亦师亦友,共同学习,共同进步; 期待有志者的加入!

    展开全文
  • 基于活度系数模型气液平衡计算

    千次阅读 2017-10-18 09:01:33
    Property Analysis计算优点是便于生成p-xy,T-xy和描述混合物不同组分吉布斯自由能之间函数关系的图表。也可用于研究共沸物形成可能性。缺点是只能用于双组分混合物计算,只用于固定温度或固定压力计算...

    Aspen中气液平衡的分析方法

    VLE:气液平衡

    在aspen中,VLE计算有三种方法。

    Property Analysis计算的优点是便于生成p-xy,T-xy和描述混合物不同组分与吉布斯自由能之间的函数关系的图表。也可用于研究共沸物形成的可能性。缺点是只能用于双组分混合物计算,只用于固定温度或固定压力下的计算,而不能用于仅知道最终压力和焓值条件下的计算。

    Simulation适用范围更广,但必须依赖流程图和分离单元模块,建立计算模型。



    气液平衡的热力学基础

    参考:《化工原理(下)》 清华大学出版社 P109 蒸馏部分

    根据相率,平衡体系自由度为:

    F=C-P+2

    其中F为自由度,即在不引起相变的条件下可以变动的独立变量的书目。C为独立组分数。P为相数。对于两组分混合系统,其自由度为2,即给定压强和温度,其平衡状态就确定了。因此可以用平面二维图像:T-X图,P-X图,X-Y图来刻画两组分混合系统平衡状态。

    对T-x图,x-y图的分析详见《化工原理》。

    混合物组分越多,自由度越多,所需参数越多,因此对于多组分体系不能用简单的等温图和等压图表示。



    判断混合物已达到平衡有三条原则:1. 两相温度相同 2. 两相压力相同 3.各组分两相中的化学势必须相等(或者各组分两相中的逸度必须相等)

    其中第三条比较重要,即为相平衡准则,可从热力学第二定律导出。

    化学势μ的定义为:偏摩尔自由焓,即在温度、压力以及除i组分外其余任意组分j的摩尔数不变时,自由焓G对i组分摩尔数ni的偏导数;其中G=U+pV-TS,U为内能,S为熵,V为体积。

    逸度的定义为:,其中μ为化学势,f为逸度。

    应注意,逸度本身由化学势定义,而对于理想气体,可化简为P/P'。





    展开全文
  • 研究结果表明:在保持围压不变情况下,煤样孔隙率随轴压增大呈线性趋势减小,并且不同方向孔隙率压力影响系数表现出明显各向异性;煤样声波速度随轴压增大呈线性趋势增大,且煤岩超声波速具有明显层...

空空如也

空空如也

1 2 3 4 5 ... 8
收藏数 159
精华内容 63
关键字:

压力系数与压力的关系