engine2.py

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D





def P_total(P_environment, T_heat):
    
    #P_environment = float(input("the environment pressure is: ")) #Pa
    T_environment = 298 #K
    R = 287 #J/kg*k
    gas_density = P_environment / (R * T_environment) #kg/m^3
    
    A = 1e-3 #m^2 (general)
    u = 150 #m/s (general)
    gas_massflow_rate = gas_density * A * u #kg/s
    
    N = 15000 #rpm
    N_rps = N / 60
    N_frequency = N_rps / 120 #Hz
    z = 6  #Nums of cylinder
    m_cyl = gas_massflow_rate / (z * N_frequency) #kg (cycle trapped gas mass)
    
    P_env_ref = 1e5 #Pa (reference)
    T_ref = T_environment #K
    gas_density_ref = P_env_ref / (R * T_ref) #kg/m^3
    gas_massflow_rate = gas_density_ref * A * u
    m_ref = gas_density_ref / (z * N_frequency) #kg (reference trapped gas mass)
    
    
    #T_heat = float(input("the burning temperature is: "))
    V_cylinder = 2.66 * 1e-4 #m^3
    r = 18 #Compression ratio
    V2 = 1.48 * 1e-5    #V_cylinder / r
    nc = 1.3 #bullish index
   
    
    def e_eta(T_heat):  # Otto efficiency fomular:  eta = 1 - 1/r^(gamma-1)
         gamma = 1.4 - 0.00005 * (T_heat - 300)
         eta = 1 - (1 / np.power(r, gamma - 1))
         return eta
    def heat_loss_factor(T_ideal):
        T_wall = 470      
        k = 0.0005          # Heat loss intensity
        delta_T = max(T_ideal - T_wall, 0.0)
        factor = np.exp(-k * delta_T)
        return factor
    
    
    
    def phasing_eff(T_heat):
    # 最佳温度对应最佳相位（可调）
        T_opt = 2800  
        width = 300    # 控制峰值宽度
        eff = np.exp(-((T_heat - T_opt) / width)**2)
        return eff

    
    
    T_2 = T_environment * (r ** (nc - 1)) # the temperature before burning
    def V_theta(theta):
        if 0 <= theta < 180:
            return V_cylinder
        elif 180 <= theta < 360:
            return V_cylinder - (V_cylinder - V2) * ((theta - 180) / 180)
        elif 360 <= theta < 540:
            return V2 + (V_cylinder - V2) * ((theta - 360) / 180)
        elif 540 <= theta <= 720:
            return V_cylinder
    
    def P_ref(theta):
        if 0 <= theta < 180:
            P_ref = P_environment
            Vt = V_theta(theta)
        elif 180 <= theta < 360:
            Vt = V_theta(theta)
            P_ref = P_environment * ((V_cylinder / Vt) ** nc)
            global P2
            P2 = P_ref 
        elif theta == 360:
            PPR = 6.5 #Peak Pressure Ratio （We do not consider the Wiebe Function)
            P_ref = (T_heat / T_2) * P2 * PPR
            Vt = V_theta(theta)
            global P3
            P3 = P_ref
        elif 360 < theta < 540:
            Vt = V_theta(theta)
            P_ref = P3 * ((V2 / Vt) ** nc)
        elif 540 <= theta <= 720:
            P_ref = 1.1 * P_environment
            Vt = V_cylinder
            
            
        R_ideal_gas = 8.314
        M_air = 0.029
        T_ideal = (P_ref * Vt) / (m_cyl * (R_ideal_gas / M_air))
        
        
        P_real = (m_cyl / m_ref) * P_ref * heat_loss_factor(T_ideal) #Pa
        return P_real
    

    
    
    
    
    theta_array = np.arange(0, 721)
    V = np.array([V_theta(theta) for theta in theta_array])
    dVdtheta = np.gradient(V, theta_array)
    P = np.array([P_ref(theta) for theta in theta_array])
    W_net = z * np.trapezoid(P * dVdtheta, x=theta_array)
    
    nums_cycle_per_second = N_rps / 2
    P_power = (W_net * nums_cycle_per_second) * e_eta(T_heat) * phasing_eff(T_heat)
    return P_power
    





Pressure_env_values = np.linspace(80000, 130000, 50)   # Pa
T_heat_values = np.linspace(2000, 3500, 50)       # K
Pressure_env_grid, T_heat_grid = np.meshgrid(Pressure_env_values, T_heat_values)

Power_net_grid = np.zeros_like(Pressure_env_grid)
for i in range(Pressure_env_grid.shape[0]):
    for j in range(Pressure_env_grid.shape[1]):
        Pressure_env = Pressure_env_grid[i, j]
        T_heat = T_heat_grid[i, j]
        Power_net_grid[i, j] = P_total(Pressure_env, T_heat)

fig = plt.figure(figsize=(12, 8)) #create a figure
ax = fig.add_subplot(111, projection='3d') # create 3D

surf = ax.plot_surface(
    Pressure_env_grid,
    T_heat_grid,
    Power_net_grid,
    cmap='viridis',
    edgecolor='none'
)     # create a 3D curved surface

ax.set_xlabel('Environment Pressure P_env (Pa)', fontsize=12)
ax.set_ylabel('combustion Temperature T_heat (K)', fontsize=12)
ax.set_zlabel('Power (W)', fontsize=12)
ax.set_title('Engine Power Surface vs P_env and T_heat', fontsize=14)
# create the axies labels

#fig.colorbar(surf, shrink=0.6, aspect=10) #show the change of colors

plt.show()