engine1.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
    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
        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 = V2 + (V_cylinder - V2) * ((theta - 360) / 180)
            global P3
            P3 = P_ref
        elif 360 < theta < 540:
            Vt = V2 + (V_cylinder - V2) * ((theta - 360) / 180)
            P_ref = P3 * ((V2 / Vt) ** nc)
        elif 540 <= theta <= 720:
            P_ref = 1.1 * P_environment
            Vt = V_cylinder
        P_real = (m_cyl / m_ref) * P_ref #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
    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()