# Internal Combustion Engine Modeling Framework in Simulink: Gas Dynamics Modeling

Authors: Bradley Thompson; Hwan-Sik Yoon
Journal: Modelling and Simulation in Engineering
 (2020)

---

## 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.


---

## Body

## 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 U, density ρ, specific internal energy e, pressure p, and area A change over the differential length dx. Area is a fixed function of x, and all flow states are a function of time and x. Wall shear τw accounts for friction losses and q̇w is the wall heat flux.Figure 1
Control volume for compressible, unsteady flow.The quasi-1D differential form of the continuity equation based on the conservation of mass, can be derived as(1)∂∂tρA+∂∂xρAU=0.Also, by applying Newton’s second law of motion, conservation of momentum can be formulated as(2)∂∂tρUA+∂∂xρU2A=−A∂p∂x−τwπD.where D=4A/π1/2 is the characteristic diameter. Neglecting shear work and using specific enthalpy h for the sum of internal energy and flow work, conservation of energy leads to
(3)∂∂tρAe+U22+∂∂xρUAh+U22=q̇wπD.


### 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Δx (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 Δx as shown in Figure 2. At each cell center (i=1,2,⋯n), conservation of mass and energy laws determine the rates of change in density ρi and specific internal energy ei, which can be used to determine cell pressure pi and temperature Ti. At each cell boundary (i=1/2,3/2,⋯n+1/2), conservation of momentum determines the mass flow rate ṁ crossing each cell boundary, and energy flow rate Ė 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].Figure 2
Staggered grid discretization.The conservation of mass equation shown in Equation (1) can be converted from the differential form by substituting Δx for ∂x (finite difference form) and using ṁ=ρUA. The rate change in cell density becomes
(4)dρidt=1AiΔxṁi−1/2−ṁi+1/2.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
(5)ddtρiei+Ui22=1AiΔxĖi−1/2−Ėi+1/2+q̇w,iAsurf,i,where Asurf,i is the ith 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 coefficientCloss. The total loss coefficient defined as
(6)Closs=Cf+Cbend+Cother=ploss1/2ρU2,includes friction shear losses Cf, pipe bend losses Cbend, and any other minor losses Cother. In this simulation model, the shear losses Cf is determined by the explicit Haaland correlation [14], and Cbend is calculated following Miller’s approach [15]. By rearranging Equation (6) using cell information, the pressure loss across the ith cell ploss,i can be written as
(7)ploss,i=Closs,i12ρiUiUi,where the term UiUi accounts for the flow direction and Closs,i 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 Δx and replacing the shear loss term with the pressure loss relationship in Equation (7) for neighboring cells. For a staggered grid, the i+1/2th boundary mass flow rate is defined as
(8)dṁi+1/2dt=1ΔxṁiUi−ṁi+1Ui+1+Ai+1/2pi−pi+1−12Aiploss,i+Ai+1ploss,i+1.


### 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 Δt. To ensure numerical stability with explicit integration, the step size must be selected according to flow properties and discretization length Δx. 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:
(9)C=U+aΔtΔx<1,where C is the CFL number, U is the flow velocity, and a is the speed of sound. For a higher order ODE solver, the solution can be stable with C>1, but the relationship between stability and U, a, Δt, and Δx 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
(10)a=γy,Tpρ=γy,TRT.


### 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(11)q̇w=hcTw−T,where hc is the convection heat transfer coefficient. The heat transfer coefficient can be determined from a correlation fit to the dimensionless Reynolds Re, Nusselt (Nu), and Prandtl (Pr) numbers. Nu and Pr are defined as
(12)Nu=hcDk,(13)Pr=cpμk,where cp and k are the gas specific heat and thermal conductivity evaluated at the mean gas temperature, respectively. For laminar flow Re<2300, the Nusselt number is constant for circular pipes [16]:
(14)Nu=3.66.For turbulent internal flowRe>4000, several correlations with varying accuracy and complexity exist. To maintain computational efficiency, the Colburn analogy is used [17, 18]:
(15)Nu=f8RePr1/3.


## 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 U, density ρ, specific internal energy e, pressure p, and area A change over the differential length dx. Area is a fixed function of x, and all flow states are a function of time and x. Wall shear τw accounts for friction losses and q̇w is the wall heat flux.Figure 1
Control volume for compressible, unsteady flow.The quasi-1D differential form of the continuity equation based on the conservation of mass, can be derived as(1)∂∂tρA+∂∂xρAU=0.Also, by applying Newton’s second law of motion, conservation of momentum can be formulated as(2)∂∂tρUA+∂∂xρU2A=−A∂p∂x−τwπD.where D=4A/π1/2 is the characteristic diameter. Neglecting shear work and using specific enthalpy h for the sum of internal energy and flow work, conservation of energy leads to
(3)∂∂tρAe+U22+∂∂xρUAh+U22=q̇wπD.


## 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Δx (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 Δx as shown in Figure 2. At each cell center (i=1,2,⋯n), conservation of mass and energy laws determine the rates of change in density ρi and specific internal energy ei, which can be used to determine cell pressure pi and temperature Ti. At each cell boundary (i=1/2,3/2,⋯n+1/2), conservation of momentum determines the mass flow rate ṁ crossing each cell boundary, and energy flow rate Ė 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].Figure 2
Staggered grid discretization.The conservation of mass equation shown in Equation (1) can be converted from the differential form by substituting Δx for ∂x (finite difference form) and using ṁ=ρUA. The rate change in cell density becomes
(4)dρidt=1AiΔxṁi−1/2−ṁi+1/2.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
(5)ddtρiei+Ui22=1AiΔxĖi−1/2−Ėi+1/2+q̇w,iAsurf,i,where Asurf,i is the ith 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 coefficientCloss. The total loss coefficient defined as
(6)Closs=Cf+Cbend+Cother=ploss1/2ρU2,includes friction shear losses Cf, pipe bend losses Cbend, and any other minor losses Cother. In this simulation model, the shear losses Cf is determined by the explicit Haaland correlation [14], and Cbend is calculated following Miller’s approach [15]. By rearranging Equation (6) using cell information, the pressure loss across the ith cell ploss,i can be written as
(7)ploss,i=Closs,i12ρiUiUi,where the term UiUi accounts for the flow direction and Closs,i 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 Δx and replacing the shear loss term with the pressure loss relationship in Equation (7) for neighboring cells. For a staggered grid, the i+1/2th boundary mass flow rate is defined as
(8)dṁi+1/2dt=1ΔxṁiUi−ṁi+1Ui+1+Ai+1/2pi−pi+1−12Aiploss,i+Ai+1ploss,i+1.


## 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 Δt. To ensure numerical stability with explicit integration, the step size must be selected according to flow properties and discretization length Δx. 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:
(9)C=U+aΔtΔx<1,where C is the CFL number, U is the flow velocity, and a is the speed of sound. For a higher order ODE solver, the solution can be stable with C>1, but the relationship between stability and U, a, Δt, and Δx 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
(10)a=γy,Tpρ=γy,TRT.


## 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(11)q̇w=hcTw−T,where hc is the convection heat transfer coefficient. The heat transfer coefficient can be determined from a correlation fit to the dimensionless Reynolds Re, Nusselt (Nu), and Prandtl (Pr) numbers. Nu and Pr are defined as
(12)Nu=hcDk,(13)Pr=cpμk,where cp and k are the gas specific heat and thermal conductivity evaluated at the mean gas temperature, respectively. For laminar flow Re<2300, the Nusselt number is constant for circular pipes [16]:
(14)Nu=3.66.For turbulent internal flowRe>4000, several correlations with varying accuracy and complexity exist. To maintain computational efficiency, the Colburn analogy is used [17, 18]:
(15)Nu=f8RePr1/3.


## 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 U [20]. Starting at a reference velocity U0, pressure p0, and acoustic velocity a0, Earnshaw showed that
(16)U−U0=2γ−1a0pp0γ−1/2γ−1.The acoustic velocity,a, is governed by a fluid’s stiffness and density, and for an ideal gas, the reference acoustic velocity can be defined as
(17)a0=γRT0,where R is the ideal gas constant and T0 is the reference temperature. To represent the acoustic wave in Equation (16), Blair defined a pressure amplitude ratio X as
(18)X=pp0γ−1/2γ.By assumingU0=0 and substituting Equation (18) into Equation (16), Earnshaw’s theory becomes
(19)U=2γ−1a0X−1.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 Figure3, two pressure amplitude ratios are present in the 1D pipe model: a leftward XL and rightward XR 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 XS relates to the flow states. According to Equation (18), the superposition pressure pS is defined as
(20)pS=p0XS2γγ−1.Figure 3
Acoustic waves traveling in pipe.The superposition acoustic velocityaS and temperature TS can be derived as
(21)TS=T0XS2,(22)aS=a0XS,assuming that the state changes from the reference conditions p0 and T0 to the superposition conditions pS and TS to be isentropic. Using Earnshaw’s theory in the form of Equation (19) and the relationship in Equation (22), it can be shown that
(23)XS=XR+XL−1,(24)Us=2γ−1a0XR−XL.For the model, the reference pressurep0 will be assumed constant and equal to the ambient absolute pressure. The reference temperature T0 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 Figure4. Left and right cell information is determined by the 1D flow model described earlier. Therefore, mass and energy flow rates, ṁ1, ṁ2, Ė1, and Ė2, 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:
(25)ṁ1=ṁ2=ṁt,where ṁt is the mass flow rate through the restriction throat. The flow must contract to pass through the throat area At, and for real gas flow, a discharge coefficient CD is typically introduced to model the contraction and velocity losses. With the discharge coefficient, the effective throat area Ateff can be defined as
(26)Ateff=CDAt.Figure 4
Schematic of adjoined pipes with restrictive orifice.From Equations (25) and (26), the following can be concluded:
(27)ṁ1−ṁ2=ρ1U1A1−ρ2U2A2=0,(28)ṁ1−ṁt=ρ1U1A1−ρtUtCDAt=0.Similarly, the conservation of energy leads to(29)h1+U122−h2+U222=0,(30)h1+U122−ht+Ut22=0.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 Figure4, 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 XR1 and XL2 and cell reference temperatures T0L and T0R. From Equation (20), the superposed pressure amplitude ratios for the left XSL and right XSR cells are defined as
(31)XSL=pLp0γL−1/2γL,XSR=pRp0γR−1/2γR,where specific heat ratios γL and γR are calculated by each respective cell’s temperature and mass fractions. From Equation (21), the reference temperature for the left T0L and right T0R cells are defined as
(32)T0L=TLXSL2,T0R=TRXSR2.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 XR1 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,
(33)XR1=12ULγL−12γLT0LRL+XSL+1.Similarly, the leftward-traveling waveXL2 can be determined by extrapolating the right cell wave to station 2, giving
(34)XL2=12−URγR−12γRT0RRR+XSR+1.With the incoming pressure waves known, the reflected wavesXL1 and XR2 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
(35)T1=T01XR1+XL1−12,T2=T02XR2+XL2−12,Tt=T0tXt2,where Xt is the superposed pressure amplitude ratio at the throat.Similar issues arise when evaluating the boundary specific heat ratioγa. Therefore, γa is evaluated at a mean mass fraction ya and the mean temperature Ta defined as
(36)Ta=TL+TR2.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:
(37)p1=p0XR1+XL1−12γa/γa−1,p2=p0XR2+XL2−12γa/γa−1,pt=p0Xt2γa/γa−1.Density at each station can be derived from Equations (35) and (37) and the ideal gas law, giving
(38)ρ1=p0RaT01XR1+XL1−12/γa−1,ρ2=p0RaT02XR2+XL2−12/γa−1,ρt=p0RaT0tXt2/γa−1,where Ra is the ideal gas constant evaluated with ya. 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
(39)U1=2γa−1γaT01RaXR1−XL1,U2=2γa−1γaT02RaXR2−XL2.The throat variables are expressed in terms of the superposed pressure amplitude ratioXt because neither the rightward nor leftward traveling waves are known. Therefore, Ut 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,XL1, XR2, Xt, T01, T0t, T02, and Ut 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 flowUt>0, and flow from the right cell to station 2 is assumed isentropic for reverse flow Ut<0. By definition, the reference temperature remains constant for an isentropic process, thus providing two constraints for the reference temperatures:
(40)ifUt≥0T0t=T01=T0L,elseT0t=T02=T0R.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(41)ifUt≥0A2pt−p2+ṁ1Ut−U2=0,elseA1pt−p1+ṁ1Ut−U1=0.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.Figure 5
(a) Forward and (b) reverse parameter constraints for adjoined pipe boundary (unknowns are circled).
(a)(b)The particle velocity at the throatUt cannot exceed the local acoustic velocity. However, the flow analysis discussed previously does not restrict Ut, and depending on conditions, Ut 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, Ut can be equated to the local acoustic velocity at:
(42)Ut=at.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:
(43)A1T01γaT01RaXR1−XL1XR1+XL1−12/γa−1−A2T02γaT02RaXR2−XL2XR2+XL2−12/γa−1=0,(44)2A1T01γa−1γaT01RaXR1−XL1XR1+XL1−12/γa−1−CDAtT0tUtXt2/γa−1=0.For conservation of energy, enthalpy must be calculated with the adjoined pipe mass fractionsya and the local temperature. The conservation of energy equations in Equations (29) and (30) can be expressed as
(45)h1ya,T1=T01XR1+XL1−12−h2ya,T2=T02XR2+XL2−12+2γaRaγa−12T01XR1−XL12−T02XR2−XL22=0,(46)h1ya,T1=T01XR1+XL1−12−htya,Tt=T0tXt2+2γaT01Raγa−12XR1−XL12−Ut22=0,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
(47)A2Xt2γa/γa−1−XR2+XL2−12γa/γa−1+2A1γaT01RaRaT01γa−1XR1−XL1XR1+XL1−12/γa−1×Ut−2γa−1γaT02RaXR2−XL2=0,and for reverse flow, Equations (41) gives
(48)A1Xt2γa/γa−1−XR1+XL1−12γa/γa−1+2A1γaT01RaRaT01γa−1XR1−XL1XR1+XL1−12/γa−1×Ut−2γa−1γaT01RaXR1−XL1=0.For chocked flow, the throat acoustic velocity can be expressed in terms of throat reference temperatureT0t and pressure amplitude ratio Xt, providing the relationship
(49)at−XtγaT0tRa=0.Finally, the limit for choked flow, Equation (42), gives
(50)Ut−XtγaT0tRa=0.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 Ut 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 Ut=0, Ut=at, or Ut=−at. As a result, each flow condition is evaluated separately. The previously determined Ut dictates the solution method, and after iterating, the next solution method is determined by the updated Ut. Therefore, the overall solver can alternate between subsonic forward, subsonic reverse, choked forward, and choked reverse equations without encountering a derivative discontinuity.Table 1
Solution summary for adjoined pipes.
Subsonic forwardSubsonic reverseChoked forwardChoked reverseVelocity range0≤Ut<at−at<Ut<0Ut=atUt=−atUnknownsXL1,XR2,Ut,Xt,andT02XL1,XR2,Ut,Xt,andT01XL1,XR2,Ut,Xt,andT02XL1,XR2,Ut,Xt,andT01Direct constraintsT01=T0LT0t=T0LT02=T0RT0t=T0RT01=T0LT0t=T0LT02=T0RT0t=T0RConstraint equations43, 44, 45, 46, and 4743, 44, 45, 46, and 4843, 44, 45, 47, and 5043, 44, 45, 48, and 50


### 4.4. Boundary Mass and Energy Flow Rates

After calculating the unknown variables listed in Table1, boundary mass and energy flow rates, ṁ1 and Ė1, are evaluated from the acoustic wave relationships. With ṁ1 and Ė1 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 ṁa 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 ṁa. 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
(51)dṁadt=2ΔxLṁLUL−ṁ1U1+A1pL−p1−12ALploss,L,(52)dṁadt=2ΔxRṁ2U2−ṁRUR+A2p2−pR−12ARploss,R.Figure 6
Mass and energy flow rate across adjoined pipe boundary.Combining Equations (51) and (52), the rate of change of mass flow rate through the boundary is determined by
(53)dṁadt=2ΔxL+ΔxR×ṁLUL−ṁRUR+ṁ1U2−U1+A1pL−p1+A2p2−pR−12ALploss,L+ARploss,R.Note that it is not implied thatṁa=ṁ1, but ṁ1, p1, p2, U1, and U2 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 Ėa can be determined based on ṁa using the fact that energy is conserved between stations 1 and 2. For either flow direction, the boundary energy flow rate becomes
(54)Ėa=ṁah1+U122,where enthalpyh1 is calculated with ya and T1.


## 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 Figure4. Left and right cell information is determined by the 1D flow model described earlier. Therefore, mass and energy flow rates, ṁ1, ṁ2, Ė1, and Ė2, 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:
(25)ṁ1=ṁ2=ṁt,where ṁt is the mass flow rate through the restriction throat. The flow must contract to pass through the throat area At, and for real gas flow, a discharge coefficient CD is typically introduced to model the contraction and velocity losses. With the discharge coefficient, the effective throat area Ateff can be defined as
(26)Ateff=CDAt.Figure 4
Schematic of adjoined pipes with restrictive orifice.From Equations (25) and (26), the following can be concluded:
(27)ṁ1−ṁ2=ρ1U1A1−ρ2U2A2=0,(28)ṁ1−ṁt=ρ1U1A1−ρtUtCDAt=0.Similarly, the conservation of energy leads to(29)h1+U122−h2+U222=0,(30)h1+U122−ht+Ut22=0.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 Figure4, 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 XR1 and XL2 and cell reference temperatures T0L and T0R. From Equation (20), the superposed pressure amplitude ratios for the left XSL and right XSR cells are defined as
(31)XSL=pLp0γL−1/2γL,XSR=pRp0γR−1/2γR,where specific heat ratios γL and γR are calculated by each respective cell’s temperature and mass fractions. From Equation (21), the reference temperature for the left T0L and right T0R cells are defined as
(32)T0L=TLXSL2,T0R=TRXSR2.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 XR1 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,
(33)XR1=12ULγL−12γLT0LRL+XSL+1.Similarly, the leftward-traveling waveXL2 can be determined by extrapolating the right cell wave to station 2, giving
(34)XL2=12−URγR−12γRT0RRR+XSR+1.With the incoming pressure waves known, the reflected wavesXL1 and XR2 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
(35)T1=T01XR1+XL1−12,T2=T02XR2+XL2−12,Tt=T0tXt2,where Xt is the superposed pressure amplitude ratio at the throat.Similar issues arise when evaluating the boundary specific heat ratioγa. Therefore, γa is evaluated at a mean mass fraction ya and the mean temperature Ta defined as
(36)Ta=TL+TR2.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:
(37)p1=p0XR1+XL1−12γa/γa−1,p2=p0XR2+XL2−12γa/γa−1,pt=p0Xt2γa/γa−1.Density at each station can be derived from Equations (35) and (37) and the ideal gas law, giving
(38)ρ1=p0RaT01XR1+XL1−12/γa−1,ρ2=p0RaT02XR2+XL2−12/γa−1,ρt=p0RaT0tXt2/γa−1,where Ra is the ideal gas constant evaluated with ya. 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
(39)U1=2γa−1γaT01RaXR1−XL1,U2=2γa−1γaT02RaXR2−XL2.The throat variables are expressed in terms of the superposed pressure amplitude ratioXt because neither the rightward nor leftward traveling waves are known. Therefore, Ut 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,XL1, XR2, Xt, T01, T0t, T02, and Ut 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 flowUt>0, and flow from the right cell to station 2 is assumed isentropic for reverse flow Ut<0. By definition, the reference temperature remains constant for an isentropic process, thus providing two constraints for the reference temperatures:
(40)ifUt≥0T0t=T01=T0L,elseT0t=T02=T0R.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(41)ifUt≥0A2pt−p2+ṁ1Ut−U2=0,elseA1pt−p1+ṁ1Ut−U1=0.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.Figure 5
(a) Forward and (b) reverse parameter constraints for adjoined pipe boundary (unknowns are circled).
(a)(b)The particle velocity at the throatUt cannot exceed the local acoustic velocity. However, the flow analysis discussed previously does not restrict Ut, and depending on conditions, Ut 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, Ut can be equated to the local acoustic velocity at:
(42)Ut=at.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:
(43)A1T01γaT01RaXR1−XL1XR1+XL1−12/γa−1−A2T02γaT02RaXR2−XL2XR2+XL2−12/γa−1=0,(44)2A1T01γa−1γaT01RaXR1−XL1XR1+XL1−12/γa−1−CDAtT0tUtXt2/γa−1=0.For conservation of energy, enthalpy must be calculated with the adjoined pipe mass fractionsya and the local temperature. The conservation of energy equations in Equations (29) and (30) can be expressed as
(45)h1ya,T1=T01XR1+XL1−12−h2ya,T2=T02XR2+XL2−12+2γaRaγa−12T01XR1−XL12−T02XR2−XL22=0,(46)h1ya,T1=T01XR1+XL1−12−htya,Tt=T0tXt2+2γaT01Raγa−12XR1−XL12−Ut22=0,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
(47)A2Xt2γa/γa−1−XR2+XL2−12γa/γa−1+2A1γaT01RaRaT01γa−1XR1−XL1XR1+XL1−12/γa−1×Ut−2γa−1γaT02RaXR2−XL2=0,and for reverse flow, Equations (41) gives
(48)A1Xt2γa/γa−1−XR1+XL1−12γa/γa−1+2A1γaT01RaRaT01γa−1XR1−XL1XR1+XL1−12/γa−1×Ut−2γa−1γaT01RaXR1−XL1=0.For chocked flow, the throat acoustic velocity can be expressed in terms of throat reference temperatureT0t and pressure amplitude ratio Xt, providing the relationship
(49)at−XtγaT0tRa=0.Finally, the limit for choked flow, Equation (42), gives
(50)Ut−XtγaT0tRa=0.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 Ut 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 Ut=0, Ut=at, or Ut=−at. As a result, each flow condition is evaluated separately. The previously determined Ut dictates the solution method, and after iterating, the next solution method is determined by the updated Ut. Therefore, the overall solver can alternate between subsonic forward, subsonic reverse, choked forward, and choked reverse equations without encountering a derivative discontinuity.Table 1
Solution summary for adjoined pipes.
Subsonic forwardSubsonic reverseChoked forwardChoked reverseVelocity range0≤Ut<at−at<Ut<0Ut=atUt=−atUnknownsXL1,XR2,Ut,Xt,andT02XL1,XR2,Ut,Xt,andT01XL1,XR2,Ut,Xt,andT02XL1,XR2,Ut,Xt,andT01Direct constraintsT01=T0LT0t=T0LT02=T0RT0t=T0RT01=T0LT0t=T0LT02=T0RT0t=T0RConstraint equations43, 44, 45, 46, and 4743, 44, 45, 46, and 4843, 44, 45, 47, and 5043, 44, 45, 48, and 50


## 4.4. Boundary Mass and Energy Flow Rates

After calculating the unknown variables listed in Table1, boundary mass and energy flow rates, ṁ1 and Ė1, are evaluated from the acoustic wave relationships. With ṁ1 and Ė1 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 ṁa 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 ṁa. 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
(51)dṁadt=2ΔxLṁLUL−ṁ1U1+A1pL−p1−12ALploss,L,(52)dṁadt=2ΔxRṁ2U2−ṁRUR+A2p2−pR−12ARploss,R.Figure 6
Mass and energy flow rate across adjoined pipe boundary.Combining Equations (51) and (52), the rate of change of mass flow rate through the boundary is determined by
(53)dṁadt=2ΔxL+ΔxR×ṁLUL−ṁRUR+ṁ1U2−U1+A1pL−p1+A2p2−pR−12ALploss,L+ARploss,R.Note that it is not implied thatṁa=ṁ1, but ṁ1, p1, p2, U1, and U2 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 Ėa can be determined based on ṁa using the fact that energy is conserved between stations 1 and 2. For either flow direction, the boundary energy flow rate becomes
(54)Ėa=ṁah1+U122,where enthalpyh1 is calculated with ya and T1.


## 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 Figure7, 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 ṁ2 and energy flow rate Ė2 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
(55)ṁt−ṁ2=ρtUtCDAt−ρ2U2A2=0.Figure 7
Schematic of 1D pipe boundary condition.Energy must also be conserved between the throat and station 2:(56)ht+Ut22−h2+U222=0.For flow into the pipe, conservation of energy states that the control volume enthalpyh1 is equivalent to the total energy per unit mass at station 2, which can be described by
(57)h1−h2+U222=0.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 waveXi2 derived from cell information and the reflected wave Xr2 shown in Figure 7. The incident pressure amplitude ratio Xi2 is derived from the cell information. According to Equation (20), the superposed pressure amplitude ratio XSC is defined as
(58)XSC=pcp0γc−1/2γc,where γC is the specific heat ratio calculated with the cell’s temperature and mass fractions. The incident wave relates to the cell reference temperature T0C defined as
(59)T0C=TCXSC2.The superposed pressure amplitudeXSC can be split into two oppositely traveling acoustic waves based on cell velocity. Referring to Figure 7, the incident wave Xi2 is determined by extrapolating the leftward traveling wave from the cell center to station 2 using the velocity relationship defined in Equation (24); hence,
(60)Xi2=12−UCγC−12γCT0CRC+XSC+1.Depending on flow direction, the 0D volume pressure amplitude ratioX1 and reference temperature T01 are required for applying constraints. From the acoustic wave relationships, the following can be concluded:
(61)X1=p1p0γ1−1/2γ1(62)T01=T1X12,where γ1 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 ratioXr2 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
(63)T2=T02Xr2+Xi2−12,Tt=T0tXt2,where Xt is the superposed pressure amplitude ratio at the throat.Similar issues arise when evaluating the boundary specific heat ratioγb. Therefore, γb is evaluated at yb, and the mean temperature Tb is defined as
(64)Tb=T1+TC2.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
(65)p2=p0Xi2+Xr2−12γb/γb−1,pt=p0Xt2γb/γb−1.Density at each station can be derived from Equations (63) and (65) and the ideal gas law as
(66)ρ2=p0RbT02Xr2+Xi2−12/γb−1,ρt=p0RbT0tXt2/γb−1,where Rb is the ideal gas constant evaluated with yb. The velocity at station 2, U2, is determined by the incident and reflected pressure amplitude ratios. Assuming the inflow to be positive,
(67)U2=2γb−1γbT02RbXr2−Xi2.The throat variables are expressed in terms of the superposed pressure amplitude ratioXt because neither the rightward nor leftward traveling waves are known. Therefore, Ut 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.,Ut>0. 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 Xr2, Xt, T0t, T02, and Ut 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 areaAt. 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,
(68)T0t=T01.Exiting the throat, the gas expands to the areaA2, 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, T02. Using conservation of momentum, flow information at the throat and station 2 can be related. The downstream momentum equation can be expressed as
(69)A2pt−p2+ṁ2Ut−U2=0.For subsonic inflow, the isentropic contraction assumption expressed in Equation (68) gives a direct solution to T0t. Referring to Figure 8(a), the remaining variables, Xr2, Xt, T02, and Ut, 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
(70)CDAtT0tUtXt2/γb−1−2A2T02γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1.Figure 8
Boundary constraints and unknown parameters for (a) subsonic and (b) choked inflow (unknowns are circled).
(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 heatCp. For an ideal gas, the change in enthalpy between the throat and station 2 can be expressed as
(71)ht−h2=CpTt−T2=γbRbγb−1Tt−T2=at2−a22γb−1.Now, the conservation of energy given in Equation (56) becomes
(72)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+Ut22−2γbT02Rbγb−12Xr2−Xi22=0.Similarly, the conservation of energy from the 0D volume to the throat in Equation (57) becomes
(73)γbRbγb−1T01X12−T02Xr2+Xi2−12−2γbT02Rbγb−12Xr2−Xi22=0.The final constraint defined in Equation (69), conservation of momentum, can be expressed as
(74)A2Xt2γb/γb−1−Xr2+Xi2−12γb/γb−1+2A2γbT02RbRbT02γb−1Xr2−Xi2Xr2+Xi2−12/γb−1×Ut−2γb−1γbT02RbXr2−Xi2=0.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 throatUt 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
(75)ptp1=2γb+1γb/γb−1.By substituting pressure relationships into Equation (75), the throat pressure amplitude ratio Xt can be directly determined by
(76)Xt=X12γb+1.Now, the particle velocity at the throat can be calculated directly by(77)Ut=at=XtγbT0tRb.Choked inflow allowsXt and Ut to be calculated directly using Equations (76) and (77). Visualized in Figure 8(b), the remaining unknown parameters, Xr2 and T02, 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 byUt<0, requires the same number of constraints as the pipe inflow. The unknown variables, Xr2, Xt, T0t, T02, and Ut, 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
(78)T0t=T02=T0C.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 implyingp1=pt, the throat pressure amplitude ratio Xt is defined as
(79)Xt=p1p0γb−1/2γb.The assumptions in Equations (78) and (79) provide direct solutions to T0t, T02, and Xt. Shown in Figure 9, the remaining unknown variables, Ut and Xr2, 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
(80)CDAtUtXt2/γb−1−2A2γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1=0.Figure 9
Boundary constraints and unknown parameters for (a) subsonic and (b) choked outflow (unknowns circled).
(a)(b)Using the principles discussed for pipe inflow, the conservation of energy given in Equation (56) can be expressed as
(81)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+Ut22−2γbT02Rbγb−12Xr2−Xi22=0.When solving forUt and Xr2 with the constraint equations, the throat velocity Ut 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
(82)Ut=−at=−XtγbT0tRb.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
(83)−CDAtγbT0tRbXtγb+1/γb−1−2A2γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1=0.Substituting Equation (82) into Equation (81), conservation of energy for the choked outflow then becomes
(84)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+γbT0tRbXt22−2γbT02Rbγb−12Xr2−Xi22=0.By definingUt directly and dropping the pressure recovery constraint, Xt and Xr2 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:Xr2, Xt, T02, and Ut. 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.Table 2
Solution summary for boundary conditions.
Subsonic inflowChoked inflowSubsonic outflowChoked outflowVelocity range0≤Ut<atUt=at−at<Ut<0Ut=−atUnknownsXr2, Xt, T02, and UtXr2 and T02Xr2 and UtXr2 and XtDirect constraintsT0t=T01T0t=T01Xt=X12γb+1Ut=XtγbT0tRbT02=T0CT0t=T0CXt=p1p0γb−1/2γbT02=T0CT0t=T0CUt=XtγbT0tRbConstraint equations70, 72, 73, and 7470 and 7380 and 8183 and 84As mentioned previously, the constraint equations listed in Table2 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 Ut=0, Ut=at, or Ut=−at. As a result, each flow condition is evaluated separately. The previously determined Ut dictates the solution method, and after each iteration, the next solution method is determined by the updated Ut. 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 Table1, boundary mass and energy flow rates, ṁ2 and Ė2, can be calculated with T02, Xi2, and Xr2. Note that an abrupt change in volume pressure p1 or temperature T1 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 ṁb shown in Figure 10 must be regulated. The rate of change of ṁb 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 ṁb is regulated by a time constant τb. Using ṁ2 as the target mass flow rate, the rate change in ṁb is defined as
(85)dṁbdt=1τbṁ2−ṁb.Figure 10
Mass and energy flow rate across boundary.The time constantτb can be selected based on the simulation time step or the CFL condition.


## 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 Figure7, 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 ṁ2 and energy flow rate Ė2 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
(55)ṁt−ṁ2=ρtUtCDAt−ρ2U2A2=0.Figure 7
Schematic of 1D pipe boundary condition.Energy must also be conserved between the throat and station 2:(56)ht+Ut22−h2+U222=0.For flow into the pipe, conservation of energy states that the control volume enthalpyh1 is equivalent to the total energy per unit mass at station 2, which can be described by
(57)h1−h2+U222=0.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 waveXi2 derived from cell information and the reflected wave Xr2 shown in Figure 7. The incident pressure amplitude ratio Xi2 is derived from the cell information. According to Equation (20), the superposed pressure amplitude ratio XSC is defined as
(58)XSC=pcp0γc−1/2γc,where γC is the specific heat ratio calculated with the cell’s temperature and mass fractions. The incident wave relates to the cell reference temperature T0C defined as
(59)T0C=TCXSC2.The superposed pressure amplitudeXSC can be split into two oppositely traveling acoustic waves based on cell velocity. Referring to Figure 7, the incident wave Xi2 is determined by extrapolating the leftward traveling wave from the cell center to station 2 using the velocity relationship defined in Equation (24); hence,
(60)Xi2=12−UCγC−12γCT0CRC+XSC+1.Depending on flow direction, the 0D volume pressure amplitude ratioX1 and reference temperature T01 are required for applying constraints. From the acoustic wave relationships, the following can be concluded:
(61)X1=p1p0γ1−1/2γ1(62)T01=T1X12,where γ1 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 ratioXr2 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
(63)T2=T02Xr2+Xi2−12,Tt=T0tXt2,where Xt is the superposed pressure amplitude ratio at the throat.Similar issues arise when evaluating the boundary specific heat ratioγb. Therefore, γb is evaluated at yb, and the mean temperature Tb is defined as
(64)Tb=T1+TC2.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
(65)p2=p0Xi2+Xr2−12γb/γb−1,pt=p0Xt2γb/γb−1.Density at each station can be derived from Equations (63) and (65) and the ideal gas law as
(66)ρ2=p0RbT02Xr2+Xi2−12/γb−1,ρt=p0RbT0tXt2/γb−1,where Rb is the ideal gas constant evaluated with yb. The velocity at station 2, U2, is determined by the incident and reflected pressure amplitude ratios. Assuming the inflow to be positive,
(67)U2=2γb−1γbT02RbXr2−Xi2.The throat variables are expressed in terms of the superposed pressure amplitude ratioXt because neither the rightward nor leftward traveling waves are known. Therefore, Ut 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.,Ut>0. 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 Xr2, Xt, T0t, T02, and Ut 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 areaAt. 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,
(68)T0t=T01.Exiting the throat, the gas expands to the areaA2, 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, T02. Using conservation of momentum, flow information at the throat and station 2 can be related. The downstream momentum equation can be expressed as
(69)A2pt−p2+ṁ2Ut−U2=0.For subsonic inflow, the isentropic contraction assumption expressed in Equation (68) gives a direct solution to T0t. Referring to Figure 8(a), the remaining variables, Xr2, Xt, T02, and Ut, 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
(70)CDAtT0tUtXt2/γb−1−2A2T02γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1.Figure 8
Boundary constraints and unknown parameters for (a) subsonic and (b) choked inflow (unknowns are circled).
(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 heatCp. For an ideal gas, the change in enthalpy between the throat and station 2 can be expressed as
(71)ht−h2=CpTt−T2=γbRbγb−1Tt−T2=at2−a22γb−1.Now, the conservation of energy given in Equation (56) becomes
(72)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+Ut22−2γbT02Rbγb−12Xr2−Xi22=0.Similarly, the conservation of energy from the 0D volume to the throat in Equation (57) becomes
(73)γbRbγb−1T01X12−T02Xr2+Xi2−12−2γbT02Rbγb−12Xr2−Xi22=0.The final constraint defined in Equation (69), conservation of momentum, can be expressed as
(74)A2Xt2γb/γb−1−Xr2+Xi2−12γb/γb−1+2A2γbT02RbRbT02γb−1Xr2−Xi2Xr2+Xi2−12/γb−1×Ut−2γb−1γbT02RbXr2−Xi2=0.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 throatUt 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
(75)ptp1=2γb+1γb/γb−1.By substituting pressure relationships into Equation (75), the throat pressure amplitude ratio Xt can be directly determined by
(76)Xt=X12γb+1.Now, the particle velocity at the throat can be calculated directly by(77)Ut=at=XtγbT0tRb.Choked inflow allowsXt and Ut to be calculated directly using Equations (76) and (77). Visualized in Figure 8(b), the remaining unknown parameters, Xr2 and T02, 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 byUt<0, requires the same number of constraints as the pipe inflow. The unknown variables, Xr2, Xt, T0t, T02, and Ut, 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
(78)T0t=T02=T0C.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 implyingp1=pt, the throat pressure amplitude ratio Xt is defined as
(79)Xt=p1p0γb−1/2γb.The assumptions in Equations (78) and (79) provide direct solutions to T0t, T02, and Xt. Shown in Figure 9, the remaining unknown variables, Ut and Xr2, 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
(80)CDAtUtXt2/γb−1−2A2γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1=0.Figure 9
Boundary constraints and unknown parameters for (a) subsonic and (b) choked outflow (unknowns circled).
(a)(b)Using the principles discussed for pipe inflow, the conservation of energy given in Equation (56) can be expressed as
(81)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+Ut22−2γbT02Rbγb−12Xr2−Xi22=0.When solving forUt and Xr2 with the constraint equations, the throat velocity Ut 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
(82)Ut=−at=−XtγbT0tRb.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
(83)−CDAtγbT0tRbXtγb+1/γb−1−2A2γb−1γbT02RbXr2−Xi2Xr2+Xi2−12/γb−1=0.Substituting Equation (82) into Equation (81), conservation of energy for the choked outflow then becomes
(84)γbRbγb−1T0tXt2−T02Xr2+Xi2−12+γbT0tRbXt22−2γbT02Rbγb−12Xr2−Xi22=0.By definingUt directly and dropping the pressure recovery constraint, Xt and Xr2 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:Xr2, Xt, T02, and Ut. 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.Table 2
Solution summary for boundary conditions.
Subsonic inflowChoked inflowSubsonic outflowChoked outflowVelocity range0≤Ut<atUt=at−at<Ut<0Ut=−atUnknownsXr2, Xt, T02, and UtXr2 and T02Xr2 and UtXr2 and XtDirect constraintsT0t=T01T0t=T01Xt=X12γb+1Ut=XtγbT0tRbT02=T0CT0t=T0CXt=p1p0γb−1/2γbT02=T0CT0t=T0CUt=XtγbT0tRbConstraint equations70, 72, 73, and 7470 and 7380 and 8183 and 84As mentioned previously, the constraint equations listed in Table2 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 Ut=0, Ut=at, or Ut=−at. As a result, each flow condition is evaluated separately. The previously determined Ut dictates the solution method, and after each iteration, the next solution method is determined by the updated Ut. 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 Table1, boundary mass and energy flow rates, ṁ2 and Ė2, can be calculated with T02, Xi2, and Xr2. Note that an abrupt change in volume pressure p1 or temperature T1 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 ṁb shown in Figure 10 must be regulated. The rate of change of ṁb 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 ṁb is regulated by a time constant τb. Using ṁ2 as the target mass flow rate, the rate change in ṁb is defined as
(85)dṁbdt=1τbṁ2−ṁb.Figure 10
Mass and energy flow rate across boundary.The time constantτb 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 Figure11 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 ε=0.046 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.Figure 11
Simulated pipe parameters and boundary conditions.As shown in Figure12(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.Figure 12
(a) GT-POWER and (b) Simulink block diagrams of the 1D flow model.
(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 Figure13. 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.Figure 13
Simulated inlet and outlet mass flow rates from GT-POWER and Simulink.


## 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 Figure11 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 ε=0.046 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.Figure 11
Simulated pipe parameters and boundary conditions.As shown in Figure12(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.Figure 12
(a) GT-POWER and (b) Simulink block diagrams of the 1D flow model.
(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 Figure13. 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.Figure 13
Simulated inlet and outlet mass flow rates from GT-POWER and Simulink.


## 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.



---

*Source: 6787408-2020-09-03.xml*