Abstract
With advancements in computer-aided design, simulation of internal combustion engines has become a vital tool for product development and design innovation. Among the simulation software packages currently available, MATLAB/Simulink is widely used for automotive system simulations, but does not contain a comprehensive engine modeling toolbox. To leverage MATLAB/Simulink’s capabilities, a Simulink-based 1D flow engine modeling framework has been developed. The framework allows engine component blocks to be connected in a physically representative manner in the Simulink environment, reducing model build time. Each component block, derived from physical laws, interacts with other blocks according to block connection. In this Part 1 of series papers, a comprehensive gas dynamics model is presented and integrated in the engine modeling framework based on MATLAB/Simulink. Then, the gas dynamics model is validated with commercial engine simulation software by conducting a simple 1D flow simulation.
1. Introduction
Simulation-based system design, optimization, and controller development have been an integral part of automotive research. Among many areas that have benefitted from numerical simulation, engine and drivetrain simulation has received a great deal of attention. Engine simulation allows designers to predict performance gains resulting from changes in engine geometries or control strategies. Using the simulation result, designs can be optimized for fuel economy, power, and emissions without collecting extensive experimental data. With the steady advancement in the related technology and available computing power, the impact of engine simulation will increase as well.
Depending on the desired accuracy and model complexity, engine models can be classified into several different categories: lookup table-based, mean-value, 1D physics-based, and 3D CFD models. While simple lookup table-based models rely on extensive test data, mean-value models combine the overall effect of engine flow and combustion phenomena [1–3]. Because such models require little physical detail, some of the parameters must be derived from experimental testing, and some characteristics of the engine performance cannot be accurately determined. As the opposite extremum, multidimensional Computational Fluid Dynamics (CFD) can be employed to simulate engine flow and combustion [4–7]. Multidimensional CFD models require detailed geometric parameters, which in turn provide detailed performance information. However, this highly predictive approach comes at the cost of long simulation time. Therefore, multidimensional CFD models cannot be efficiently applied to simulating multiple engine cycles.
As a compromise between simulation accuracy and computational time, a reduced dimension modeling approach (1D manifold flow and 0D cylinder) is used in many commercial software packages [8, 9]. However, most of the current 1D physics-based engine simulation models are not compatible to or directly implemented in today’s predominant MATLAB/Simulink environment. In an effort to develop a MATLAB/Simulink-based engine modeling and simulation framework, a high-fidelity 1D physics-based engine simulation model has been developed by the authors based on the work of Blair [10, 11].
This paper is an extension of the authors’ previous work by including comprehensive component models for the 1D flow. From a modeling standpoint, the 1D flow model predicts mass transfer into and out of the cylinder during intake and exhaust, while the cylinder combustion model predicts piston force and engine torque. In this Part 1 of series papers, component blocks for engine breathing simulation are presented including restrictions, adjoined pipes, and boundary conditions. Then, a simple 1D pipe system consisting of boundary conditions, flow sections, and an abrupt change in flow area is simulated, and the results are compared with those of the commercial software GT-POWER.
2. Quasi-One-Dimensional Unsteady Flow
Gas flow associated with engine intake and exhaust systems is unsteady; internal energy, density, pressure, and velocity vary with time. Although flow within each duct is best described in three dimensions, the nature of internal flow restricts gases to flow primarily in the axial direction of the duct. Therefore, to reduce model complexity, flow states are defined along a single dimension resulting in a 1D model. By including a variable cross-sectional area, the model can be considered quasi-1D.
2.1. Conservation Laws
To predict flow behavior through an engine duct, rate of changes in flow states need to be determined by conservation laws [9]. For the control volume in Figure 1, flow velocity , density , specific internal energy , pressure , and area change over the differential length . Area is a fixed function of , and all flow states are a function of time and . Wall shear accounts for friction losses and is the wall heat flux.
The quasi-1D differential form of the continuity equation based on the conservation of mass, can be derived as
Also, by applying Newton’s second law of motion, conservation of momentum can be formulated as where is the characteristic diameter. Neglecting shear work and using specific enthalpy for the sum of internal energy and flow work, conservation of energy leads to
2.2. Spatial Discretization
The nonlinear partial differential equations derived from the conservation of mass, momentum, and energy cannot be solved analytically. The equations can be converted to ordinary differential equations by replacing the infinitesimal length with a finite length (finite difference) and integrating with a proper ODE solver (e.g., Runge-Kutta and Euler method). To discretize the flow duct, a staggered grid approach is utilized. A staggered grid approach divides the pipe or duct into sections with an equal length as shown in Figure 2. At each cell center ), conservation of mass and energy laws determine the rates of change in density and specific internal energy , which can be used to determine cell pressure and temperature . At each cell boundary (), conservation of momentum determines the mass flow rate crossing each cell boundary, and energy flow rate can be derived using upstream cell information. The staggered grid approach was chosen over a collocated method such as the Lax-Wendroff method to improve stability and simplify Simulink block communication at the boundaries [12, 13].
The conservation of mass equation shown in Equation (1) can be converted from the differential form by substituting for (finite difference form) and using . The rate change in cell density becomes
Following the same procedure as the conservation of mass, the rate of change of energy at each cell center can be derived from Equation (3). In the finite difference form, conservation of energy becomes where is the cell’s wall surface area.
Conservation of momentum equation directly governs the boundary mass flow rate based on adjacent cell pressures, momentums, and minor losses. For generality, all pressure losses are represented by a single pressure loss coefficient . The total loss coefficient defined as includes friction shear losses , pipe bend losses , and any other minor losses . In this simulation model, the shear losses Cf is determined by the explicit Haaland correlation [14], and is calculated following Miller’s approach [15]. By rearranging Equation (6) using cell information, the pressure loss across the cell can be written as where the term accounts for the flow direction and describes pressure losses between two cells. The momentum Equation (2) can now be converted into a finite difference form by replacing differential terms with the finite length and replacing the shear loss term with the pressure loss relationship in Equation (7) for neighboring cells. For a staggered grid, the boundary mass flow rate is defined as
2.3. Numerical Integration
The conservation laws produce three ODE equations, Equations (4), (5), and (8), that can be solved numerically with an ODE solver. Starting from a specified initial condition, flow variables are updated every time step . To ensure numerical stability with explicit integration, the step size must be selected according to flow properties and discretization length . According to the Courant-Friedrichs-Lewy (CFL) condition, stability is related to the propagation velocity, time step size, and discretization length. Assuming a first-order accurate explicit ODE solver, the CFL condition states that the system will be stable if the following condition is met: where C is the CFL number, is the flow velocity, and is the speed of sound. For a higher order ODE solver, the solution can be stable with , but the relationship between stability and , , , and remains. The speed of sound relates to the stiffness of the gas and thus relates to the mass fractions and temperature. For an ideal gas, acoustic velocity is defined as
2.4. Heat Transfer
Conduction, convection, and radiative heat transfer contribute to the overall heat transfer. Among these, convection heat transfer plays an important role in accurately simulating engine flow characteristics. Heat transferred to the air entering the cylinder decreases the air density and therefore the amount of oxygen available for combustion. On the exhaust side, a significant amount of energy is transferred from the exhaust gasses to the exhaust valves, runners, and manifold. The significant energy transfer lowers the flow temperature and therefore the acoustic wave velocity.
Assuming that the heat transfer from the wall to the gas is positive, the forced-convection heat transfer relates to the wall and gas temperatures according to the relationship where is the convection heat transfer coefficient. The heat transfer coefficient can be determined from a correlation fit to the dimensionless Reynolds , Nusselt (), and Prandtl () numbers. and are defined as where and are the gas specific heat and thermal conductivity evaluated at the mean gas temperature, respectively. For laminar flow , the Nusselt number is constant for circular pipes [16]:
For turbulent internal flow , several correlations with varying accuracy and complexity exist. To maintain computational efficiency, the Colburn analogy is used [17, 18]:
3. Pressure Wave Motion
Although the mass and energy flow rate between neighboring cells can be determined with the momentum equation and cell relationships, boundary conditions and flow restrictions require another modeling approach. Acoustic pressure waves, which are small amplitude pressure expansions or contractions, travel in the 1D engine duct, and the combination of the incoming and outgoing waves dictates flow characteristics. The conservation laws derived previously capture the superposition effect of pressure wave propagation. For an abrupt change in flow area or at a flow boundary, however, mass flow must be determined from the incoming wave amplitude. The incoming boundary pressure wave amplitude can be extracted from cell states. Based on the incoming wave, boundary conditions, and geometry, the reflected wave can be derived from conservation laws. The incoming and reflected acoustic waves then dictate the boundary mass flow rate.
Blair developed a method based on acoustic wave propagation [19]. First predicted by Earnshaw, the amplitude of an acoustic pressure wave relates a fluid’s particle velocity [20]. Starting at a reference velocity , pressure , and acoustic velocity , Earnshaw showed that
The acoustic velocity, , is governed by a fluid’s stiffness and density, and for an ideal gas, the reference acoustic velocity can be defined as where is the ideal gas constant and is the reference temperature. To represent the acoustic wave in Equation (16), Blair defined a pressure amplitude ratio as
By assuming and substituting Equation (18) into Equation (16), Earnshaw’s theory becomes
The pressure wave propagates at the acoustic velocity relative to the gas particle velocity. Therefore, in reference to a fixed coordinate, the propagation velocity is the sum of the acoustic and particle velocities.
As shown in Figure 3, two pressure amplitude ratios are present in the 1D pipe model: a leftward and rightward traveling pressure amplitude ratio. The waves propagate in opposite directions according to the acoustic and particle velocities. By superimposing the two waves, a superposition pressure amplitude ratio relates to the flow states. According to Equation (18), the superposition pressure is defined as
The superposition acoustic velocity and temperature can be derived as assuming that the state changes from the reference conditions and to the superposition conditions and to be isentropic. Using Earnshaw’s theory in the form of Equation (19) and the relationship in Equation (22), it can be shown that
For the model, the reference pressure will be assumed constant and equal to the ambient absolute pressure. The reference temperature can fluctuate based on nonisentropic flow behavior.
Since the introduction of Blair’s method, numerous references have outlined boundary condition, engine valve, and flow junction models [21–23]. Additionally, the model has been thoroughly validated [21–24]. Therefore, Blair’s method is implemented at the boundaries in this study.
4. Flow Restrictions and Adjoined Pipes
The conservation of momentum equation for 1D flow Equation (8) does not consider cell boundary flow restrictions or area discontinuities. A flow restriction can be used to model an orifice, throttle valve, or any other 1D restriction not described by a loss coefficient. Flow area discontinuities can be formed by adjoined two pipes that do not have the same cross-section. Because of the abrupt change in area at a cell boundary due to restriction or adjoined pipe, pressure waves entering either side of the boundary get reflected. The amplitudes of the reflected waves are dictated by the change in area and conservation laws and, in turn, govern the cell boundary mass flow rate.
4.1. Model Setup and Conservation Laws
With the staggered grid approach, two pipes collinearly joined form a common cell boundary that may include a pipe area discontinuity and/or a flow restriction as depicted in Figure 4. Left and right cell information is determined by the 1D flow model described earlier. Therefore, mass and energy flow rates, , , , and , must be determined based on the connection geometry and adjoining cell information. By definition, the cell boundary is not a volume, and the rate of change in density and energy between stations 1 and 2 becomes zero. According to the conservation of mass, the mass flow rates across the boundary can be equated: where is the mass flow rate through the restriction throat. The flow must contract to pass through the throat area , and for real gas flow, a discharge coefficient is typically introduced to model the contraction and velocity losses. With the discharge coefficient, the effective throat area can be defined as
From Equations (25) and (26), the following can be concluded:
Similarly, the conservation of energy leads to
Parameters at station 1, station 2, and the throat are not explicitly available and must be calculated from cell information using acoustic wave theory.
4.2. Boundary Parameter Relationships
Referring to Figure 4, thermodynamic state variables and velocities at station 1, the throat, and station 2 do not hold a direct relationship to the left and right cell states but relate to incoming acoustic waves. States defined at the cell centers can provide the incoming pressure waves and and cell reference temperatures and . From Equation (20), the superposed pressure amplitude ratios for the left and right cells are defined as where specific heat ratios and are calculated by each respective cell’s temperature and mass fractions. From Equation (21), the reference temperature for the left and right cells are defined as
Now, the superposed pressure amplitude ratios defined in Equation (31) can to be split into opposite traveling acoustic waves. Referring to Figure 4, the incoming wave can be determined by extrapolating the rightward traveling wave from the left cell center to station 1 using the velocity relationship defined in Equation (24); hence,
Similarly, the leftward-traveling wave can be determined by extrapolating the right cell wave to station 2, giving
With the incoming pressure waves known, the reflected waves and and reference temperatures are determined by conservation laws and flow characteristics. First, each state and velocity must be expressed in a convenient form, and to express boundary states, thermodynamic properties must be evaluated at station temperatures and mass fractions. The temperatures at the boundary stations can be expressed as where is the superposed pressure amplitude ratio at the throat.
Similar issues arise when evaluating the boundary specific heat ratio . Therefore, is evaluated at a mean mass fraction and the mean temperature defined as
How properties are evaluated will become more apparent when discussing the overall solution method.
Using the temperature and property information, defining the remaining boundary states becomes straightforward. According to Equation (20), the pressure relationships can be derived:
Density at each station can be derived from Equations (35) and (37) and the ideal gas law, giving where is the ideal gas constant evaluated with . Velocities at stations 1 and 2 vary based on the incoming and reflected pressure amplitude ratios. According to Equation (24), velocities at stations 1 and 2 can be calculated as
The throat variables are expressed in terms of the superposed pressure amplitude ratio because neither the rightward nor leftward traveling waves are known. Therefore, must be solved iteratively and does not require a relationship similar to Equation (39).
4.3. Solution Overview
The boundary state variables defined in terms of reference temperatures and pressure amplitude ratios, when substituted into conservation laws, form a set of constraint equations. The equations relate adjoining pipe cell reference temperatures and incoming waves to the boundary reference temperatures and reflected waves. Examining the conservation laws and the relationships presented in the previous section, , , , , , , and are unknown. Therefore, solving for the unknown variables requires seven constraints. Conservation of mass and energy provide four constraints, while the remaining constraints are derived from entropy and momentum relationships. For a 1D flow model, Benson suggested modeling flow through a sudden change in area as an isentropic process [25]. Based on experience, Blair claims the assumption to be accurate for only certain situations [22]. Therefore, to ensure accuracy for all configurations, the more complete non-isentropic model proposed by Blair is used.
Flow must contract in order to pass through the junction throat. Gas contraction does not create flow separation or significant turbulence, and therefore, flow contraction is assumed isentropic. Likewise, flow from the left cell to station 1 is assumed isentropic for forward flow , and flow from the right cell to station 2 is assumed isentropic for reverse flow . By definition, the reference temperature remains constant for an isentropic process, thus providing two constraints for the reference temperatures:
Flow exiting the throat expands to the downstream cross-section area, giving rise to particle flow separation and turbulent vortices. The flow separation implies a nonisentropic process, and another relationship must be used for calculating the downstream reference temperature. Using conservation of momentum, flow information at the throat and station 2 can be related for forward flow, and throat and station 1 can be related for reverse flow. The downstream momentum equation is given by
The relationships in Equation (40) provide direct solutions to two of the unknown variables, reducing the number of unknown variables to five. Referring to Figure 5, the circled unknown variables must be obtained from the conservation mass equations, Equations (27) and (28); conservation of energy equations, Equations (29) and (30); and the momentum equation, Equation (41). After substituting velocities, densities, pressures, and temperatures expressed in terms of acoustic waves, the five nonlinear constraint equations cannot be solved analytically but must be solved iteratively. Based on experience, Blair found the Newton-Raphson method to be stable, accurate, and fast for solving the equations [22]. At the start of simulation, unknown variables are approximated based on cell initial conditions, and for subsequent time steps, the initial iterative guesses are taken from the previous time step.
(a)
(b)
The particle velocity at the throat cannot exceed the local acoustic velocity. However, the flow analysis discussed previously does not restrict , and depending on conditions, can be found to reach or exceed the throat acoustic velocity. Therefore, a new relationship must be introduced for the velocity limit known as choked or critical flow. For chocked flow, can be equated to the local acoustic velocity :
Velocity at other stations is assumed subsonic. The choked flow constraint must replace one of previously defined equations when solving for the five unknowns. The isentropic contraction assumption is still valid, mass and energy must be conserved between stations 1 and 2, and the momentum equation in Equation (41) must be retained to account for pressure recovery. Therefore, the intermediate energy equation, Equation (30), is chosen to be replaced. For chocked flow, the unknown variables shown in Figure 5 are solved from Equations (27), (28), (29), (41), and (42).
According to the model equations, four possible situations can potentially be encountered: subsonic forward, subsonic reverse, choked forward, and choked reverse flows. Each situation requires five equations to solve for five unknowns. Before finding the unknown variables, the equations must be expressed in terms of the incoming pressure amplitude ratios, cell reference temperatures, thermodynamic properties, flow areas, throat discharge coefficient, and unknown variables. Substituting density and velocity in terms of acoustic variables, Equations (38) and (39), into the conservation of mass equations, Equations (27) and (28), produces the following:
For conservation of energy, enthalpy must be calculated with the adjoined pipe mass fractions and the local temperature. The conservation of energy equations in Equations (29) and (30) can be expressed as using the acoustic relationships for temperature and velocity. The conservation of momentum equation defined in Equation (41) can be split into two constraints depending on flow direction. For forward flow, Equation (41) gives and for reverse flow, Equations (41) gives
For chocked flow, the throat acoustic velocity can be expressed in terms of throat reference temperature and pressure amplitude ratio , providing the relationship
Finally, the limit for choked flow, Equation (42), gives
With the isentropic relationships given in Equation (40), the unknown variables can be solved iteratively with Equation (43) through Equation (50) according to the flow condition (subsonic forward, subsonic reverse, choked forward, and choked reverse). For each flow situation, the velocity range, unknown variables, directly applied constraints, and iteration equation numbers are summarized in Table 1. Note that the throat velocity dictates the solution method but is also an unknown parameter. The equations for a given flow condition are continuously differentiable and can be solved using a Newton-type solver. However, when considering the solution as a whole, equations are not continuously differentiable at , , or . As a result, each flow condition is evaluated separately. The previously determined dictates the solution method, and after iterating, the next solution method is determined by the updated . Therefore, the overall solver can alternate between subsonic forward, subsonic reverse, choked forward, and choked reverse equations without encountering a derivative discontinuity.
4.4. Boundary Mass and Energy Flow Rates
After calculating the unknown variables listed in Table 1, boundary mass and energy flow rates, and , are evaluated from the acoustic wave relationships. With and as the left and right cell boundary flow rates, the rate of change of momentum between the left and right cells is neglected, breaking the form of the staggered grid approach. To account for momentum changes, a new boundary mass flow rate is introduced for the left and right cell boundary mass flow rate. Referring to Figure 6, conservation of momentum can be applied between the right and left cell centers to determine the rate of change of mass flow through the boundary . Because of the discontinuity at the boundary, momentum conservation must be applied between consecutive stations. Derived similarly to Equation (8), conservation of momentum from the left cell center to station 1 and station 2 to the right cell center is given by
Combining Equations (51) and (52), the rate of change of mass flow rate through the boundary is determined by
Note that it is not implied that , but , , , , and calculated from the acoustic wave relationships provide a way to determine the changes in momentum due to the area discontinuity. The boundary energy flow rate can be determined based on using the fact that energy is conserved between stations 1 and 2. For either flow direction, the boundary energy flow rate becomes
where enthalpy is calculated with and .
5. Boundary Condition Element
With the 1D staggered grid, the rate of change of mass flow rate cannot be determined by the momentum equation, Equation (8), at a pipe boundary—the interface between a 1D cell and a 0D volume. Therefore, mass and energy flow must be established based on external conditions and pipe boundary geometry. In some cases, the mass and energy flow rates can be explicitly defined. However, engine pipes most often connect to engine cylinders or ambient conditions, where mass and energy flow are not explicitly available. A 0D boundary (e.g., ambient boundary and engine cylinder) does not have flow velocity and is typically defined by a pressure, temperature, and flow area. The flow area can be fixed to represent the interface between a pipe and ambient conditions or vary to represent a poppet valve.
5.1. Model Setup and Conservation Laws
Determining mass and energy flow rates at the interface between a 1D duct cell and ambient conditions or a control volume (e.g., engine cylinder, crankcase, and tank) requires special considerations. The velocity of ambient conditions or a control volume can be best described in three dimensions. However, ambient velocity is typically assumed zero because a control volume is considered sufficiently large, which means that flow into the volume has very little influence on the volume particle velocity. As a result, ambient conditions and large volumes are modeled as 0D. A 0D volume does not have a velocity field and can be defined by mass fractions and two thermodynamic state variables.
Referring to Figure 7, station 1 represents a 0D volume, and cell parameters, denoted with subscript “C,” are governed by the 1D flow model discussed previously. Flow from the 0D volume to the 1D cell is assumed positive, but depending on pipe flow convention, signs of flow rates can be switched without loss of generality. Using the cell and boundary information, flow rate and energy flow rate are determined from thermodynamic constraints. Mass must be conserved between the throat and station 2, and by introducing a discharge coefficient CD to represent the effective throat area, mass conservation gives
Energy must also be conserved between the throat and station 2:
For flow into the pipe, conservation of energy states that the control volume enthalpy is equivalent to the total energy per unit mass at station 2, which can be described by
The throat and station 2 parameters are not explicitly available and must be calculated based on incoming and reflected acoustic waves. With the conservation relationships expressed in terms of pressure amplitude ratios and reference temperatures, the boundary mass and energy flow rates can be evaluated.
5.2. Boundary Parameter Relationships
The thermodynamic state variables at station 2 do not hold a direct relationship with the 0D volume or 1D cell. Instead, station 2 parameters relate to the incident acoustic wave derived from cell information and the reflected wave shown in Figure 7. The incident pressure amplitude ratio is derived from the cell information. According to Equation (20), the superposed pressure amplitude ratio is defined as where is the specific heat ratio calculated with the cell’s temperature and mass fractions. The incident wave relates to the cell reference temperature defined as
The superposed pressure amplitude can be split into two oppositely traveling acoustic waves based on cell velocity. Referring to Figure 7, the incident wave is determined by extrapolating the leftward traveling wave from the cell center to station 2 using the velocity relationship defined in Equation (24); hence,
Depending on flow direction, the 0D volume pressure amplitude ratio and reference temperature are required for applying constraints. From the acoustic wave relationships, the following can be concluded: where is the specific heat ratio calculated with the 0D volume’s temperature and mass fractions.
With the incoming wave and reference temperatures known, the reflected pressure amplitude ratio can be determined by conservation laws and flow characteristics. First, each state and velocity must be expressed in a convenient form, and to express boundary states, thermodynamic properties must be evaluated at station temperatures and mass fractions. The temperatures at station 2 and the throat are expressed as where is the superposed pressure amplitude ratio at the throat.
Similar issues arise when evaluating the boundary specific heat ratio . Therefore, is evaluated at , and the mean temperature is defined as
Before solving for unknown parameters, boundary states must be expressed in terms of pressure amplitude ratios and reference temperatures. According to Equation (20), the throat and station 2 pressures are defined as
Density at each station can be derived from Equations (63) and (65) and the ideal gas law as where is the ideal gas constant evaluated with . The velocity at station 2, , is determined by the incident and reflected pressure amplitude ratios. Assuming the inflow to be positive,
The throat variables are expressed in terms of the superposed pressure amplitude ratio because neither the rightward nor leftward traveling waves are known. Therefore, must be solved iteratively.
5.3. Pipe Inflow Constraints
Pipe boundary inflow and outflow require distinctly different solution approaches and therefore are discussed in separate sections. Pipe inflow, flow from 0D volume into a pipe, is assumed positive, i.e., . Similar to the adjoined pipe solution, pressure amplitude ratios and reference temperatures at the boundary are unknown and must be solved for using isentropic relationships and conservation laws. Examining the conservation laws and relationships presented in the previous section reveals that , , , , and are unknown. As a result, the boundary solution requires five constraints. Some constraints result in a direct solution to specific variables, while the remaining constraint equations must be solved iteratively. For inflow, conservation of mass and energy equations provide three constrains. Isentropic assumptions and conservation of momentum provide the remaining relationships.
During inflow, the gas must contract to pass through the boundary throat area . The contraction does not create turbulence and can be considered an isentropic process. According to the definition of the reference temperature, the 0D volume and throat reference temperate can be equated; thus, for inflow,
Exiting the throat, the gas expands to the area , giving rise to particle flow separation and turbulent vortices. The flow separation implies a nonisentropic process, and another relationship must be used to determine the downstream reference temperature, . Using conservation of momentum, flow information at the throat and station 2 can be related. The downstream momentum equation can be expressed as
For subsonic inflow, the isentropic contraction assumption expressed in Equation (68) gives a direct solution to . Referring to Figure 8(a), the remaining variables, , , , and , are determined by simultaneously solving the equations of conservation of mass, Equation (55); energy, Equations (56) and (57); and momentum, Equation (69). Before solving, equations must be expressed in terms of acoustic variables and boundary information. Substituting acoustic wave variables into Equation (55), conservation of mass provides the constraint
(a)
(b)
The conservation of energy equations contain terms for enthalpy, and from preliminary testing, calculating enthalpy at the local temperature has convergence issues. To provide stability, enthalpy changes are assumed to have constant slope, implying a constant specific heat . For an ideal gas, the change in enthalpy between the throat and station 2 can be expressed as
Now, the conservation of energy given in Equation (56) becomes
Similarly, the conservation of energy from the 0D volume to the throat in Equation (57) becomes
The final constraint defined in Equation (69), conservation of momentum, can be expressed as
Using the constraint equations, the throat velocity may be found to reach or exceed the local acoustic velocity depending on boundary conditions and throat flow area. However, the particle velocity at the throat cannot exceed the local acoustic velocity. Therefore, new relationships must be derived for the choked flow. The ratio between the throat and 0D volume pressures for the chocked flow, known as the critical pressure ratio, can be derived from the conservation of energy equations. The critical pressure ratio is derived as
By substituting pressure relationships into Equation (75), the throat pressure amplitude ratio can be directly determined by
Now, the particle velocity at the throat can be calculated directly by
Choked inflow allows and to be calculated directly using Equations (76) and (77). Visualized in Figure 8(b), the remaining unknown parameters, and , are evaluated by solving the equations of conservation of mass, Equation (55), and energy, Equation (57), assuming velocity at station 2 is subsonic.
5.4. Pipe Outflow Constraints
Pipe outflow, defined by , requires the same number of constraints as the pipe inflow. The unknown variables, , , , , and , are determined from conservation laws and isentropic relationships. As before, flow contraction is assumed isentropic. Thus, for pipe outflow, the throat and station 2 reference temperatures can be defined as
As flow exits the throat, the gas expands into the open space of the volume, creating significant turbulence. The dissipation of energy due to turbulence has traditionally been assumed to not produce pressure recovery. With no pressure recovery from the throat to the volume implying , the throat pressure amplitude ratio is defined as
The assumptions in Equations (78) and (79) provide direct solutions to , , and . Shown in Figure 9, the remaining unknown variables, and , can be determined from the conservation of mass and energy equations. Like pipe inflow, mass and energy must be conserved from station 2 to the boundary throat. Substituting acoustic wave variables into Equation (55), conservation of mass provides the constraint
(a)
(b)
Using the principles discussed for pipe inflow, the conservation of energy given in Equation (56) can be expressed as
When solving for and with the constraint equations, the throat velocity may be found to reach or exceed the local acoustic velocity. New relationships must be introduced for the choked flow. The isentropic relationships defined in Equation (78) are still valid for choked outflow but velocity is restricted by
With the constraint, the pressure recovery assumption defined in Equation (79) can be dropped for the choked flow, and Equation (82) is substituted into the conservation of mass and energy equations. Substituting Equation (82) into Equation (80), conservation of mass for the choked outflow becomes
Substituting Equation (82) into Equation (81), conservation of energy for the choked outflow then becomes
By defining directly and dropping the pressure recovery constraint, and become unknown variables which can be solved with Equations (83) and (84) assuming velocity at station 2 is subsonic.
5.5. Solution Overview
When determining the mass and energy flow rates at the interface between a 0D volume and a 1D cell, four distinct situations can be encountered: subsonic inflow, choked inflow, subsonic outflow, and choked outflow. Each situation requires different solution approaches, and because boundary state variables are not available, constraints must be formed using acoustic wave information. After expressing boundary variables in terms of reference temperatures and pressure amplitude ratios, four unknown variables must be determined: , , , and . Some constraints provide direct solutions to some of the unknowns, while the remaining variables must be solved iteratively. After solving for the unknown variables, the mass and energy flow rates at the boundary can be determined. Velocity range, unknown variables, directly applied constraints, and iteration equation numbers are summarized in Table 2.
As mentioned previously, the constraint equations listed in Table 2 cannot be reduced and must be solved iteratively. Blair found the Newton-Raphson method to be stable, accurate, and fast for solving boundary constraints [22]. To implement a Newton-type solver, equations must be continuously differentiable. For a given flow condition, equations meet the requirement. When considering the solution as a whole, however, equations are not continuously differentiable at , , or . As a result, each flow condition is evaluated separately. The previously determined dictates the solution method, and after each iteration, the next solution method is determined by the updated . Therefore, the overall solver can alternate between subsonic inflow, choked inflow, subsonic outflow, and choked outflow equations without encountering a derivative discontinuity.
5.6. Boundary Mass and Energy Flow Rates
After calculating the unknown variables listed in Table 1, boundary mass and energy flow rates, and , can be calculated with , , and . Note that an abrupt change in volume pressure or temperature inevitably causes an abrupt change in boundary mass and energy flow rates and, as a result, causes stability issues. To prevent numerical instability, the rate of change of boundary mass flow rate shown in Figure 10 must be regulated. The rate of change of can be determined by applying the conservation of momentum from the cell center to station 2. However, the resulting formulation would include half the cell length, and according to the CFL condition defined in Equation (9), the stable time step would be halved as a result. Additionally, the formulation does not follow the central difference scheme utilized by the staggered grid approach. Considering the issues with applying the momentum equation, the rate of change of is regulated by a time constant . Using as the target mass flow rate, the rate change in is defined as
The time constant can be selected based on the simulation time step or the CFL condition.
6. Model Validation
Comparing the proposed Simulink architecture and models to similar commercial software provides understanding of accuracy and usability. A simple 1D pipe is simulated, and the results are compared to those of the commercial engine simulation software GT-POWER. Similar to the proposed model, GT-POWER allows users to connect engine components in a physically representative manner, and engine performance can be predicted based on user inputs. Unlike GT-POWER, however, the proposed model allows users to completely customize flow and combustion models for research purposes. The model setup and results discussed in this section demonstrate the performance similarities of the two flow models.
6.1. Model Setup
In general, a 1D pipe system consists of boundary conditions, flow sections, and a possible abrupt change in flow area. To compare the proposed Simulink model to GT-POWER, the pipe system shown in Figure 11 is simulated. Temperature at both boundaries are held at 300 K, and the inlet pressure starts at 1 bar and increases to a steady 1.1 bar after 0.001 s, while the outlet pressure remains at 1 bar. Due to the increasing inlet pressure, flow enters the 25 mm pipe, and at the pipe exit, the gas must restrict to pass through the 20 mm pipe, creating a pressure drop at the pipe interface. The outlet boundary has a 15 mm orifice to represent a boundary restriction loss. Both the GT-POWER and Simulink models assume adiabatic flow, surface roughness of mm, and thermodynamic properties of dry air. Initially, the flow in the pipes is at rest with pressure at 1 bar and temperature at 300 K.
As shown in Figure 12(a), the GT-POWER represents each type of component with blocks: “EndEnvironment,” “OrificeConn,” and “PipeRound.” The blocks handle data logging and inlet boundary pressure internally. The developed model shown in Figure 12(b) represents the pipe in a similar manner within the Simulink environment. However, mass flow rates are measured by optional mass flow rate sensors, and inlet pressure is provided by an external Simulink block. The mass flow rates are logged by the Simulink “Scope,” and inlet pressure is provided by a source block. In general, outputs from the engine model blocks can be connected to any Simulink block, and inputs can be provided by any traditional block. To accept the inputs and provide outputs, the S-function contained in the “Engine Model Control” block remotely communicates to the “Inlet Boundary” and mass flow rate sensors.
(a)
(b)
6.2. Results and Discussion
The simulation results obtained by the Simulink model closely match those provided by GT-POWER as shown in Figure 13. At the start of simulation, the increasing pressure at the inlet boundary causes a gradual rise in inlet mass flow rate, and after reaching a steady boundary pressure, mass flow becomes steady until a pressure wave reflects back to the boundary. As expected, the outlet flow rate does not increase until the initial acoustic wave reaches the 15 mm orifice at the exit. The step changes in mass flow rates during unsteady flow are a result of the initial pressure wave propagating and reflecting at the pipe interface, inlet boundary, and outlet restriction. Note that the Simulink and GT-POWER produce nearly identical results at the start of simulation, but as time progresses, the wave front produced by the Simulink model tends to lag behind the GT-POWER model due to a difference in wave propagation velocity. The difference in wave velocity can be attributed to minor differences in thermodynamic properties or model assumptions. Variation in the steady flow rates are likely a result of differences in friction factor models.
7. Conclusions
As an extension of a previous work to develop a Simulink-based modeling framework for internal combustion engines, this paper presents a comprehensive set of component models including restrictions, adjoined pipes, boundary conditions, and junctions. The framework allows engine component blocks to be connected in a physically representative manner in the Simulink environment, which significantly reduces model build time.
To examine the performance of the modeling and simulation framework, a simple 1D pipe model was created using the developed model components and simulated in Simulink. Simulation results show that the Simulink-based 1D gas dynamics model can produce comparable results to those produced by commercially available software, GT-POWER. In Part 2 of this series papers, an engine combustion model will be presented together with crank dynamics.
Data Availability
This paper describes developing an engine modeling and simulation tool in MATLAB/Simulink and presents a basic simulation result compared with commercial software. Therefore, there is no meaningful data to share produced in the research.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This article has been adapted from a Ph.D. dissertation of the first author, Bradley Thompson.