# ===================================================================== # 2019A 高压油管的压力控制 (CUMCM 2019 Problem A) -- standalone source # Self-contained physics: compressibility (注1), orifice flow (注2), # cam/plunger pump (问题2), dual-nozzle + relief valve (问题3). # Reads official appendix tables from assets/problems/: # app1_cam.csv (凸轮边缘曲线), app2_needle.csv (针阀升程), app3_emod.csv (弹性模量) # All numbers used by figures / papers / appendix code come from here. # ===================================================================== import os, math, csv as _csv ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) OUT = os.path.join(ROOT, "assets", "problems", "data") os.makedirs(OUT, exist_ok=True) # ---- geometry & operating constants (from problem statement) ---- A19_L_TUBE = 500.0 # mm A19_D_TUBE = 10.0 # mm A19_V_TUBE = math.pi * (A19_D_TUBE/2)**2 * A19_L_TUBE # mm^3 A19_A_IN = math.pi * (1.4/2)**2 # inlet / relief hole area (d=1.4 mm), mm^2 A19_P_SUPPLY = 160.0 # MPa A19_P0 = 100.0 # initial tube pressure, MPa A19_C = 0.85 # flow coefficient A19_T_CLOSE = 10.0 # ms valve closed after each open A19_T_INJ_PERIOD = 100.0 # ms (10 / s) A19_T_INJ = 2.4 # ms injection duration A19_INJ_PEAK = 20.0 # mm^3/ms A19_D_PLUNGER = 5.0 A19_A_PLUNGER = math.pi * (A19_D_PLUNGER/2)**2 A19_V_RES = 20.0 # mm^3 residual volume at top A19_P_LOW = 0.5 # MPa low-pressure supply / relief return A19_RHO_REF_P = 100.0 A19_RHO_REF = 0.850 # mg/mm^3 at 100 MPa def _load_app(fn): rows = [] with open(os.path.join(ROOT, "assets", "problems", fn)) as f: rd = _csv.reader(f); next(rd) for r in rd: if r[0] is None: continue rows.append((float(r[0]), float(r[1]))) return rows # ---- elastic modulus E(p) and density rho(p) from 附件3 ---- _A19_EMOD = _load_app("app3_emod.csv") A19_E_P = [r[0] for r in _A19_EMOD] A19_E_E = [r[1] for r in _A19_EMOD] def a19_E_of(p): if p <= A19_E_P[0]: return A19_E_E[0] if p >= A19_E_P[-1]: return A19_E_E[-1] lo, hi = 0, len(A19_E_P)-1 while hi - lo > 1: mid = (lo+hi)//2 if A19_E_P[mid] <= p: lo = mid else: hi = mid t = (p - A19_E_P[lo])/(A19_E_P[hi]-A19_E_P[lo]) return A19_E_E[lo] + t*(A19_E_E[hi]-A19_E_E[lo]) # integrate d rho/rho = dp/E(p) (注1) def _a19_dInt(a, b): n = max(2, int((b-a)/0.05)+1); s = 0.0; pa = a for k in range(1, n+1): pb = a + (b-a)*k/n; s += (pb-pa)/a19_E_of(0.5*(pa+pb)); pa = pb return s _A19_P, _A19_RHO = [A19_RHO_REF_P], [A19_RHO_REF] p = A19_RHO_REF_P + 0.05 while p <= 200.0 + 1e-9: _A19_RHO.append(_A19_RHO[-1]*math.exp(_a19_dInt(p-0.05, p))); _A19_P.append(p); p += 0.05 p = A19_RHO_REF_P - 0.05 while p >= 0.0 - 1e-9: _A19_RHO.append(_A19_RHO[0]*math.exp(-_a19_dInt(p, p+0.05))); _A19_P.append(p); p -= 0.05 _A19_PR = sorted(zip(_A19_P, _A19_RHO)); A19_P_ARR = [x[0] for x in _A19_PR]; A19_RHO_ARR = [x[1] for x in _A19_PR] def a19_rho_of(p): if p <= A19_P_ARR[0]: return A19_RHO_ARR[0] if p >= A19_P_ARR[-1]: return A19_RHO_ARR[-1] lo, hi = 0, len(A19_P_ARR)-1 while hi - lo > 1: mid = (lo+hi)//2 if A19_P_ARR[mid] <= p: lo = mid else: hi = mid t = (p - A19_P_ARR[lo])/(A19_P_ARR[hi]-A19_P_ARR[lo]) return A19_RHO_ARR[lo] + t*(A19_RHO_ARR[hi]-A19_RHO_ARR[lo]) def a19_p_of_rho(rho): if rho <= A19_RHO_ARR[0]: return A19_P_ARR[0] if rho >= A19_RHO_ARR[-1]: return A19_P_ARR[-1] lo, hi = 0, len(A19_RHO_ARR)-1 while hi - lo > 1: mid = (lo+hi)//2 if A19_RHO_ARR[mid] <= rho: lo = mid else: hi = mid t = (rho - A19_RHO_ARR[lo])/(A19_RHO_ARR[hi]-A19_RHO_ARR[lo]) return A19_P_ARR[lo] + t*(A19_RHO_ARR[hi]-A19_RHO_ARR[lo]) # ---- orifice flow (注2, literal competition formula) ---- # Q[mm^3/ms] = C*A[mm^2]*sqrt(2*dP[MPa]/rho[mg/mm^3]) def a19_q_orifice(A, p_high, p_low, rho_high_mgmm3): dP = p_high - p_low if dP <= 0: return 0.0 return A19_C * A * math.sqrt(2.0 * dP / rho_high_mgmm3) def a19_q_inj(t): if t < 0 or t > A19_T_INJ: return 0.0 if t <= 0.2: return A19_INJ_PEAK*(t/0.2) if t <= 2.2: return A19_INJ_PEAK return A19_INJ_PEAK*(1.0 - (t-2.2)/0.2) def a19_valve_open(t, topen): cyc = topen + A19_T_CLOSE return (t % cyc) < topen # ============================ 问题1 ============================ def a19_sim_p1(topen, T_total=4000.0, dt=0.02): n = int(T_total/dt); p = A19_P0; ts, ps = [], [] for k in range(n+1): t = k*dt qout = a19_q_inj(t % A19_T_INJ_PERIOD) if a19_valve_open(t, topen) and p < A19_P_SUPPLY: rho_s = a19_rho_of(A19_P_SUPPLY); rho = a19_rho_of(p) qin = a19_q_orifice(A19_A_IN, A19_P_SUPPLY, p, rho_s) dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(rho_s*qin - rho*qout) else: rho = a19_rho_of(p); dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(0 - rho*qout) p = p + dpdt*dt if k % max(1,int(20/dt)) == 0: ts.append(t); ps.append(p) return ts, ps def a19_mean_std(ps, tail=200): tp = ps[-tail:]; m = sum(tp)/len(tp) var = sum((x-m)**2 for x in tp)/len(tp) return m, math.sqrt(var) def a19_find_topen(target, lo=0.05, hi=3.0, T_total=4000.0, dt=0.02, iters=30): for _ in range(iters): mid = 0.5*(lo+hi) _, ps = a19_sim_p1(mid, T_total, dt) m, _ = a19_mean_std(ps) if m < target: lo = mid else: hi = mid return 0.5*(lo+hi) def a19_sim_p1_schedule(schedule, T_total=12000.0, dt=0.02): n = int(T_total/dt); p = A19_P0; ts, ps = [], []; seg = 0 for k in range(n+1): t = k*dt while seg < len(schedule)-1 and t > schedule[seg][0]: seg += 1 topen = schedule[seg][1] qout = a19_q_inj(t % A19_T_INJ_PERIOD) if a19_valve_open(t, topen) and p < A19_P_SUPPLY: rho_s = a19_rho_of(A19_P_SUPPLY); rho = a19_rho_of(p) qin = a19_q_orifice(A19_A_IN, A19_P_SUPPLY, p, rho_s) dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(rho_s*qin - rho*qout) else: rho = a19_rho_of(p); dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(0 - rho*qout) p = p + dpdt*dt if k % max(1,int(20/dt)) == 0: ts.append(t); ps.append(p) return ts, ps def a19_time_to_reach_150(topen, T_total=20000.0, dt=0.05): n = int(T_total/dt); p = A19_P0 for k in range(1, n+1): t = k*dt qout = a19_q_inj(t % A19_T_INJ_PERIOD) if a19_valve_open(t, topen) and p < A19_P_SUPPLY: rho_s = a19_rho_of(A19_P_SUPPLY); rho = a19_rho_of(p) qin = a19_q_orifice(A19_A_IN, A19_P_SUPPLY, p, rho_s) dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(rho_s*qin - rho*qout) else: rho = a19_rho_of(p); dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(0 - rho*qout) p = p + dpdt*dt if p >= 150.0: return t return None def a19_crossing(schedule, T_total=15000, dt=0.05): n = int(T_total/dt); p = A19_P0; seg = 0 for k in range(1, n+1): t = k*dt while seg < len(schedule)-1 and t > schedule[seg][0]: seg += 1 topen = schedule[seg][1] qout = a19_q_inj(t % A19_T_INJ_PERIOD) if a19_valve_open(t, topen) and p < A19_P_SUPPLY: rho_s = a19_rho_of(A19_P_SUPPLY); rho = a19_rho_of(p) qin = a19_q_orifice(A19_A_IN, A19_P_SUPPLY, p, rho_s) dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(rho_s*qin - rho*qout) else: rho = a19_rho_of(p); dpdt = (a19_E_of(p)/(rho*A19_V_TUBE))*(0 - rho*qout) p = p + dpdt*dt if p >= 150.0: return t return None def a19_find_charge(T_target, topen150, iters=30): lo, hi = topen150, 3.0 for _ in range(iters): mid = 0.5*(lo+hi) tc = a19_crossing([(T_target, mid), (1e9, topen150)]) if tc is None or tc > T_target: lo = mid else: hi = mid return 0.5*(lo+hi) # ============================ 问题2/3: cam & plunger ============================ _A19_CAM = _load_app("app1_cam.csv") A19_CAM_TH = [r[0] for r in _A19_CAM] A19_CAM_R = [r[1] for r in _A19_CAM] A19_CAM_RMIN = min(A19_CAM_R); A19_CAM_H = max(A19_CAM_R) - A19_CAM_RMIN A19_CAM_PERIOD = A19_CAM_TH[-1] def a19_L_of(phi): phi_eff = phi % A19_CAM_PERIOD idx = phi_eff/0.01; i0 = int(idx); frac = idx - i0 if i0 >= len(A19_CAM_R)-1: return A19_CAM_R[-1] - A19_CAM_RMIN r = A19_CAM_R[i0] + frac*(A19_CAM_R[i0+1]-A19_CAM_R[i0]) return r - A19_CAM_RMIN def a19_sim_p2(omega, T_total=3000.0, dt=0.02, relief=None, nozzles=1): # chamber: dead volume V_RES stays at tube pressure; swept volume # V_sweep = A_PLUNGER*(H - L(phi)) filled with low-pressure fuel on intake, # expelled through the check valve on compression (volume-balanced delivery). n = int(T_total/dt); p = A19_P0; phi = 0.0 L0 = a19_L_of(phi); Vsw = A19_A_PLUNGER*(A19_CAM_H - L0); m_sweep = 0.0 ts, ps, pcs = [], [], [] for k in range(n+1): t = k*dt phi_new = (phi + omega*dt/1000.0) % A19_CAM_PERIOD L1 = a19_L_of(phi_new); dL = L1 - L0; Vsw = A19_A_PLUNGER*(A19_CAM_H - L1) if dL < 0: # intake: swept volume fills with low-pressure fuel m_sweep = a19_rho_of(A19_P_LOW)*Vsw; pc = A19_P_LOW; q_in_ch = 0.0 else: # compression pc_closed = a19_p_of_rho(m_sweep/Vsw) if Vsw > 1e-9 else A19_P_SUPPLY q_in_ch = 0.0; pc = pc_closed if pc_closed >= p and p < A19_P_SUPPLY: dLdt = dL/dt q_in_ch = A19_A_PLUNGER * dLdt pc = p + 0.5*a19_rho_of(p)*(q_in_ch/(A19_C*A19_A_IN))**2 m_sweep = m_sweep - a19_rho_of(pc)*q_in_ch*dt q_relief = 0.0 if relief is not None and p > relief["p_high"]: q_relief = a19_q_orifice(A19_A_IN, p, relief.get("p_return", A19_P_LOW), a19_rho_of(p)) qout = nozzles * a19_q_inj(t % A19_T_INJ_PERIOD) rho_t = a19_rho_of(p) dm = (a19_rho_of(pc) if q_in_ch > 0 else rho_t)*q_in_ch - rho_t*qout - rho_t*q_relief dpdt = (a19_E_of(p)/(rho_t*A19_V_TUBE)) * dm p = p + dpdt*dt L0 = L1; phi = phi_new if k % max(1,int(20/dt)) == 0: ts.append(t); ps.append(p); pcs.append(pc) return ts, ps, pcs def a19_mean_p(ps, tail=200): tp = ps[-tail:]; return sum(tp)/len(tp) def a19_find_omega(target=100.0, lo=10.0, hi=90.0, T_total=3000.0, dt=0.02, iters=28, relief=None, nozzles=1): for _ in range(iters): mid = 0.5*(lo+hi) _, ps, _ = a19_sim_p2(mid, T_total, dt, relief=relief, nozzles=nozzles) m = a19_mean_p(ps) if m < target: lo = mid else: hi = mid return 0.5*(lo+hi) def gen_2019a(seed=2019, out_csv="cumcm2019a.csv"): # ---- Problem 1 ---- topen100 = a19_find_topen(100.0) topen150 = a19_find_topen(150.0) t1_100 = a19_sim_p1(topen100, 4000, 0.02) t1_150 = a19_sim_p1(topen150, 4000, 0.02) t1_wrong = a19_sim_p1(0.50, 4000, 0.02) # wrong open time -> drift T_nat = a19_time_to_reach_150(topen150) # topen sensitivity sweep topen_sweep = [] for tp in [0.10,0.15,0.20,0.25,0.29,0.35,0.5,0.7,0.816,1.0,1.5,2.0]: _, ps = a19_sim_p1(tp, 4000, 0.02); m, s = a19_mean_std(ps) topen_sweep.append((round(tp,3), round(m,3), round(s,4))) # ramp schedules tc2 = a19_find_charge(2000.0, topen150) sch2 = [(2000.0, tc2), (1e9, topen150)] sch5 = [(5000.0 - T_nat, topen100), (1e9, topen150)] sch10 = [(10000.0 - T_nat, topen100), (1e9, topen150)] t_r2 = a19_sim_p1_schedule(sch2, 6000, 0.05) t_r5 = a19_sim_p1_schedule(sch5, 12000, 0.05) t_r10 = a19_sim_p1_schedule(sch10, 15000, 0.05) cr2 = a19_crossing(sch2); cr5 = a19_crossing(sch5); cr10 = a19_crossing(sch10) # injection profile inj_t = [i*0.02 for i in range(0, 130)] inj_rate = [a19_q_inj(x % A19_T_INJ_PERIOD) for x in inj_t] # E(p) & rho(p) curves ep_p = list(range(0, 201, 10)); ep_E = [a19_E_of(x) for x in ep_p] rp_p = list(range(0, 201, 10)); rp_rho = [a19_rho_of(x) for x in rp_p] # ---- Problem 2 ---- omega2 = a19_find_omega(100.0, 28, 36, 3000, 0.02) t2 = a19_sim_p2(omega2, 3000, 0.02) # cam / plunger profile cam_phi = [i*0.05 for i in range(0, 126)] # 0..6.3 cam_L = [a19_L_of(ph) for ph in cam_phi] # omega sensitivity omega_sweep = [] for w in [20,24,28,30,32,32.66,34,36,40,45,50]: _, ps, _ = a19_sim_p2(w, 3000, 0.02); omega_sweep.append((w, round(a19_mean_p(ps),3))) # ---- Problem 3: two nozzles + relief valve ---- # twin nozzles double the drain; the relief valve caps the PEAK at 100 MPa, # so the achievable mean asymptotes just below 100 (trough injection-limited # at ~95.2 MPa). Target a mean of 98.5 MPa (comfortably inside 100+-2). omega3 = a19_find_omega(98.5, 55, 100, 3000, 0.02, relief={"p_high":100.0}, nozzles=2) t3 = a19_sim_p2(omega3, 3000, 0.02, relief={"p_high":100.0}, nozzles=2) t3_norelief = a19_sim_p2(omega3, 3000, 0.02, relief=None, nozzles=2) omega_sweep3 = [] for w in [50,55,60,65,68,70,75,80]: _, ps, _ = a19_sim_p2(w, 3000, 0.02, relief={"p_high":100.0}, nozzles=2) omega_sweep3.append((w, round(a19_mean_p(ps),3))) # ---- write CSV (key results table) ---- rows = [ ["参数", "数值", "单位"], ["高压油管内腔长度 L", A19_L_TUBE, "mm"], ["高压油管内径 d", A19_D_TUBE, "mm"], ["高压油管容积 V", round(A19_V_TUBE,2), "mm^3"], ["入口A小孔直径", 1.4, "mm"], ["入口小孔面积 A_in", round(A19_A_IN,4), "mm^2"], ["高压油泵供油压力", A19_P_SUPPLY, "MPa"], ["初始压力 p0", A19_P0, "MPa"], ["流量系数 C", A19_C, "-"], ["单向阀关闭时间", A19_T_CLOSE, "ms"], ["喷油频率", 10, "Hz"], ["喷油峰值速率", A19_INJ_PEAK, "mm^3/ms"], ["柱塞腔直径", A19_D_PLUNGER, "mm"], ["柱塞腔残余容积", A19_V_RES, "mm^3"], ["低压燃油压力", A19_P_LOW, "MPa"], ["100MPa处密度", A19_RHO_REF, "mg/mm^3"], ["100MPa处弹性模量", round(a19_E_of(100),1), "MPa"], ["凸轮升程幅值 H", round(A19_CAM_H,4), "mm"], ["问题1 稳定100MPa 开阀时长", round(topen100,4), "ms"], ["问题1 稳定150MPa 开阀时长", round(topen150,4), "ms"], ["问题1 100->150 自然调整时间", round(T_nat,1), "ms"], ["问题1 2s工况 充电开阀时长", round(tc2,4), "ms"], ["问题1 2s工况 到达150时间", round(cr2,1), "ms"], ["问题1 5s工况 到达150时间", round(cr5,1), "ms"], ["问题1 10s工况 到达150时间", round(cr10,1), "ms"], ["问题2 凸轮角速度(稳100MPa)", round(omega2,3), "rad/s"], ["问题2 凸轮角速度(折合转速)", round(omega2*60/(2*math.pi),2), "rpm"], ["问题3 双喷嘴凸轮角速度", round(omega3,3), "rad/s"], ["问题3 减压阀开启阈值", 100.0, "MPa"], ] with open(os.path.join(OUT, out_csv), "w", newline="") as f: w = _csv.writer(f); w.writerows(rows) return dict( topen100=topen100, topen150=topen150, T_nat=T_nat, t1_100=t1_100, t1_150=t1_150, t1_wrong=t1_wrong, topen_sweep=topen_sweep, sch2=sch2, sch5=sch5, sch10=sch10, t_r2=t_r2, t_r5=t_r5, t_r10=t_r10, cr2=cr2, cr5=cr5, cr10=cr10, tc2=tc2, inj_t=inj_t, inj_rate=inj_rate, ep_p=ep_p, ep_E=ep_E, rp_p=rp_p, rp_rho=rp_rho, omega2=omega2, t2=t2, cam_phi=cam_phi, cam_L=cam_L, omega_sweep=omega_sweep, omega3=omega3, t3=t3, t3_norelief=t3_norelief, omega_sweep3=omega_sweep3, V_tube=A19_V_TUBE, A_in=A19_A_IN, P_supply=A19_P_SUPPLY, P0=A19_P0, cam_H=A19_CAM_H, cam_RMIN=A19_CAM_RMIN, E100=a19_E_of(100), rho100=A19_RHO_REF, csv=out_csv, ) if __name__ == "__main__": D = gen_2019a() print("topen100=%.4f topen150=%.4f T_nat=%.1f ms" % (D["topen100"], D["topen150"], D["T_nat"])) print("omega2=%.3f rad/s (%.2f rpm) omega3=%.3f rad/s" % (D["omega2"], D["omega2"]*60/(2*math.pi), D["omega3"])) print("ramp crossings: 2s=%.1f 5s=%.1f 10s=%.1f ms" % (D["cr2"], D["cr5"], D["cr10"])) print("cam H=%.4f mm E(100)=%.1f MPa V_tube=%.1f mm^3" % (D["cam_H"], D["E100"], D["V_tube"])) print("CSV ->", D["csv"])