JP2005330931A - Control device of engine - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a shaft torque calculation device capable of estimating a shaft torque generated by an engine without directly detecting a pressure in a combustion chamber by a pressure sensor. <P>SOLUTION: This control device of the engine comprises an operating state detection means 31 detecting the operating state of the engine, an estimating means 31 estimating, based on the detected operating state, either one of a combustion mass at crank angles during the combustion of a gas in a combustion chamber, and a combustion mass ratio at crank angles during the combustion of the gas in the combustion chamber and a heat generating pattern, and a shaft torque estimating means 31 estimating the shaft torque occurring from the engine based on either one of the estimated values. <P>COPYRIGHT: (C)2006,JPO&NCIPI

Description

  The present invention relates to an engine (internal combustion engine) shaft torque calculation device and an engine control device.

There is one that calculates the indicated torque based on the output of a pressure sensor that detects the pressure in the combustion chamber (see Patent Document 1).
JP-A-8-31426

  By the way, in the technique described in Patent Document 1, a pressure sensor for detecting the pressure in the combustion chamber is required, which is expensive.

  Therefore, the present invention provides an axial torque calculation device capable of estimating the axial torque generated by the engine without directly detecting the pressure in the combustion chamber by the pressure sensor, and an engine control device using the estimated axial torque. The purpose is to provide.

  The present invention detects the operating state of the engine, and based on the detected operating state, the combustion mass at the crank angle during combustion of the gas in the combustion chamber, the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber, One of the heat generation patterns, which is a value related to the combustion mass or the combustion mass ratio, is estimated, and the shaft torque generated by the engine is estimated based on the estimated value.

  In addition to this, the present invention is also configured to perform torque control based on the estimated shaft torque.

  According to the present invention, the engine operating state is detected, and based on the detected operating state, the combustion mass ratio (X1 ′, X2 ′) at the crank angle during combustion of the gas in the combustion chamber is estimated, and this estimation is performed. Since it is configured to estimate the shaft torque generated by the engine based on the combustion mass ratio, the shaft torque generated by the engine can be estimated without directly detecting the pressure in the combustion chamber by the pressure sensor. A device can be realized.

  Hereinafter, embodiments of the present invention will be described with reference to the drawings.

  FIG. 1 is a schematic diagram for explaining the system of the present invention.

  The air is stored in the intake collector 2 and then introduced into the combustion chamber 5 of each cylinder via the intake manifold 3. Fuel (gasoline) is injected and supplied from a fuel injector 21 arranged in the intake port 4 of each cylinder. The fuel injected into the air is vaporized and mixed with the air to form a gas (air mixture) and flows into the combustion chamber 5. This air-fuel mixture is confined in the combustion chamber 5 when the intake valve 15 is closed, and is compressed by the rise of the piston 6.

  In order to ignite this compressed air-fuel mixture with a high-pressure spark, an ignition device 11 of an electronic power distribution system is provided in which an ignition coil with a built-in power transistor is arranged in each cylinder. That is, the ignition device 11 includes an ignition coil 13 that stores electrical energy from the battery, a power transistor that energizes and shuts off the primary side of the ignition coil 13, and a primary current of the ignition coil 13 that is provided on the ceiling of the combustion chamber 5. It includes a spark plug 14 that receives a high voltage generated on the secondary side of the ignition coil 13 due to interruption of the spark coil 13 and performs spark discharge.

  When a spark is blown off by the spark plug 14 slightly before the compression top dead center and the compressed mixture is ignited, the flame spreads and then explosively burns, and the gas pressure by this combustion works to push down the piston 6. This work is taken out as the rotational force of the crankshaft 7. The combusted gas (exhaust gas) is discharged into the exhaust passage 8 when the exhaust valve 16 is opened.

A three-way catalyst 9 is provided in the exhaust passage 8. When the air-fuel ratio of the exhaust gas is in a narrow range (window) centered on the stoichiometric air-fuel ratio, the three-way catalyst 9 can efficiently remove harmful three components such as HC, CO, and NOx contained in the exhaust gas simultaneously. Since the air-fuel ratio is the ratio of the intake air amount and the fuel amount, the intake air amount introduced into the combustion chamber 5 per one cycle of the engine (crank angle 720 ° section in a four-cycle engine) and the fuel injector 21 The engine controller 31 uses the intake air flow rate signal from the air flow meter 32 and the fuel from the fuel injector 21 based on the signals from the crank angle sensors (33, 34) so that the ratio to the fuel injection amount becomes the stoichiometric air-fuel ratio. The injection amount is determined, and the air-fuel ratio is feedback controlled based on a signal from an O ₂ sensor 35 provided upstream of the three-way catalyst 9.

  A so-called electronically controlled throttle 22 in which a throttle valve 23 is driven by a throttle motor 24 is provided upstream of the intake collector 2. Since the torque required by the driver appears in the amount of depression of the accelerator pedal 41 (accelerator opening), the engine controller 31 determines a target torque based on a signal from the accelerator sensor 42, and a target air for realizing this target torque. The amount is determined, and the opening degree of the throttle valve 23 is controlled via the throttle motor 24 so as to obtain this target air amount.

  Cam sprockets and crank sprockets are attached to the front portions of the intake valve camshaft 25, the exhaust valve camshaft 26, and the crankshaft 7, respectively, and a timing chain (not shown) is hung around these sprockets so that the camshaft 25 and 26 are driven by the crankshaft 7 of the engine, and are interposed between the cam sprocket and the intake valve camshaft 25 to continuously adjust the phase of the intake valve cam with a constant operating angle. A variable intake valve timing control mechanism (hereinafter referred to as “intake VTC mechanism”) 27 that can be controlled, and a cam sprocket and an exhaust valve camshaft 26 are interposed between the camshaft 26 and the exhaust valve cam. Variable exhaust valve timing control mechanism (hereinafter referred to as “exhaust”) that can control the phase continuously. Provided with that.) 28 and TC mechanism ". When the opening / closing timing of the intake valve 15 and the opening / closing timing of the exhaust valve 16 are changed, the amount of the inert gas remaining in the combustion chamber 5 changes. As the amount of the inert gas in the combustion chamber 5 increases, the pumping loss decreases and the fuel consumption improves. Therefore, the target intake valve closing timing and the target exhaust gas indicate how much inert gas should remain in the combustion chamber 5 depending on the operating conditions. The valve closing timing is determined in advance, and the engine controller 31 determines the target intake valve closing timing and the target exhaust valve closing timing from the operating conditions (engine load and rotation speed) at that time, so that these target values can be obtained. The intake valve closing timing and the exhaust valve closing timing are controlled via the actuators of the intake VTC mechanism 27 and the exhaust VTC mechanism 28.

  An intake air temperature signal from the intake air temperature sensor 43, an intake air pressure signal from the intake air pressure sensor 44, an exhaust gas temperature signal from the exhaust air temperature sensor 45, and an exhaust gas pressure signal from the exhaust air pressure sensor 46 are output from the water temperature sensor 37. The engine controller 31 that is input together with the coolant temperature signal controls the ignition timing that is the primary current cutoff timing of the spark plug 14 via the power transistor 13.

  In this case, in the present invention, in order to estimate the shaft torque generated by the engine without directly detecting the pressure in the combustion chamber by the pressure sensor, the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber is calculated. Based on the estimated combustion mass ratio, an estimated value of shaft torque generated by the engine is calculated.

  Furthermore, in the present invention, torque control (intake air amount increase and ignition timing retardation correction) is performed based on this estimated shaft torque value.

  The torque control here will be briefly explained. The ignition is basically performed at the ignition timing at which MBT is obtained, but the ignition timing may be retarded from the ignition timing at which MBT is obtained. The shaft torque generated at the time is lower than the shaft torque when ignition is performed at the ignition timing at which MBT is obtained. Therefore, the intake air amount is increased in order to increase the shaft torque at the ignition timing retarded from the ignition timing at which MBT is obtained and to approach the shaft torque at the ignition timing at which MBT is obtained.

  On the other hand, the required shaft torque requested by the driver is determined from the accelerator opening and the engine speed, and if the shaft torque at the ignition timing obtained by MBT exceeds this required shaft torque, the shaft torque is excessively output. This can cause a sense of incongruity. Therefore, when the shaft torque at the ignition timing at which MBT is obtained exceeds the required shaft torque, the ignition timing is retarded.

  FIG. 2 is a flowchart showing the overall flow of ignition timing control. This flow shows a flow of operations, not a flow executed at regular intervals.

  In step 1, the basic ignition timing MBTCAL [degBTDC] is calculated, and in steps 2 to 7, the ignition timing retard amount RET (new) [deg] is calculated.

  Here, the calculation of the basic ignition timing MBTCAL will be described first. First, the ignition timing control based on the combustion analysis will be outlined (the basic concept is described in Japanese Patent Laid-Open No. 2003-148236).

  As shown in FIG. 3, the crank angle at which the combustion pressure of the air-fuel mixture reaches the maximum value Pmax when the air-fuel mixture is ignited with MBT (minimum advance angle value at which maximum torque can be obtained) is defined as a reference crank angle θPMAX [degATDC]. The reference crank angle θPMAX is substantially constant regardless of the combustion method, and is generally in the range of 12 to 15 degrees after compression top dead center, and at most 10 to 20 degrees after compression top dead center.

  FIG. 4 shows changes in the combustion mass ratio BR (combustion gas mass ratio) obtained by the combustion analysis in the combustion chamber in the spark ignition engine. The combustion mass ratio BR representing the ratio of the combustion mass to the fuel supplied to the combustion chamber is 0% at the time of ignition, and reaches 100% by complete combustion. Experiments have confirmed that the combustion mass ratio at the reference crank angle θPMAX is constant and about 60%.

  The combustion period corresponding to the change allowance until the combustion mass ratio BR reaches about 60% corresponding to the reference crank angle θPMAX from 0% is a period in which there is almost no change in both the combustion mass ratio and the combustion pressure immediately after the start of combustion. It is divided into an initial combustion period and a main combustion period in which the combustion mass ratio and the combustion pressure increase rapidly. The initial combustion period is a stage from the start of combustion until flame nuclei are formed, and the flame nuclei are formed when the combustion mass ratio changes from 0% to 2% to 10%. During this initial combustion period, the rate of increase in combustion pressure and combustion temperature is small, and the initial combustion period is long with respect to changes in the combustion mass ratio. The length of the initial combustion period is susceptible to changes in temperature and pressure in the combustion chamber.

  On the other hand, in the main combustion period, the flame propagates from the flame kernel to the outside, and the flame speed (that is, the combustion speed) increases rapidly. Therefore, the change in the combustion mass ratio during the main combustion period is larger than the change in the combustion mass ratio during the initial combustion period.

  In the engine controller 31, the period until the combustion mass ratio reaches 2% (changes) is set as the initial combustion period BURN1 [deg], and the period from the end of the initial combustion period BURN1 to the reference crank angle θPMAX (combustion chamber volume ratio) In other words, the main combustion period BURN2 [deg] is distinguished from 2% to about 60%. Then, a combustion period BURN [deg] that is the sum of the initial combustion period BURN1 and the main combustion period BURN2 is calculated, a reference crank angle θPMAX [degATDC] is subtracted from the combustion period BURN, and an ignition dead time equivalent crank described later is further calculated. The crank angle position to which the angle IGNDEAD [deg] is added is set as the basic ignition timing MBTCAL [degBTDC], which is the ignition timing at which MBT is obtained.

  The pressure and temperature in the combustion chamber 5 during the initial combustion period in which flame nuclei are formed are almost equivalent to the pressure and temperature at the time of ignition, but the ignition timing is calculated from this, but accurate ignition is performed from the beginning. The time cannot be set. Therefore, as shown in FIG. 13, the previous value of the basic ignition timing is calculated as the previous combustion start timing MBTCYCL [degBTDC] (step 44), and this value is used to calculate the initial combustion period as shown in FIG. (Step 162), the calculation of the initial combustion period is cyclically repeated to obtain a highly accurate result without time delay.

  Next, the calculation of the basic ignition timing MBTCAL executed by the engine controller 31 will be described in detail with reference to the following flowchart.

  FIG. 5 is for calculating various physical quantities necessary for calculating the ignition timing, and is executed at regular time intervals (for example, every 10 msec).

  First, in step 11, the intake valve closing timing IVC [degBTDC], the collector internal temperature TCOL [K] detected by the temperature sensor 43, the collector internal pressure PCOL [Pa] detected by the pressure sensor 44, and the temperature sensor 45 are detected. Exhaust temperature TEXH [K], internal inert gas rate MRESFR [%], cooling water temperature TWK [K] detected by temperature sensor 37, target equivalent ratio TFBYA, engine rotational speed NRPM [rpm] detected by crank angle sensor ].

  Here, the crank angle sensor includes a position sensor 33 for detecting the position of the crankshaft 7 and a phase sensor 34 for detecting the position of the intake camshaft 25. The engine is based on signals from the two sensors 33 and 34. The rotational speed NRPM [rpm] is calculated.

  The intake valve closing timing IVC is known from the command value given to the intake VTC mechanism 27. Alternatively, the actual intake valve closing timing may be detected by the phase sensor 34.

  The internal inert gas ratio MRESFR is a value obtained by dividing the amount of inert gas remaining in the combustion chamber by the total amount of gas in the combustion chamber, and the calculation thereof will be described later.

  The target equivalent ratio TFBYA is calculated in a fuel injection amount calculation flow (not shown). The target equivalent ratio TFBYA is an unnamed number, and is a value represented by the following expression when the theoretical air-fuel ratio is 14.7.

TFBYA = 14.7 / target air-fuel ratio (1)
For example, from equation (1), when the target air-fuel ratio is the stoichiometric air-fuel ratio, TFBYA = 1.0, and when the target air-fuel ratio is a lean value such as 22.0, TFBYA is a positive value less than 1.0. is there.

In step 12, the volume of the combustion chamber 5 at the intake valve closing timing IVC (that is, the volume at the compression start timing) VIVC [m ³ ] is calculated. The volume VIVC of the combustion chamber 5 when the intake valve is closed is determined by the stroke position of the piston 6. The stroke position of the piston 6 is determined by the crank angle position of the engine.

  Referring to FIG. 6, consider a case where the rotation center 72 of the crankshaft 71 of the engine is offset from the center axis 73 of the cylinder. Assume that the connecting rod 74, the joint point 75 between the connecting rod 74 and the crankshaft 71, and the piston pin 76 that connects the connecting rod 74 and the piston are in the relationship shown in the figure. At this time, the volume VIVC of the combustion chamber 5 at the closing timing of the intake valve can be expressed by the following equations (2) to (6).

VIVC = f1 (θivc) = Vc + (π / 4) D ² · Hivc (2)
Vc = (π / 4) D ² · Hx / (ε−1) (3)
^{Hivc = {(CND + ST 2} /2) - (CRoff-PISoff) 2} 1/2
− {(ST / 2) · cos (θivc + θoff)}
+ (CND ² −X ² ) ^1/2 (4)
X = (ST / 2) · sin (θivc + θoff) −CRoff + PISoff (5)
θoff = arcsin {(CRoff−PISoff) / (CND · (ST / 2))} (6)
Where Vc: gap volume [m ³ ],
ε: compression ratio,
D: cylinder bore diameter [m],
ST: Full piston stroke [m],
Hivc: Piston pin 76 at the intake valve closing timing
Distance from TDC [m],
Hx: The maximum and minimum values of the distance from the TDC of the piston pin 76
Difference [m],
CND: length of connecting rod 74 [m],
CRoff: Offset distance of the nodal point 75 from the cylinder center axis 73
[M],
PISoff: From the center axis 73 of the crankshaft rotation center 72
Offset distance [m],
θivc: Intake valve closing timing crank angle [degATDC],
θoff: Piston pin 76 and crankshaft rotation center 72 are connected
Angle [deg] between the line and the vertical line in TDC,
X: horizontal distance [m] between the nodal point 75 and the piston pin 76,
As described above, the crank angle θivc at the time of closing the intake valve is known because it is determined by the command signal from the engine controller 31 to the intake VTC mechanism 27. By substituting the crank angle θivc (= IVC) at this time into the equations (2) to (6), the volume VIVC of the combustion chamber 5 at the closing timing of the intake valve can be calculated. Therefore, in practice, the volume VIVC of the combustion chamber 5 at the intake valve closing timing is set by a table using the intake valve closing timing IVC as a parameter. When the intake VTC mechanism 27 is not provided, a constant value can be given.

  In step 13, the temperature (that is, the compression start timing temperature) TINI [K] of the combustion chamber 5 at the intake valve closing timing IVC is calculated. The temperature of the gas flowing into the combustion chamber 5 is a temperature of a gas in which the fresh air flowing into the combustion chamber 5 and the inert gas remaining in the combustion chamber 5 are mixed. The temperature of the fresh air flowing into the combustion chamber 5 is Since the temperature of the inert gas equal to the fresh air temperature TCOL in the intake collector 2 and remaining in the combustion chamber 5 can be approximated by the exhaust temperature TEXH in the vicinity of the exhaust port, the temperature of the combustion chamber 5 at the intake valve closing timing IVC. TINI is the following from the fresh air temperature TCOL in the intake collector 2, the exhaust gas temperature TEXH, and the internal inert gas ratio MRESFR that is the ratio of the inert gas remaining in the combustion chamber 5 at the timing when the intake valve closing timing IVC is reached. It can be obtained by an expression.

TINI = TEXH × MRESFR + TCOL × (1−MRESFR) (7)
In step 14, the pressure (that is, compression start timing pressure) PINI [Pa] at the intake valve closing timing IVC of the combustion chamber 5 is calculated. That is, the collector internal pressure PCOL at the timing when the intake valve closing timing IVC is reached is taken in as the pressure PINI at the intake valve closing timing IVC.

  In step 15, a reaction probability RPROBA [%] representing the flammability of the air-fuel mixture in the combustion chamber 5 is calculated. The reaction probability RPROBA is a dimensionless value and depends on the three parameters of the residual inert gas ratio MRESFR, the cooling water temperature TWK [K], and the target equivalent ratio TFBYA, and can be expressed by the following equation.

RPROBA = f3 (MRESFR, TWK, TFBYA) (8)
More specifically, the maximum value of the reaction probability obtained by the combination of the three parameters MRESFR, TWK, and TFBYA is set to 100%, the relationship between these parameters and the reaction probability RPROBA is experimentally obtained, and the obtained reaction probability RPROBA is obtained. Are stored in advance in the memory of the engine controller 31 as a table corresponding to the parameters. In step 14, the reaction probability RPROBA is obtained by searching this table according to the parameters.

  Specifically, a table of water temperature correction coefficients having characteristics as shown in FIG. 7 according to the cooling water temperature TWK, a table of internal inert gas rate correction coefficients (not shown) set similarly, and a target equivalent ratio A table of equivalence ratio correction coefficients having characteristics as shown in FIG. 8 according to TFBYA is stored in the memory in advance. The maximum value of each correction coefficient is 1.0, and the reaction probability RPROBA is calculated by multiplying the product of the three types of correction coefficients by the maximum value of 100% of the reaction probability.

  Explaining each table, the water temperature correction coefficient shown in FIG. 7 becomes larger as the cooling water temperature TWK is higher, and becomes 1.0 when the cooling water temperature TWK is 80 ° C. or higher. The equivalence ratio correction coefficient shown in FIG. 8 is the maximum value of 1.0 when the target equivalence ratio TFBYA is 1.0, that is, the stoichiometric air-fuel ratio. The ratio correction factor decreases. Although the internal inert gas rate correction coefficient is not shown, it is 1.0 when the internal inert gas rate MRESFR is zero.

  In step 16, a reference crank angle θPMAX [degATDC] is calculated. As described above, the reference crank angle θPMAX does not fluctuate very much, but it still tends to advance as the engine speed NRPM increases. Therefore, the reference crank angle θPMAX can be expressed by the following equation as a function of the engine speed NRPM. it can.

θPMAX = f4 (NRPM) (9)
Specifically, the reference crank angle θPMAX is obtained by searching a table of characteristics shown in FIG. 9 stored in advance in the memory of the engine controller 31 from the engine speed NRPM. In order to facilitate calculation, the reference crank angle θPMAX can be regarded as constant.

  FIG. 10 is for calculating the initial combustion period BURN1 [deg], and FIG. 12 is for calculating the main combustion period BURN2 [deg], which is executed at regular intervals (for example, every 10 msec). 10 and 12 are executed following FIG. Either of FIGS. 10 and 12 may be executed first.

First, referring to FIG. 10, in step 161, the previous combustion start timing MBTCYCL [degBTDC], the volume VIVC [m ³ ] at the intake valve closing timing of the combustion chamber 5 calculated in step 12 of FIG. 5, and the step of FIG. The temperature TINI [K] at the intake valve closing timing of the combustion chamber 5 calculated in 13, the pressure PINI [Pa] at the intake valve closing timing of the combustion chamber 5 calculated in Step 14 of FIG. 5, and the engine speed NRPM [Rpm], the reaction probability RPROBA [%] calculated in step 15 of FIG. 5 is read.

  Here, the previous combustion start timing MBTCYCL is a value one cycle before [degBTDC] of the basic ignition timing MBTCAL, and the calculation thereof will be described later with reference to FIG.

In step 162, the volume V0 [m ³ ] at the combustion start timing of the combustion chamber 5 is calculated. As described above, the ignition timing (combustion start timing) here is not the basic ignition timing MBTCAL calculated in the current cycle but a value one cycle before the basic ignition timing. That is, the volume V0 at the combustion start timing of the combustion chamber 5 is calculated from MBTCYCL, which is a value one cycle before the basic ignition timing, by the following equation.

V0 = f6 (MBTCYCL) (11)
Specifically, the volume V0 of MBTCYCL in the combustion chamber 5 is calculated from the stroke position of the piston 6 at the previous combustion start timing MBTCYCL and the bore diameter of the combustion chamber 5. In step 12 of FIG. 5, the volume VIVC of the combustion chamber 5 at the intake valve closing timing IVC is obtained by searching a table of intake valve closing timing volumes using the intake valve closing timing as a parameter. Here, MBTCYCL is set as a parameter. The volume V0 of the combustion chamber 5 at the previous combustion start time MBTCYCL may be obtained by searching the table of the previous combustion start time volume.

  In step 163, an effective compression ratio Ec at the combustion start timing is calculated. The effective compression ratio Ec is a dimensionless value, and is a value obtained by dividing the volume V0 of the combustion chamber 5 at the combustion start timing by the volume VIVC of the combustion chamber 5 at the intake valve closing timing, as shown in the following equation.

Ec = f7 (V0, VIVC)
= V0 / VIVC (12)
In step 164, the temperature increase rate TCOMP in the combustion chamber 5 from the intake valve closing timing IVC to the combustion start timing is calculated based on the effective compression ratio Ec as shown in the following equation.

TCOMP = f8 (Ec) = Ec ^ (κ−1) (13)
Where κ: specific heat ratio,
Equation (13) is an equation for the rate of temperature rise of the adiabatic compressed gas. Note that “^” on the right side of the equation (13) represents power calculation. This symbol is also used in the formula described later.

  κ is a value obtained by dividing the constant pressure specific heat of the gas adiabatically compressed by the constant volume specific heat. If the gas adiabatically compressed is air, κ = 1.4, and this value may be used simply. However, the calculation accuracy can be further improved by experimentally determining the value of κ for the air-fuel mixture.

  FIG. 11 illustrates equation (13). Therefore, it is possible to obtain the temperature increase rate TCOMP by storing a table of such characteristics in advance in the memory of the engine controller 31 and searching the table based on the effective compression ratio Ec.

In step 165, the temperature T0 [K] at the combustion start timing of the combustion chamber 5 is multiplied by the temperature TINI at the intake valve closing timing of the combustion chamber 5 by the temperature increase rate TCOMP, that is, T0 = TINI × TCOMP (14)
It is calculated by the following formula.

  Steps 166 and 167 are the same as steps 164 and 165. That is, at step 166, the pressure increase rate PCOMP in the combustion chamber 5 from the intake valve closing timing IVC to the combustion start timing is calculated based on the effective compression ratio Ec as shown in the following equation.

PCOMP = f9 (Ec) = Ec ^ κ (41)
Where κ: specific heat ratio,
The equation (41) is also an equation for the rate of increase in pressure of the gas that is adiabatically compressed as in the equation (13). “^” On the right side of the equation (41) represents power calculation as in the equation (13).

  κ is the same as the value used in the above equation (13). If the gas to be adiabatically compressed is air, κ = 1.4, and this value may be used simply. However, the calculation accuracy can be further improved by obtaining the value of κ from the composition and temperature of the air-fuel mixture.

  It is also possible to obtain a pressure increase rate PCOMP by storing a table having the same characteristics as in FIG. 11 in advance in the memory of the engine controller 31 and searching the table based on the effective compression ratio Ec.

In step 167, the pressure P0 [Pa] at the combustion start timing of the combustion chamber 5 is multiplied by the pressure PINI at the intake valve closing timing of the combustion chamber 5 by the pressure increase rate PCOMP, that is, P0 = PINI × PCOMP (42)
It is calculated by the following formula.

  In step 168, the laminar combustion speed SL1 [m / sec] in the initial combustion period is calculated by the following equation (known).

SL1 = f10 (T0, P0)
= SLstd × (T0 / Tstd) ^2.18 × (P0 / Pstd) ^−0.16 (15)
Where Tstd: reference temperature [K],
Pstd: Reference pressure [Pa]
SLstd: Reference laminar combustion at reference temperature Tstd and reference pressure Pstd
Speed [m / sec],
T0: temperature [K] at the combustion start timing of the combustion chamber 5;
P0: pressure [Pa] at the combustion start timing of the combustion chamber 5;
The laminar combustion velocity (laminar flame velocity) is the flame propagation velocity in the absence of gas flow, and is independent of the compression velocity in the combustion chamber 5 and the intake air flow velocity in the combustion chamber 5. 5 is known to be a function of the temperature and pressure of 5. Therefore, the laminar combustion speed in the initial combustion period is a function of the combustion start temperature T0 and the combustion start pressure P0, and as described later, the main combustion period. Is a function of the compression top dead center temperature TTDC and the compression top dead center pressure PTDC. This is because the laminar combustion speed generally varies depending on the engine load, the inert gas ratio in the combustion chamber 5, the intake valve closing timing, the specific heat ratio, and the intake air temperature. Since it is a factor that affects T and pressure P, it is assumed that the laminar combustion speed can be finally defined by the temperature T and pressure P in the combustion chamber 5.

  In the above equation (15), the reference temperature Tstd, the reference pressure Pstd, and the reference laminar combustion speed SLstd are values determined in advance by experiments.

Under a pressure of 2 bar or more, which is a normal pressure in the combustion chamber 5, the pressure term (P0 / Pstd) ^−0.16 in the equation (15) becomes a small value. Accordingly, it is also possible to define the reference laminar combustion speed SLstd only by the reference temperature Tstd with the pressure term (P0 / Pstd) ^−0.16 being a constant value.

  Accordingly, when the reference temperature Tstd is 550 [K], the reference laminar combustion rate SLstd is 1.0 [m / sec], and the pressure term is 0.7, the temperature T0 and the laminar combustion rate at the combustion start timing The relationship with SL1 can be approximately defined by the following equation.

SL1 = f11 (T0)
= 1.0 × 0.7 × (T0 / 550) ^2.18 (16)
In step 169, the gas flow turbulence intensity ST1 in the initial combustion period is calculated. The turbulence intensity ST1 of the gas flow is a dimensionless value, and depends on the flow rate of fresh air flowing into the combustion chamber 5 and the penetration of injected fuel from the fuel injector 21.

  The flow rate of fresh air flowing into the combustion chamber 5 depends on the shape of the intake passage, the operating state of the intake valve 15, and the shape of the intake port 4 where the intake valve 15 is provided. The penetration of the injected fuel depends on the injection pressure of the fuel injector 21, the fuel injection period, and the combustion injection timing.

  Finally, the turbulence strength ST1 of the gas flow during the initial combustion period can be expressed by the following equation as a function of the engine speed NRPM.

ST1 = f12 (NRPM) = C1 × NRPM (17)
Where C1: constant,
It is also possible to obtain the turbulence intensity ST1 from a table using the rotational speed NRPM as a parameter.

  In step 170, the gas combustion speed FLAME1 [m / sec] in the initial combustion period is calculated from the laminar combustion speed S1 and the turbulence intensity ST1 by the following equation.

FLAME1 = SL1 × ST1 (18)
If there is gas turbulence in the combustion chamber 5, the gas combustion speed changes. Equation (18) takes into account the contribution (influence) to the combustion speed associated with this gas turbulence.

  In step 171, the initial combustion period BURN1 [deg] is calculated by the following equation.

BURN1 = {(NRPM × 6) × BR1 × V0}
/ (RPROBA × AF1 × FLAME1) (19)
However, AF1: Reaction area (fixed value) of flame kernel [m ² ],
The equation (19) and the equation (22) to be described later are derived from the following basic equation that the combustion period is obtained by dividing the mass of the combustion gas by the combustion speed. The equations (19) and (22) The numerator and denominator on the right-hand side do not immediately indicate the mass of combustion gas or the burning rate.

Combustion period [sec] = total mass in cylinder [g] / (unburned gas density [g / m ³ ]
× Flame surface area [m ² ] × Flame speed [m / sec])
... (Supplement 1)
Since the unburned gas density in the right side denominator of (Supplement 1) is a value obtained by dividing the unburned gas mass [g] by the unburned gas volume [m ³ ], the charging efficiency ITAC corresponding to the mass as in the conventional device. It cannot be said that the unburned gas density can be calculated accurately with the function of only. Therefore, experimental formulas shown in the above formula (19) and formula (22) to be described later are obtained for the first time by introducing a predetermined approximation to the formula (complement 1) while checking the experimental results.

  Here, BR1 on the right side of the equation (19) is a change amount of the combustion mass ratio from the combustion start timing to the end timing of the initial combustion period BURN1, and here BR1 = 2% is set. (NRPM × 6) on the right side of the equation (19) is a process for converting the unit from rpm to crank angle (deg). The reaction area AF1 of the flame kernel is set experimentally.

  Further, it can be assumed that the combustion chamber volume does not change during the initial combustion period. Therefore, when calculating the initial combustion period BURN1, the combustion chamber volume V0 at the start of combustion, which is the first combustion chamber volume, is employed.

Next, the flow of FIG. 12 is followed. In step 181, similarly to step 161 of FIG. 10, the volume VIVC [m ³ ] at the intake valve closing timing of the combustion chamber 5 calculated in step 12 of FIG. The temperature TINI [K] at the intake valve closing timing of the combustion chamber 5 calculated in step 13 of FIG. 5, the pressure PINI [Pa] at the intake valve closing timing of the combustion chamber 5 calculated in step 14 of FIG. The speed NRPM [rpm], the reaction probability RPROBA [%] calculated in step 15 of FIG. 5 is read, and the cylinder fresh air amount MACYL [g], the target equivalent ratio TFBYA, the internal inert gas amount MRES [g], The external inert gas amount MEGR [g] is read.

  Here, although an external EGR device is not shown in FIG. 1, as far as FIG. 12 is concerned, an explanation will be given on the premise of an engine equipped with an external EGR device. In this case, the external inert gas amount MEGR may be calculated using, for example, a known method (see Japanese Patent Laid-Open No. 10-141150). It should be noted that when an engine that does not include an external EGR device as in the present embodiment shown in FIG. 1 is used, it is sufficient to handle the external inert gas amount MEGR = 0. The calculation of the cylinder fresh air amount MACYL and the internal inert gas amount MRES will be described later with reference to FIG.

Steps 182 and 183 are the same as steps 163 and 164 in FIG. That is, in step 182, the effective compression ratio Ec at the compression top dead center time. 2 is calculated. Effective compression ratio Ec 2 is also a dimensionless value like the effective compression ratio Ec of the above equation (12), and the volume VTDC at the time of compression top dead center of the combustion chamber 5 is calculated at the closing timing of the intake valve of the combustion chamber 5 as shown in the following equation. The value divided by the volume VIVC.

Ec 2 = f13 (VTDC, VIVC) = VTDC / VIVC
... (43)
In equation (43), the volume VTDC at the time of compression top dead center of the combustion chamber 5 is constant regardless of the operating conditions, and may be stored in the memory of the engine controller 31 in advance.

In step 183, the temperature increase rate TCOMP due to adiabatic compression in the combustion chamber 5 during the period from the intake valve closing timing IVC to the compression top dead center. 2 is an effective compression ratio Ec as shown in the following equation: 2 is calculated.

TCOMP 2 = f14 (Ec 2)
= Ec 2 ^ (κ-1) (44)
Where κ: specific heat ratio,
A table having the same characteristics as in FIG. 11 is stored in advance in the memory of the engine controller 31 and the effective compression ratio Ec. By searching the table from 2, the temperature rise rate TCOMP 2 can also be obtained.

  In step 184, the total gas mass MGAS [g] in the combustion chamber 5 is calculated from the cylinder fresh air amount MACYL, the target equivalent ratio TFBYA, the internal inert gas amount MRES, and the external inert gas amount MEGR by the following equation.

MGAS = MACYL × (1 + TFBYA / 14.7) + MRES + MEGR
... (45)
“1” in parentheses on the right side of the equation (45) is a fresh air, and “TFBYA / 14.7” is a fuel.

  In step 185, the total gas mass MGAS of the combustion chamber 5, the cylinder fresh air amount MACYL, and the target equivalent ratio TFBYA are used to calculate the temperature increase amount (combustion increase temperature) TBURN [K] due to the combustion of the air-fuel mixture using the following equation. .

TBURN = {MACYL × TFBYA / 14.7 × BRk × Q}
/ (Cv × MGAS) (46)
Where Q is the constant calorific value of the fuel,
BRk: Combustion mass ratio of fuel in cylinder,
Cv: constant volume specific heat,
The numerator on the right side of the equation (46) means the total heat generated by the fuel in the cylinder [J], and the denominator means the temperature increase rate per unit generated heat [J / K]. That is, the equation (46) is an approximate equation applied to the thermodynamic formula.

  Here, the combustion mass ratio BRk of the in-cylinder fuel is adapted in advance through experiments or the like. For example, 60% / 2 = 30% is set. In this embodiment, since the combustion mass ratio BR reaches about 60% as the combustion period, the intermediate 30% is set as BRk.

  Since the constant calorific value Q of fuel varies depending on the type of fuel, it is obtained in advance by experiments or the like according to the type of fuel. The constant volume specific heat Cv is a value of 2 to 3, and the representative value is adapted beforehand by an experiment or the like. However, the calculation accuracy can be further improved by obtaining the value of the constant volume specific heat Cv from the composition and temperature of the air-fuel mixture.

In step 186, the temperature TTDC [K] at the compression top dead center of the combustion chamber 5 is changed from the temperature TINI at the intake valve closing timing of the combustion chamber 5 to the temperature increase rate TCOMP to the compression top dead center. Multiply by 2 and add the combustion rise temperature TBURN to the multiplication value, that is, the following equation is used.

TTDC = TINI × TCOMP 2 + TBURN
... (47)
In step 187, the temperature TTDC and the volume VTDC at the compression top dead center of the combustion chamber 5 and the pressure PINI, the volume VIVC and the temperature TINI at the closing timing of the intake valve of the combustion chamber 5 are The pressure PTDC [K] is calculated.

PTDC = PINI × VIVC × TTDC / (VTDC × TINI)
... (48)
Equation (48) is obtained using the equation of state. That is, the following equation of state is established using the pressure, volume and temperature (PINI, VIVC, TINI) at the intake valve closing timing.

PINI x VIVC = n · R · TINI (Supplement 2)
Where n is the number of moles
R: gas constant,
Since the volume is almost equal in the vicinity of the compression top dead center, the following equation of state is established using the pressure, volume and temperature (PTDC, VTDC, TTDC) at the compression top dead center.

PTDC × VTDC = n · R · TTDC (Supplement 3)
When n · R is eliminated from both (complement 3) and (complement 2) and PTPT is solved, the above equation (48) is obtained.

  In step 188, similarly to step 168 of FIG. 10, a laminar combustion speed SL2 [m / sec] in the main combustion period is calculated by the following equation (known).

SL2 = f15 (TTDC, PTDC)
= SLstd × (TTDC / Tstd) ^2.18 × (PTDC / Pstd) ^−0.16
... (49)
Where Tstd: reference temperature [K],
Pstd: Reference pressure [Pa]
SLstd: reference laminar burning velocity at reference temperature Tstd and reference pressure Pstd
[M / sec],
TTDC: temperature [K] at the compression top dead center of the combustion chamber 5;
PTDC: pressure [Pa] at the compression top dead center of the combustion chamber 5
The explanation of the equation (49) is the same as the equation (16). That is, the reference temperature Tstd, the reference pressure Pstd, and the reference laminar combustion speed SLstd in the equation (49) are values determined in advance by experiments. Under a pressure of 2 bar or more, which is a normal pressure in the combustion chamber 5, the pressure term (PTDC / Pstd) ^−0.16 in the equation (49) becomes a small value. Accordingly, it is also possible to define the reference laminar combustion speed SLstd only by the reference temperature Tstd with the pressure term (PTDC / Pstd) ^−0.16 being a constant value. Therefore, when the reference temperature Tstd is 550 [K], the reference laminar combustion speed SLstd is 1.0 [m / sec], and the pressure term is 0.7, the temperature TTDC at the compression top dead center and the laminar combustion The relationship with the speed SL2 can be approximately defined by the following equation.

SL2 = f16 (TTDC)
= 1.0 × 0.7 × (TTDC / 550) ^2.18
... (50)
In step 189, the turbulence intensity ST2 of the gas flow during the main combustion period is calculated. The turbulence intensity ST2 of the gas flow can be expressed by the following equation as a function of the engine speed NRPM, similarly to the turbulence intensity ST1 of the gas flow during the initial combustion period.

ST2 = f17 (NRPM) = C2 × NRPM (20)
Where C2 is a constant,
It is also possible to obtain the turbulence intensity ST2 from a table using the rotation speed as a parameter.

  In step 190, from the laminar combustion speed SL2 [m / sec] and the turbulence intensity ST2 of the gas flow in the main combustion period, the combustion speed FLAME2 [m / sec] in the main combustion period is calculated by the following equation.

FLAME2 = SL2 × ST2 (21)
However, SL2: Laminar burning velocity [m / sec],
Equation (21) considers the contribution to the combustion speed associated with gas turbulence, as in Equation (18).

  In step 191, the main combustion period BURN2 [deg] is calculated by the following equation similar to the equation (19).

BURN2 = {(NRPM × 6) × (BR2 × VTDC)}
/ (RPROBA × AF2 × FLAME2) (22)
However, AF2: reaction area [m ² ] of the flame kernel,
Here, BR2 on the right side of the equation (22) is a change amount of the combustion mass ratio from the start timing to the end timing of the main combustion period. Since the combustion mass ratio BR becomes 2% at the end of the initial combustion period, and then the main combustion period starts and the combustion mass ratio BR reaches 60% and the main combustion period ends, BR2 = 60 % -2% = 58% is set. AF2 is an average reaction area in the growth process of the flame kernel, and is set to a fixed value experimentally determined in advance, like AF1 in the equation (19).

  During the main combustion period, the combustion chamber volume changes with the compression top dead center interposed therebetween. That is, it can be considered that the compression top dead center position exists at the approximate center between the start timing of the main combustion period and the end timing of the main combustion period. In addition, the combustion chamber volume does not change much in the vicinity of the compression top dead center even if the crank angle changes. Therefore, the combustion chamber volume in the main combustion period is represented by the combustion chamber volume VTDC at the compression top dead center.

  FIG. 13 is for calculating the basic ignition timing MBTCAL [degBTDC], and is executed at regular intervals (for example, every 10 msec). This is executed following the flow that is executed later in FIGS.

  In step 41, the initial combustion period BURN1 [deg] calculated in step 171 of FIG. 10, the main combustion period BURN2 [deg] calculated in step 191 of FIG. 12, and calculated in step 17 of FIG. The reference crank angle θPMAX [degATDC] and the ignition timing dead time equivalent crank angle IGNDEAD [deg] are read.

  The calculation of the ignition dead time equivalent crank angle IGNDEAD will be described with reference to the flow of FIG. 29 (refer to Japanese Patent Application No. 2003-109040 for details). The flow of FIG. 29 is executed before the execution of the flow of FIG. 13 after calculating the combustion periods BURN1 and BURN2 shown in the flow of FIGS.

  Here, the ignition dead time DEADTIME [μsec] is first calculated, and the ignition dead time DEADTIME is converted into a crank angle. The dead ignition time DEADTIME is a delay time from when the ignition command is received until when combustion actually starts in the combustion chamber 5. There are three factors that influence the ignition dead time, namely, the pressure and temperature of the combustion chamber 5 and the air-fuel ratio of the air-fuel mixture, and the influence of these factors on the ignition dead time is calculated as independent of each other.

29, in step 201, the basic ignition timing MBTCAL [degBTDC] calculated in step 43 of FIG. 13, the volume VIVC [m ³ ] at the intake valve closing timing of the combustion chamber 5 calculated in step 12 of FIG. The temperature TINI [K] at the intake valve closing timing of the combustion chamber 5 calculated in step 13 of FIG. 5, the pressure PINI [Pa] at the intake valve closing timing of the combustion chamber 5 calculated in step 14 of FIG. The target equivalent ratio TFBYA and the engine speed NRPM [rpm] are read.

  In steps 202 to 205, the effective compression ratio Ec (= V0 / VIVC) at the combustion start timing is calculated in the same manner as steps 162 to 167 in FIG. 10, and the temperature of the combustion chamber 5 at the combustion start timing is calculated using the compression ratio Ec. T0 and the pressure P0 at the combustion start timing of the combustion chamber 5 are calculated.

T0 = TINI × Ec ^ (κ−1) (10-1)
P0 = PINI × Ec ^ κ (10-2)
Where κ: specific heat ratio,
In step 206, the ignition dead time basic value DEADTIME0 [μsec] is calculated by the following equation as a function of the pressure P0 at the combustion start timing of the combustion chamber 5.

DEADTIME0 = f21 (P0) (10-3)
For example, the ignition dead time basic value DEADTIME0 may be obtained by searching a table having the contents shown in FIG. 30 from the pressure P0 at the combustion start timing of the combustion chamber 5. The ignition dead time basic value DEADTIME0 is a value when the theoretical air-fuel ratio and the temperature at the combustion start timing of the combustion chamber 5 are the reference temperature, and increases as the pressure P0 increases as shown in FIG. The reason for the increased ignition dead time at high pressure is that the air-fuel mixture is less likely to vaporize at high pressure.

  Steps 207 to 209 are parts for calculating two correction coefficients for the ignition dead time. That is, when the temperature at the combustion start timing of the combustion chamber 5 falls below the reference temperature or the air-fuel ratio of the air-fuel mixture deviates from the stoichiometric air-fuel ratio, the ignition dead time basic value DEADTIME0 does not match the actual value, so the temperature correction coefficient Kdtmt An excess air ratio correction coefficient Kdtmlmb is introduced.

  First, at step 207, a temperature correction coefficient Kdtmt for the ignition dead time is calculated by the following equation as a function of the temperature T0 at the combustion start timing of the combustion chamber 5.

Kdtmt = f22 (T0) (10-4)
For example, the temperature correction coefficient Kdtmt may be obtained by searching a table containing FIG. 31 from the temperature T0 at the combustion start timing of the combustion chamber 5. The temperature correction coefficient Kdtmt is 1.0 when the temperature is equal to or higher than the reference temperature as shown in FIG. That is, as the temperature T0 is lower, the ignition dead time becomes longer.

  In step 208, the reciprocal of the target equivalent ratio TFBYA is determined as the target excess air ratio TLAMBDA, that is, the target excess air ratio TLAMBDA is obtained by the following equation.

TLAMBDA = 1 / TFBYA (10-5)
In step 209, the excess air ratio correction coefficient Kdtmlmb for the ignition dead time is calculated by the following equation as a function of the target excess air ratio TLAMBDA.

Kdtmlmb = f23 (TLAMBDA) (10-6)
For example, the excess air ratio correction coefficient Kdtmlmb may be obtained by searching a table containing the contents of FIG. 32 from the target excess air ratio TLAMBDA. The excess air ratio correction coefficient Kdtmlmb is a minimum value of 1.0 when the target excess air ratio TLAMBDA is 1 (that is, the stoichiometric air-fuel ratio) as shown in FIG. 32, and is a value that is larger or larger than this. That is, the ignition dead time becomes longer even when the air-fuel ratio deviates from the stoichiometric air-fuel ratio to the rich side or to the lean side.

  In step 210, the ignition dead time basic value DEADTIME0 is corrected by two correction coefficients Kdtmt and Kdtmlmb, that is, the ignition dead time DEADTIME [μsec] is calculated by the following equation.

DEADTIME = DEADTIME0 × Kdtmt × Kdtmlmb
... (10-7)
Here, the pressure P0 is corrected to the basic value (DEADTIME0) by the temperature T0 and the air-fuel ratio (TLAMBDA). However, the present invention is not limited to this. The pressure P0 and the empty value are determined based on the temperature T0. It may be corrected by the fuel ratio (TLAMBDA) or may be corrected by the temperature T0 and the pressure P0 based on the air / fuel ratio (TLAMBDA). In addition, the factors affecting the ignition dead time are not limited to this, and there are charging efficiency and inert gas rate (internal inert gas rate MRES and external inert gas rate). Thus, the ignition dead time can be calculated. However, it is not always necessary to consider all five factors, and the ignition dead time may be calculated using at least one factor as a parameter.

  Finally, in step 211, the crank angle IGNDEAD [deg] corresponding to the ignition dead time is calculated. The ignition dead time equivalent crank angle IGNDEAD is a crank angle section from the timing at which the engine controller 31 outputs a signal for cutting off the primary current of the ignition coil 13 until the ignition plug 14 actually ignites, and can be expressed by the following equation. .

IGNDEAD = f5 (DEADTIME, NRPM) (10-8)
The expression (10-8) is for calculating the ignition dead time equivalent crank angle IGNDEAD, which is the crank angle corresponding to the ignition dead time DEADTIME, from the engine speed NRPM.

  Returning to FIG. 13, in step 42, the sum of the initial combustion period BURN1 and the main combustion period BURN2 is calculated as the combustion period BURN [deg].

  In step 43, the basic ignition timing MBTCAL [degBTDC] is calculated by the following equation.

MBTCAL = BURN−θPMAX + IGNDEAD (23)
In step 44, the value obtained by subtracting the ignition dead time equivalent crank angle IGNDEAD from the basic ignition timing MBTCAL is calculated as the previous combustion start timing MBTCYCL [degBTDC].

  Assuming that the basic ignition timing MBTCAL calculated in step 43 is used as the ignition timing command value for the current cycle, the previous combustion start timing MBTCYCL calculated in step 44 is calculated until the ignition timing for the next cycle is reached. Used in step 162.

  FIG. 14 is for calculating the internal inert gas ratio MRESFR in the combustion chamber 5 and is executed at regular intervals (for example, every 10 msec). This flow is executed prior to the flow of FIG.

  In step 51, the output of the air flow meter 32 and the target equivalent ratio TFBYA are read. In step 52, based on the output of the air flow meter 32, a new air amount (cylinder fresh air amount) MACYL flowing into the combustion chamber 5 is calculated. As a method for calculating the cylinder fresh air amount MACYL, a known method may be used (see JP 2001-50091 A).

  In step 53, an internal inert gas amount MRES in the combustion chamber 5 is calculated. The calculation of the internal inert gas amount MRES will be described with reference to the flow of FIG.

  In FIG. 15 (subroutine of step 53 in FIG. 14), in step 61, an inert gas amount MRESCYL at the exhaust valve closing timing EVC in the combustion chamber 5 is calculated. The calculation of the inert gas amount MRESCYL will be further described with reference to the flowchart of FIG.

  In FIG. 16 (subroutine of step 61 in FIG. 15), in step 71, the exhaust valve closing timing EVC [degBTDC], the exhaust temperature TEXH [K] detected by the temperature sensor 45, the exhaust pressure PEXH [kPa] detected by the pressure sensor 46 ].

  Here, just as the intake valve closing timing IVC is known from the command value given to the intake VTC mechanism 27, the exhaust valve closing timing EVC is also known from the command value given to the exhaust VTC mechanism 28.

  In step 72, the volume VEVC of the combustion chamber 5 at the exhaust valve closing timing EVC is calculated. This may be obtained by searching a table using the exhaust valve closing timing as a parameter, similarly to the volume VIVC at the intake valve closing timing IVC. That is, when the exhaust valve VTC mechanism 28 is provided, the volume VEVC at the exhaust valve closing timing EVC of the combustion chamber 5 may be obtained by searching the table shown in FIG. 23 from the exhaust valve closing timing EVC. When the exhaust VTC mechanism 28 is not provided, a constant value can be given.

  Although not shown, when a mechanism for changing the compression ratio is provided, the combustion chamber volume VEVC at the exhaust valve closing timing corresponding to the amount of change in the compression ratio is obtained from the table. When a mechanism for changing the compression ratio in addition to the exhaust VTC mechanism 28 is provided, the combustion chamber volume at the exhaust valve closing timing is obtained by searching a map corresponding to the exhaust valve closing timing and the compression ratio change amount.

  In step 73, the gas constant REX of the inert gas in the combustion chamber 5 is obtained by searching the table shown in FIG. 24 from the target equivalent ratio TFBYA. As shown in FIG. 24, the gas constant REX of the inert gas is the smallest when the target equivalent ratio TFBYA is 1.0, that is, the stoichiometric air-fuel ratio, and becomes larger whether it is larger or smaller.

  In step 74, the temperature TEVC of the combustion chamber 5 at the exhaust valve closing timing EVC is estimated based on the exhaust temperature TEXH. For simplicity, the exhaust temperature TEXH may be set as TEVC as it is. Note that the temperature TEVC of the combustion chamber 5 at the closing timing of the exhaust valve changes depending on the amount of heat corresponding to the fuel injection amount of the injector 21. Therefore, if such characteristics are taken into consideration, the calculation accuracy of TEVC is improved.

  In step 75, the pressure PEVC at the exhaust valve closing timing EVC of the combustion chamber 5 is calculated based on the exhaust pressure PEXH. Simply, the exhaust pressure PEXH may be set to PEVC.

  In step 76, from the volume VEVC at the exhaust valve closing timing EVC of the combustion chamber 5, the temperature TEVC at the exhaust valve closing timing EVC, the pressure PEVC at the exhaust valve closing timing EVC, and the inert gas gas constant REX, the exhaust valve of the combustion chamber 5 is obtained. The inert gas amount MRESCYL at the closing timing EVC is calculated by the following equation.

MRESCYL = (PEVC × VEVC) / (REX × TEVC) (24)
When the calculation of the inert gas amount MRESCYL at the exhaust valve closing timing EVC of the combustion chamber 5 is completed in this way, the processing returns to FIG. 15, and in step 62, the overlap of the intake and exhaust valves 15 and 16 (“O / L” in the figure). The amount of inactive gas MRESOL during the overlap is calculated, which is the amount of inert gas that is blown back from the exhaust side to the intake side.

  The calculation of the inert gas amount MRESOL will be described with reference to the flowchart of FIG.

  In FIG. 17 (subroutine of step 62 in FIG. 15), in step 81, the intake valve opening timing IVO [degBTDC], the exhaust valve closing timing EVC [degBTDC], and the exhaust valve of the combustion chamber 5 calculated in step 74 of FIG. The temperature TEVC at the closing timing EVC is read.

  Here, since the intake valve opening timing IVO is a timing earlier than the intake valve closing timing IVC by the opening angle of the intake valve 15, the opening angle of the intake valve 15 (which is known in advance) from the intake valve closing timing IVC. Can be sought.

  In step 82, the intake valve exhaust timing VTCOL [deg] is calculated from the intake valve opening timing IVO and the exhaust valve closing timing EVC by the following equation.

VTCOL = IVO + EVC (25)
For example, the intake valve opening timing IVO is at the intake top dead center position when the intake VTC mechanism 27 actuator is not energized, and the intake valve opening timing is advanced from the intake top dead center when the intake VTC mechanism 27 actuator is energized. The exhaust valve closing timing EVC is at the exhaust top dead center when the exhaust VTC mechanism 28 actuator is not energized, and the exhaust valve closing timing EVC is at the exhaust top dead center when the exhaust valve VTC mechanism 28 actuator is energized. In the case of more advanced characteristics, the sum of IVO and EVC becomes the overlap amount VTCOL of the intake and exhaust valves.

  In step 83, the accumulated effective area ASUMOL during the overlap is calculated by searching the table shown in FIG. 25 from the overlap amount VTCOL of the intake and exhaust valves. As shown in FIG. 25, the integrated effective area ASUMOL during the overlap is a value that increases as the overlap amount VTCOL of the intake and exhaust valves increases.

  Here, FIG. 26 is an explanatory diagram of the integrated effective area ASUMOL during the overlap of the intake and exhaust valves, where the horizontal axis indicates the crank angle, and the vertical axis indicates the respective opening areas of the intake valve 12 and the exhaust valve 15. Yes. The effective opening area at any time during the overlap is the smaller of the exhaust valve opening area and the intake valve opening area. The integrated effective area ASUMOL in the entire period during the overlap is an integral value (hatched portion in the figure) during the period in which the intake valve 15 and the exhaust valve 16 are open.

  By calculating the accumulated effective area ASUMOL during the overlap in this way, the overlap amount between the intake valve 15 and the exhaust valve 16 can be approximated as one orifice (outflow hole), and the state of the exhaust system and The gas flow rate passing through the virtual orifice can be simply calculated from the state of the intake system.

  In step 84, the specific heat ratio SHEATR of the inert gas remaining in the combustion chamber 5 is obtained by searching the map shown in FIG. 27 from the target equivalent ratio TFBYA and the temperature TEVC at the exhaust valve closing timing EVC of the combustion chamber 5. calculate. As shown in FIG. 27, the specific heat ratio SHEATR of the inert gas remaining in the combustion chamber is the smallest when the target equivalent ratio TFBYA is in the vicinity of 1.0, and it is larger or smaller than that. Moreover, under the condition where the target equivalent ratio TFBYA is constant, the temperature TEVC at the exhaust valve closing timing EVC of the combustion chamber 5 becomes smaller as the temperature becomes higher.

  In step 85, a supercharging determination flag TBCRG and a choke determination flag CHOKE are set. The setting of the supercharging determination flag TBCRG and the choke determination flag CHOKE will be described with reference to the flowchart of FIG.

  In FIG. 18 (subroutine of step 85 in FIG. 17), in step 101, the intake pressure PIN detected by the intake pressure sensor 44 and the pressure PEVC at the exhaust valve closing timing EVC of the combustion chamber 5 calculated in step 75 of FIG. Is read.

  In step 102, an intake exhaust pressure ratio PINBYEX is calculated from the intake pressure PIN and the pressure PEVC at the exhaust valve closing timing EVC of the combustion chamber 5 by the following equation.

PINBYEX = PIN / PEVC (26)
This intake / exhaust pressure ratio PINBYEX is an unknown number, and 1 is compared with this in step 103. When the intake / exhaust pressure ratio PINBYEX is 1 or less, it is determined that there is no supercharging, and the routine proceeds to step 104 where the supercharging determination flag TBCRG (initially set to zero) = 0.

  If the intake / exhaust pressure ratio PINBYEX is greater than 1, it is determined that there is supercharging, and the routine proceeds to step 105 where the supercharging determination flag TBCRG = 1.

  In Step 106, the specific heat ratio MIXAIRSHR of the mixture is obtained by searching the table shown in FIG. 28 from the target equivalent ratio TFBYA read in Step 51 of FIG. Replace with SHEATR. As shown in FIG. 28, the specific heat ratio MIXAIRSHR of the air-fuel mixture is a value that increases as the target equivalent ratio TFBYA decreases.

  In steps 106 and 107, the specific heat ratio SHEATR of the inert gas is replaced with the specific heat ratio MIXAIRSHR of the air-fuel mixture in consideration of supercharging such as turbocharging or inertial supercharging. That is, during supercharging, the gas flow during the overlap of the intake / exhaust valve is directed (blows through) from the intake system to the exhaust system. In this case, the specific heat ratio of the gas passing through the virtual orifice is the specific heat of the inert gas. By changing the ratio to the specific heat ratio of the air-fuel mixture, the amount of gas blown through is accurately estimated, and the amount of internal inert gas is accurately calculated.

  In step 108, based on the specific heat ratio SHEATR of the inert gas calculated in step 84 of FIG. 17 or steps 106 and 107 of FIG. 18, the minimum and maximum choke determination thresholds SLCHOKE and SLCHOKEH are calculated by the following equations. To do.

SLCHOKER = {2 / (SHEATR + 1)}
^ {SHEATR / (SHEATR-1)} (27a)
SLCHOKEH = {− 2 / (SHEATR + 1)}
^ {-SHEATR / (SHEATR-1)} (27b)
These choke determination threshold values SLCHOKE and SLCHOKEH calculate the limit value for choking.

  If it is difficult to calculate each power of (27a) right side and (27b) right side in step 108, the calculation results of equations (27a) and (27b) are used as the minimum choke determination threshold value SLCHOKEL table and maximum choke determination. The threshold value SLCHOKEH table may be stored in advance in the memory of the engine controller 31 and obtained by searching the table from the specific heat ratio SHEATR of the inert gas.

  In step 109, it is determined whether or not the intake / exhaust pressure ratio PINBYEX is within the range of not less than the minimum choke determination threshold value SLCHOKEEL and not more than the maximum choke determination threshold value SLCHOKEH, that is, not in the choke state. If the intake / exhaust pressure ratio PINBYEX is within the range, it is determined that there is no choke, and the routine proceeds to step 110 where the choke determination flag CHOKE (initially set to zero) = 0.

  If the intake / exhaust pressure ratio P1NBYEX is not within the range, it is determined that choke is present, and the routine proceeds to step 111 where the choke determination flag CHOKE = 1.

  When the setting of the supercharging determination flag and the choke determination flag is thus completed, the process returns to FIG. 17 and the following four cases are performed in steps 86 to 88.

<1> When supercharging determination flag TBCRG = 0 and choke determination flag CHOKE = 0 <2> When supercharging determination flag TBCRG = 0 and choke determination flag CHOKE = 0 <3> Supercharging determination flag TBCRG = 0 and choke When determination flag CHOKE = 1 <4> When supercharging determination flag TBCRG = 1 and choke determination flag CHOKE = 0 When the above <1>, the routine proceeds to step 89 and overlaps when there is no supercharging and no choke If the average blown back inert gas flow rate MRESOLtmp1 is <2>, the process proceeds to step 90, and the blown back inert gas flow rate MRESOLtmp2 during the overlap when there is no supercharging and choke is present is step 91 when <3>. Go on to overlap when there is supercharging and no choke The average blowback inert gas flow rate MRESOLtmp3 in the case of <4> is advanced to step 92 to calculate the blowback inert gas flow rate MRESOLtmp4 with supercharging and with choke, respectively. Move to active gas flow rate MRESOLtmp.

Here, calculation of the inert gas flow rate MRESOLtmp1 during the overlap when there is no supercharging and no choke will be described with reference to the flow of FIG. 19. In step 121 of FIG. 19 (subroutine of step 89 in FIG. 17), step 73 in FIG. , 75, the gas constant REX of the inert gas, and the pressure PEVC when the exhaust valve of the combustion chamber 5 is closed are read.

  In step 122, based on the gas constant REX of the inert gas and the temperature TEVC at the exhaust valve closing timing of the combustion chamber 5 read in step 81 of FIG. 17, a density term MRSOLD used in a gas flow rate calculation formula described later. Is calculated by the following equation.

MRSOLD = SQRT {1 / (REX × TEVC)} (28)
Here, “SQRT” on the right side of equation (28) is a function that calculates the square root of the value in the parenthesis on the right.

  If the square root of the density term MRSOLD is difficult to calculate, the calculation result of equation (28) is stored in advance in the memory of the engine controller 31 as a map, and the gas constant REX and the temperature of the combustion chamber 5 at the exhaust valve closing timing are stored. You may obtain | require by searching the map from TEVC.

  In step 123, based on the specific heat ratio SHEATR of the inert gas calculated in step 84 of FIG. 17 and the intake / exhaust pressure ratio PINBYEX calculated in step 102 of FIG. The pressure difference term MRSOLP used is calculated by the following equation.

MRSOLP = SQRT [SHEATR / (SHEATR-1)
× {PTNBYEX ^ (2 / SHEATR)
-PTNBYEX ^ ((SHEATR + 1) / SHEATR)}] (29)
In step 124, from the density term MRSOLD, the pressure difference term MRSOLP, and the pressure PEVC at the exhaust valve closing timing of the combustion chamber 5, the blowback inert gas flow rate MRESOLtmp1 during the overlap without supercharging and without choke is expressed by In step 125, the calculated value is transferred to the blown back inert gas flow rate MRESOLtmp.

MRESOLtmp1 = 1.4 × PEVC × MRSOLD × MRSOLP (30)
Next, the calculation of the blown back inert gas flow rate when there is no supercharging and when there is a choke will be described with reference to the flow of FIG. 20. In FIG. 20 (subroutine of step 90 of FIG. 17), in steps 131 and 132, steps 121 and 122 of FIG. In the same manner as described above, the gas constant REX of the inert gas and the pressure PEVC at the time of closing the exhaust valve of the combustion chamber 5 are read, and the density term MRSOLD is calculated from the above equation (28).

  In step 133, the choke pressure difference term MRSOLPC is calculated by the following equation based on the specific heat ratio SHEATR of the inert gas calculated in step 84 of FIG.

MRSOLPC = SQRT [SHEATR × {2 / (SHEATR + 1)} ^ {(SHEATR + 1) / [SHEATR-1)}] (31)
If it is difficult to calculate the power and square root of equation (31), the calculation result of equation (31) is stored in advance in the memory of the engine controller 31 as a table of choke pressure difference term MRSOLPC. Alternatively, it may be obtained by searching the table from the specific heat ratio SHEATR of the inert gas.

  In step 134, from the density term MRSOLD, the choke pressure difference term MRSOLPC, and the pressure PEVC at the time of closing the exhaust valve of the combustion chamber 5, the recirculating inert gas flow rate MRESOLtmp2 during the overlap without supercharging and with choke is In step 135, the calculated value is transferred to the blown back inert gas flow rate MRESOLtmp.

MRESOLtmp2 = PEVC × MRSOLD × MRSOLPC (32)
Next, calculation of the blowback gas flow rate with supercharging and without choke will be described with reference to the flow of FIG. 21. In FIG. 21 (subroutine of step 91 in FIG. 17), in step 141, the intake pressure PIN detected by the intake pressure sensor 44 is explained. Is read.

  In step 142, the supercharging pressure difference term MRSOLPT is calculated from the specific heat ratio SHEATR of the inert gas calculated in steps 106 and 107 in FIG. 18 and the intake / exhaust pressure ratio PINBYEX calculated in step 102 in FIG. Is calculated by the following equation.

MRSOLPT = SQRT [SHEATR / (SHEATR-1)
× {PINBYEX ^ (-2 / SHEATR)
-PINBYEX ^ (-(SHEATR + 1) / SHEATR)}] (33)
If the power calculation and the square root calculation of Expression (33) are difficult, the calculation result of Expression (33) is stored in advance in the memory of the engine controller 31 as a map of the supercharging pressure difference term MRSOLPT, You may obtain | require by searching the map from the specific heat ratio SHEATR of an inert gas, and the intake-exhaust pressure ratio PINBYEX.

  In step 143, based on the pressure difference term at the time of supercharging MRSOLPT and the intake pressure PIN, a recirculating inert gas flow rate MRESOLtmp3 during overlap with supercharging and without choke is calculated by the following equation, and the calculated value is calculated in step 143. At 144, the flow returns to the blown back inert gas flow rate MRESOLtmp.

MRESOLtmp3 = −0.152 × PIN × MRSOLPT (34)
Here, by setting the blown back inert gas flow rate MRESOLtmp3 in the equation (34) to a negative value, the gas flow rate of the air-fuel mixture blown from the intake system to the exhaust system during the overlap can be expressed.

Next, the calculation of the flow rate of the inert gas blown back during overlap with supercharging and with choke will be described with reference to the flow of FIG. 22. In FIG. 22 (subroutine of step 92 in FIG. 17), in steps 151 and 152, the step of FIG. Similarly to 141, the intake pressure PIN detected by the intake pressure sensor 44 is read, and the choke pressure difference term MRSOLPC is calculated by the aforementioned equation (31) as in step 132 of FIG.

  In step 153, based on the choke pressure difference term MRSOLPC and the intake pressure PIN, the overlap blowback gas flow rate MRESOLtmp4 with supercharging and with choke is calculated by the following equation, and the calculated value is exceeded in step 154. Transfer to blown inert gas flow rate MRESOLtmp during lap.

MRESOLtmp4 = −0.108 × PIN × MRSOLPC (35)
Here, the blown back inert gas flow rate MRESOLtmp4 of the equation (35) can also represent a gas flow rate of the air-fuel mixture blown from the intake side to the exhaust side during the overlap, similarly to MRESOLtmp3.

  In this way, when the calculation of the blown back inert gas flow rate MRESOLtmp during the overlap divided according to the combination of the presence or absence of supercharging and the presence or absence of choke is completed, the process returns to FIG. From the inert gas flow rate MRESOLtmp and the integrated effective area ASUMOL during the overlap period, the blown back inert gas amount MRESOL during the overlap is calculated by the following equation.

MRESOL = (MRESOLtmP × ASUMOL × 60)
/ (NRPM × 360) (36)
When the calculation of the blown back inert gas amount MRESSOL is completed in this way, the process returns to FIG. 15, and in step 63, the inert gas amount MRESCYL at the exhaust valve closing timing EVC in the combustion chamber 5 and the blowback during overlap are returned. The internal inert gas amount MRES is calculated by adding the gas amount MRESOL, that is, the following equation.

MRES = MRESCYL + MRESSOL (37)
As described above, the flow rate of the inert gas blown back during overlap (MRESOLtmp3, MRESOLtmp4) becomes negative when there is supercharging, and therefore, the amount of blown inert gas MRESOL during overlap of the above equation (36) also becomes negative. According to the equation (37), the internal inert gas amount is reduced by the amount of the blown back inert gas amount MRESOL during the overlap.

  When the calculation of the internal inert gas amount MRES is completed in this manner, the flow returns to FIG. 14, and in step 54, the internal inert gas rate MRESFR is calculated by the following equation using the internal inert gas amount MRES and the target equivalent ratio TFBYA. (The ratio of the amount of internal inert gas to the total amount of gas in the combustion chamber 5) is calculated.

MRESFR = MRES
/{MRES+MACYL×(1+TFBYA/14.7)} (38)
This completes the calculation of the internal inert gas ratio MRESFR.

  As described above, according to the present embodiment, the internal inert gas amount MRES is configured by the inert gas amount MRESCYL when the exhaust valve of the combustion chamber 5 is closed and the blowback gas amount MRESOL during the overlap of the intake and exhaust valves. In this case, the temperature TEV and pressure PEVC at the closing timing of the exhaust valve in the combustion chamber 5 are calculated (steps 74 and 75 in FIG. 16), and these temperature TEVC, pressure PEVC and inert gas are calculated. Since the inert gas amount MRESCYL at the exhaust valve closing timing of the combustion chamber 5 is calculated from the state constant (Equation (24) above) based on the gas constant REX (see step 76 of FIG. 16), Even during transient operation in which the state quantity (PEVC, VEVC, TEVC) inside the combustion chamber 5 changes every moment, the combustion chamber 5 can be accurately obtained regardless of the operating conditions. The inert gas amount MRESCYL in the exhaust valve closing timing can be calculated (estimated).

  In addition, based on the temperature TEVC and pressure PEVC of the combustion chamber 5 when the exhaust valve is closed, the gas constant REX and specific heat ratio SHEATR of the inert gas, and the intake pressure PIN, the blown back inert gas flow rate (MRESOLtmp1, MRESOLtmp2) during the overlap is determined. This was calculated (see FIGS. 19 and 20), and this gas flow rate was multiplied by the integrated effective area ASUMOL during the overlap to calculate the amount of blown back gas MRESOL during the overlap (see step 93 in FIG. 17). Therefore, the overlapped blow-back gas amount MRESOL can be calculated (estimated) with high accuracy.

  As described above, if both the inert gas amount MRESCYL and the overlapped blow-back gas amount MRESSOL in the combustion chamber 5 at the closing timing of the exhaust valve can be calculated (estimated) with high accuracy, the internal inert gas amount MRES, which is the sum of them, can also be accurately calculated. The internal inert gas ratio MRESFR calculated based on the internal inert gas amount MRES that can be estimated with high accuracy can be calculated (estimated). By making use of the temperature TINI at the intake valve closing timing IVC (see step 13 in FIG. 5), the temperature TINI at the intake valve closing timing IVC in the combustion chamber 5 can be accurately calculated. In addition, the engine can be appropriately controlled by making use of the internal inert gas amount MRES that can be accurately estimated for the fuel injection amount, valve opening / closing timing (overlap amount), and the like. .

  Further, since the gas constant REX of the inert gas and the specific heat ratio SHEATR of the inert gas are values corresponding to the target equivalent ratio TFBYA (see FIGS. 24 and 27), the operation is performed at an air-fuel ratio that deviates from the theoretical air-fuel ratio. (E.g. during lean operation where the air / fuel ratio is leaner than the stoichiometric air / fuel ratio, in order to stabilize the engine's originally unstable state, such as during cold start, the engine is on the rich side of the stoichiometric air / fuel ratio. Immediately after starting the engine operating at the air-fuel ratio, since the same large output is required, the exhaust valve closing timing of the combustion chamber 5 is also applied during full load operation where the air-fuel ratio is richer than the stoichiometric air-fuel ratio. Inert gas amount MRESCYL, overlapped blow-back gas amount MRESSOL, total internal inert gas amount MRES, and internal inert gas ratio MRESFR based on this are accurately calculated. It can be.

  Further, since the accumulated effective area ASUMOL during the overlap period is defined as the area of the virtual orifice, and it is assumed that the exhaust gas is blown through the virtual orifice from the combustion chamber 5 to the intake system, the amount of blown back inert gas amount MRESOL during the overlap is assumed. Calculation is simplified.

  This concludes the description of the calculation of the internal inert gas ratio MRESFR in the combustion chamber 5, and next, the calculation of the estimated shaft torque generated by the engine will be described.

  In Step 2 of FIG. 2, shaft torque estimated values Trquest1 and Trquest2 are calculated. Here, there are two types of estimated shaft torque values: an estimated value of shaft torque when ignited at the basic ignition timing MBTCAL, and an estimated value of shaft torque when ignited at an ignition timing retarded (retarded) from the basic ignition timing MBTCAL. Therefore, the estimated shaft torque when ignited at the basic ignition timing MBTCAL is the first shaft torque estimated value Trquest1, and the estimated shaft torque when ignited at the ignition timing retarded from the basic ignition timing MBTCAL is the second axis. It distinguishes as torque estimated value Trquest2.

  Further, the initial combustion period BURN1, the ignition dead time equivalent crank angle IGNDEAD, and the main combustion period BURN2 are calculated in the process of obtaining each shaft torque estimated value, and these are also calculated in the process of obtaining the first shaft torque estimated value Trquest1, IGNDEAD and BURN2 are given “ref1” to be “BURN1ref1”, “IGNEADref1”, and “BURN2ref1”, which are respectively referred to as a first initial combustion period, a first ignition dead time equivalent crank angle, and a first main combustion period. Further, BURN1, IGNDEAD, and BURN2 calculated in the process of obtaining the second shaft torque estimated value Trquest2 are added with “ref2” to be “BURN1ref2,” “IGNEADref2,” and “BURN2ref2,” respectively. The crank angle corresponding to the second ignition dead time is referred to as a second main combustion period.

  Further, an approximate expression of the combustion mass ratio calculated in the process of obtaining each shaft torque estimated value and an estimated positive work estimated value are calculated, respectively, and these are also the combustion mass ratios calculated in the process of obtaining the first shaft torque estimated value Trquest1. “1” is added to the approximate expression and the indicated positive work estimated value to give “X1 ′” and “Wpos1”, respectively, which are referred to as the first approximate expression of the combustion mass ratio and the first illustrated positive work estimated value, respectively. Further, an approximate expression of the combustion mass ratio calculated in the process of obtaining the second shaft torque estimated value Trquest2, and “2” are added to the illustrated positive work estimated values to give “X2 ′” and “Wpos2”, respectively, which are burned respectively. This is referred to as a second approximate expression of the mass ratio and a second illustrated positive work estimated value.

  The calculation of the first shaft torque estimated value Trquest1 will be described with reference to the flow of FIG. 33A, and the calculation of the second shaft torque Trquest2 will be described with reference to the flow of FIG. 33B. In this case, since the calculation method of the shaft torque estimated value itself does not change because the ignition timing is different, the calculation of the first shaft torque estimated value Trquest1 will be described in detail, and then the second shaft torque estimated value Trquest2 will be calculated. Mention.

  In FIG. 33 (A) (subroutine of step 2 of FIG. 2), in step 221, the basic ignition timing MBTCAL [degBTDC] calculated in step 43 of FIG. 13 is read.

  In step 222, the value of the ignition timing θign [degBTDC] is transferred to the previous ignition timing value θignz [degBTDC], and the basic ignition timing MBTCAL is set to the ignition timing θign [degBTDC].

  In step 223, the first initial combustion period BURN1ref1, the first main combustion period BURN2ref1, and the first ignition dead time equivalent crank angle IGNDEADref1 are set to the first initial combustion period previous value BURN1ref1z, the first main combustion period previous value BURN2ref1z, the first Shift to the ignition dead time equivalent crank angle previous value IGNDEADref1z.

  Here, the values of BURN1ref1, BURN2ref1, and IGNDEADref1 are values calculated in steps 224 to 226 described later. Therefore, since BURN1ref1, BURN2ref1, and IGNDEADref1 are not calculated at the first time, the value of BURN1 calculated in step 171 of FIG. 10 is set in BURN1ref1z, and the value of BURN2 calculated in step 191 of FIG. 12 is set in BURN2ref1z. The value is set to IGNDEADref1z with the value of IGNDEAD calculated in step 211 in FIG.

  In steps 224 to 226, a first initial combustion period BURN1ref1, a first ignition dead time equivalent crank angle IGNDEADref1, and a first main combustion period BURN2ref1 are calculated. These calculations will be described with reference to the flowcharts of FIGS. 34, 35, and 36.

  First, from FIG. 34, the basic concept of the calculation method of the first initial combustion period BURN1ref1 is the same as the calculation method of the initial combustion period shown in FIG.

In FIG. 34 (subroutine of step 224 of FIG. 33A), in step 241, the cylinder fresh air amount MACYL [kg] and the internal inert gas amount MRES [kg] calculated in steps 52 and 53 of FIG. Amount QINJ [kg], collector temperature TCOL [K] detected by the temperature sensor 43, exhaust temperature TEXH [K] detected by the temperature sensor 45, collector pressure PCOL [Pa] detected by the pressure sensor 44, The volume VIVC [m ³ ] at the intake valve closing timing of the combustion chamber 5 calculated in step 12 of FIG. 5, the engine speed NRPM [rpm], and the reaction probability RPROBA [%] calculated in step 15 of FIG. , Ignition timing previous value θignz [degBTDC] obtained in steps 222 and 223 in FIG. 1 initial combustion period previous value BURN1ref1z [deg], the first main combustion period previous value BURN2ref1z [deg], the first ignition dead time corresponds crank angle immediately preceding value IGNDEADRf1z read [deg].

In step 242, the volume Vmb1p [m ³ ] when the combustion mass ratio of the combustion chamber 5 is 1% is calculated by the following equation using the previous ignition timing value θignz.

Vmb1p = f25 (θignz−IGNEADrefref1z
-BURN1ref1z / 2)
... (61)
In step 243, the effective compression ratio Emb1p when the combustion mass ratio is 1% is calculated by the following equation.

Emb1p = Vmb1p / VIVC (62)
In step 244, the cylinder fresh air amount MACYL [kg] is moved to WIDRY [kg], and the internal inert gas amount MRES [kg] is moved to MASSZ [kg]. These WIDRY and MASSZ are introduced for use only in calculating the shaft torque estimation value, WIDRY is the cylinder fresh air amount, and MASSZ is the internal inert gas amount.

  In step 245, the compression start temperature TC0 [K] of the combustion chamber 5 is calculated by the following equation.

TC0 = {(WIDRY + QINJ) × TCOL + MASSZ × TEXH}
/ (WIDRY + QINJ + MASSZ) (63)
Here, the specific heat of the inert gas and fresh air is made equal to simplify the equation.

  In step 246, the compression start pressure PC0 [Pa] of the combustion chamber 5 is calculated. For this purpose, the collector internal pressure PCOL at the intake valve closing timing IVC detected by the pressure sensor 44 may be set to PC0.

  In steps 247 and 248, the average when the fuel in the combustion chamber 5 is burned using Emb1p, which is the instantaneous compression ratio when the combustion mass ratio is 1%, and BR1p, which is the combustion mass ratio when the combustion mass ratio is 1%. The temperature TC [K] and the average pressure PC [Pa] are calculated by the following equations.

TC = TC0 × Emb1p ^ 0.35
+ CF # × QINJ × BR1p / (MASSZ + WIDRY + QINJ)
... (64)
PC = PC0 * Emb1p ^ 1.35 * TC / TC0 / Emb1p ^ 0.35
... (65)
Where CF #: lower heating value of fuel,
The equations (64) and (65) are equations when it is assumed that the gas is adiabatically compressed in the combustion chamber 5 and combusted with a constant volume change. That is, the first term on the right side of equation (64) is the temperature after adiabatic compression, PC0 × Emb1p ^ 1.35 on the right side of equation (65) is the pressure after adiabatic compression, and the second item on the right side of equation (64). TC / TC0 / Emb1p ^ 0.35 on the right side of equation (65) represents the rate of increase in pressure due to combustion at constant volume change, where the term is the temperature rise due to combustion due to constant volume change.

  The operation for calculating the first initial combustion period in steps 249 to 252 is the same as the operation for calculating the initial combustion period in steps 168 to 171 of FIG. The difference is that the laminar combustion speed SL1 is calculated using the temperature T0 and the pressure P0 at the combustion start timing of the combustion chamber 5 in step 168 of FIG. 10, whereas in step 249 of FIG. , The point where the laminar combustion speed SL1 is calculated using the temperature TC and the pressure PC obtained by the equation (65), and BR1 (= 2 which is the combustion mass ratio at the end of the initial combustion period BURN1 in step 171 of FIG. %) And V0 which is the volume at the combustion start timing of the combustion chamber 5 is calculated, whereas the initial combustion period BURN1 is calculated in step 252 of FIG. The first initial combustion period BURN1ref1 is calculated using 1% (= BR1p) and Vmb1p which is the volume of the combustion chamber 5 when the combustion mass ratio is 1%. In that it is. That is, in step 249, the laminar combustion speed SL1 [m / sec] in the first initial combustion period is calculated, and in step 252 the first initial combustion period BURN1ref1 [deg] is calculated by the following equations.

SL1 = SLstd × (TC / Tstd) ^2.18 × (PC / Pstd) ^−0.16 (66)
Where Tstd: reference temperature [K],
Pstd: Reference pressure [Pa]
SLstd: Reference laminar combustion at reference temperature Tstd and reference pressure Pstd
Speed [m / sec],
BURN1ref1 = {(NRPM × 6) × BR1p × Vmb1p}
/ (RPROBA × AF1 × FLAME1) (67)
However, AF1: Reaction area (fixed value) of flame kernel [m ² ],
The flow of FIG. 35 (subroutine of step 225 of FIG. 33A) is for calculating the first ignition dead time equivalent crank angle IGNDEAD, and the calculation method is the aforementioned ignition dead time equivalent crank angle shown in FIG. It is basically the same as the IGNDEAD calculation method.

29 differs from the calculation method of the ignition dead time equivalent crank angle IGNDEAD shown in FIG. 29 in that the previous combustion start timing MBTCYCL is used to calculate the volume V0 at the combustion start timing of the combustion chamber 5 in FIG. In FIG. 35, only the ignition timing θign is used (see steps 261 and 262 in FIG. 35). Therefore, the differences will be mainly described. In step 261 in FIG. 35, the ignition timing θign [degBTDC] obtained in step 222 in FIG. 33A and the intake air of the combustion chamber 5 calculated in step 12 in FIG. The volume VIVC [m ³ ] at the valve closing timing, the temperature TINI [K] at the intake valve closing timing calculated at step 13 in FIG. 5, and the combustion chamber 5 calculated at step 14 in FIG. The pressure PINI [Pa] at the intake valve closing timing, the target equivalence ratio TFBYA, and the engine speed NRPM [rpm] are read. Among these, the ignition timing θign is used, and in step 262, the volume V0 [m ³ ] at the combustion start timing in the combustion chamber 5 is read. Is calculated.

  In steps 263 to 265, the temperature T0 [K] at the combustion start timing of the combustion chamber 5 and the pressure P0 [Pa] at the combustion start timing of the combustion chamber 5 are calculated in the same manner as in steps 203 to 205 of FIG.

  The flow of FIG. 36 (subroutine of step 226 in FIG. 33A) is for calculating the first main combustion period, and the basic concept of the calculation method of the first main combustion period here is shown in FIG. This is the same as the calculation method of the main initial combustion period shown.

In step 281, the cylinder fresh air amount MACYL [kg], the internal inert gas amount MRES [kg], the fuel amount QINJ [kg] calculated in steps 52 and 53 of FIG. The temperature TCOL [K], the exhaust temperature TEXH [K] detected by the temperature sensor 45, the collector internal pressure PCOL [Pa] detected by the pressure sensor 44, and the intake air of the combustion chamber 5 calculated in step 12 of FIG. Volume VIVC [m ³ ] at the valve closing timing, engine rotation speed NRPM [rpm], reaction probability RPROBA [%] calculated in step 15 of FIG. 5, ignition obtained in step 222 of FIG. Timing θign, ignition timing previous value θignz, first initial combustion period previous value BUR obtained in step 223 of FIG. N1ref1z [deg], first main combustion period previous value BURN2ref1z [deg], first ignition dead time equivalent crank angle previous value IGNDEADref1z [deg], first initial combustion period BURN1ref1 [deg] calculated in step 252 of FIG. ], The first ignition dead time equivalent crank angle IGNDEADref1 [deg] calculated in step 271 of FIG. 35 is read.

  Steps 282 to 285 are parts for calculating the volume Vmb31p when the combustion mass ratio of the combustion chamber 5 is 31%. First, in step 282, the previous crank angle value θmb31pz [degBTDC] when the combustion mass ratio is 31% is calculated by the following equation.

θmb31pz = θignz-IGNEADref1z-BURN1ref1z
-BURNN2ref1z / 2
... (68)
In step 283, a difference θ2pdif [deg] between the previous value of the crank angle when the combustion mass ratio is 2% and the current value of the crank angle when the combustion mass ratio is 2% is calculated by the following equation.

θ2pdif = (θignz−IGNEADref1z−BURN1ref1z)
− (Θsign−IGNEADref1-BURN1ref1) (69)
In step 284, the current value θmb31p [degBTDC] of the crank angle when the combustion mass ratio is 31% is calculated by the following equation.

θmb31p = θmb31pz−θ2pdif (70)
In step 285, the volume Vmb31p [m ³ ] at the combustion mass ratio 31% in the combustion chamber 5 is calculated from the crank angle θmb31p by the following equation. This operation is the same as the operation in step 242 of FIG.

Vmb31p = f25 (θmb31p) (71)
In step 286, the effective compression ratio Emb31p at the combustion mass ratio 31% is calculated from the volume Vmb31p at the combustion mass ratio 31% of the combustion chamber 5 by the following equation.

Emb31p = Vmb31p / VIVC (72)
The operation for calculating TC and PC in steps 287 to 291 is the same as the operation for calculating TC and PC in steps 244 to 248 in FIG. The difference is that in steps 247 and 248 of FIG. 34, Emb1p, which is the compression ratio when the combustion mass ratio is 1%, and BR1p, which is the combustion mass ratio when the combustion mass ratio is 1%, are used. Steps 290 and 291 use Emb31p which is a compression ratio at a combustion mass ratio of 31% and BR31p which is a combustion mass ratio at a combustion mass ratio of 31%. That is, in steps 290 and 291 of FIG. 36, the average temperature TC [K] and the average pressure PC [Pa] when the fuel in the combustion chamber 5 burns are calculated by the following equations.

TC = TC0 × Emb31p ^ 0.35
+ CF # × QINJ × BR31p / (MASSZ + WIDRY + QINJ)
... (73)
PC = PC0 * Emb31p ^ 1.35 * TC / TC0 / Emb31p ^ 0.35
... (74)
Where CF #: lower heating value of fuel,
The operation for calculating the first main combustion period in steps 292 to 295 is the same as the operation for calculating the main combustion period in steps 188 to 191 in FIG. The difference is that the temperature TTDC and the pressure PTDC at the compression top dead center timing of the combustion chamber 5 are used in step 188 in FIG. 12, whereas in step 292 in FIG. 36, the above equations (73) and (74) are used. The point at which the obtained temperature TC and pressure PC are used, and the combustion mass ratio BR2 (= 58%) at the end timing of the main combustion period BURN2 in step 191 of FIG. Whereas the volume VTDC is used, in step 295 of FIG. 36, the combustion mass ratio is 31% (= BR31p) in the vicinity of the intermediate period of the main combustion period, and the combustion mass ratio of the combustion chamber 5 is 31%. It is a point using Vmb31p which is the volume in time. That is, in step 292, the laminar combustion speed SL2 [m / sec] in the first main combustion period is calculated, and in step 295, the first main combustion period BURN2ref1 [deg] is calculated by the following equations.

SL2 = SLstd × (TC / Tstd) ^2.18 × (PC / Pstd) ^−0.16 (75)
Where Tstd: reference temperature [K],
Pstd: Reference pressure [Pa]
SLstd: Reference laminar combustion at reference temperature Tstd and reference pressure Pstd
Speed [m / sec],
BURN2ref1 = {(NRPM × 6) × BR31p × Vmb31p}
/ (RPROBA × AF2 × FLAME1) (76)
However, AF2: Reaction area (fixed value) of flame kernel [m ² ],
When the calculation of the first initial combustion period BURN1ref1, the first ignition dead time equivalent crank angle IGNDEADref1, and the first main combustion period BURN2ref1 is thus completed, the process returns to step 227 of FIG. Formula X1 ′ is calculated. The calculation of the first approximate expression will be described with reference to the flowchart of FIG.

  Here, before explaining the flow of FIG. 37, how the approximate expression of the combustion mass ratio is obtained will be described first.

  In general, the combustion mass ratio X [%] is often described by the following equation as an exponential function with the crank angle θ [deg] as a variable.

X = 1−exp [−a · {(θ−θs) / θB} ^ (m + 1)]
... (Supplement 4)
Where θs: heat generation start crank angle,
θB: crank angle corresponding to the combustion period,
a, m: constants that determine the profile,
The constants a and m in (Supplement 4) are values uniquely determined by the engine model, such as the shape of the combustion chamber, the position of the spark plug, and the gas flow characteristics in the combustion chamber. Therefore, the constants a and m can be determined in advance by experiments or the like. (Supplement 4) The expression is called a Wiebe function.

  In the present embodiment, the constants a and m are calculated based on predetermined parameters, and an expression obtained by substituting the calculated constants a and m into the above (Supplement 4) is obtained as an approximate expression of the combustion mass ratio. More specifically, the (complement 4) equation is modified as follows.

X = 1−exp [− {a / θB ^ (m + 1)} · (θ−θs) ^ (m + 1)]
... (Supplement 5)
In this (Supplement 5), the constant k is
k≡a / {θB ^ (m + 1)} (Supplement 6)
The following equation is obtained in place of the (complement 4) equation.

X = 1−exp [−k · (θ−θs) ^ (m + 1)] (Supplement 7)
Now, consider applying the Wiebe function of (Supplement 4) to this embodiment.

  First, the heat generation start crank angle θs [degBTDC] is given by the following equation in the present embodiment.

θs = θsign [degBTDC] −IGNEAD [deg]
... (Supplement 8)
If the crank angle when the combustion mass ratio is 2% is θ2 and the crank angle when the combustion mass ratio is 60% is θ60, these are substituted into the above (Supplement 4) equation and the following two equations Is established.

0.02 = 1 − [− a · {(θ2−θs) / θB} ^ (m + 1)]
... (Supplement 9)
0.6 = 1 − [− a · {(θ60−θs) / θB} ^ (m + 1)]
... (Supplement 10)
(Supplement 9) and (Supplement 10) are modified respectively,
−a · {(θ2−θs) / θB} ^ (m + 1) = ln (1−0.02)
= Ln (0.98)
(Supplement 11)
-A · {(θ60−θs) / θB} ^ (m + 1) = ln (1-0.6)
= Ln (0.4)
... (Supplement 12)
Divide both sides of (Supplement 11) and (Supplement 12) respectively,
{(Θ2−θs) / (θ60−θs)} ^ (m + 1)
= Ln (0.98) / ln (0.4)
... (Supplement 13)
(Supplement 13) The following equation is obtained by arranging the equation with respect to the constant m.

m = ln {ln (0.98) / ln (0.4)}
/ Ln {(θ2−θs) / (θ60−θs)} − 1
... (Supplement 14)
Here, “ln” in (Supplement 14) represents a natural logarithm.

  Further, when the above (Supplement 11) equation is modified, the following equation is obtained.

a / {θB ^ (m + 1)} = − ln (0.98) / {(θ2−θs) ^ (m + 1)}
... (Supplement 15)
From (Supplement 15) and (Supplement 6), the constant k is expressed by the following equation using the constant m.

k = a / {θB ^ (m + 1)}
= −ln (0.98) / {(θ2−θs) ^ (m + 1)}
... (Supplement 16)
To summarize the above, if the heat generation start crank angle θs obtained by the above (Supplement 8) is substituted into the above (Supplement 14), the constant m is obtained, and the obtained constant m is substituted into the (Supplement 16) equation. Then, the constant k can be obtained. When the constant m and the constant k are obtained in this way, an equation in which the obtained constants m and k are substituted into the above (Supplement 7) equation can be determined as an approximate equation for the combustion mass ratio.

  FIG. 37 (subroutine of step 227 in FIG. 33A) is for calculating the first approximate expression X1 'of the combustion mass ratio.

  In FIG. 37, in step 301, the ignition timing θsign [degBTDC] obtained in step 222 in FIG. 33A, the first ignition dead time equivalent crank angle IGNDEADref1 [deg] calculated in step 225 in FIG. ].

  In step 302, the heat generation start crank angle θs [degBTDC] is calculated by the following equation.

θs = θsign−IGNDEADref1 (77)
Equation (77) is the same as the above (Supplement 8).

  In step 303, the constant m is calculated by the following equation using the heat generation start crank angle θs.

m = ln {ln (0.98) / ln (0.4)}
/ Ln {(θ2−θs) / (θ60−θs)} − 1 (78)
Where θ2 is the crank angle when the combustion mass ratio is 2%,
θ60: crank angle when the combustion mass ratio is 60%,
The expression (78) is the same as the above (complement 14) expression. Here, θ2 and θ60 in the equation (78) are obtained from the basic ignition timing MBTCAL. Or set by conformity.

  In step 304, another constant k is calculated by the following equation using the constant m thus calculated.

k = a / {θB ^ (m + 1)}
= −ln (0.98) / {(θ2−θs) ^ (m + 1)} (79)
The expression (79) is the same as the above (complement 16) expression.

  In step 305, the first approximate expression X1 '[%] of the combustion mass ratio is calculated by the following expression using the constant k and the constant m described above.

X1 ′ = 1−exp [−k · (θ−θs) ^ (m + 1)] (80)
Expression (80) is the same as the above (complement 7).

  When the calculation of the first approximate expression X1 'for the combustion mass ratio is thus completed, the process returns to step 228 of FIG. 33A to calculate the first illustrated positive work estimated value Wpos1. The calculation of the first illustrated positive work estimated value Wpos1 will be described with reference to the flowcharts of FIGS.

  FIG. 38 and FIG. 39 (subroutine of step 228 in FIG. 33A) are for calculating the first positive work estimated value Wpos1.

38, in step 311, the ignition timing θign [degBTDC] obtained in step 222 in FIG. 33A, the first ignition dead time equivalent crank angle IGNDEADref1 [deg] calculated in step 225 in FIG. ], The fuel amount QINJ [kg], the volume VIVC [m ³ ] at the closing timing of the intake valve of the combustion chamber 5 calculated in step 12 of FIG. 5, and the intake air of the combustion chamber 5 calculated in step 13 of FIG. The temperature TINI [K] at the valve closing timing, the pressure PINI [Pa] at the intake valve closing timing calculated in step 14 of FIG. 5, and the combustion chamber wall temperature Twall [K] detected by the temperature sensor 48 The target equivalent ratio TFBYA is read.

  In step 312, the intake valve closing timing IVC [degBTDC] is set to the crank angle previous value θ (old) [degBTDC].

  In step 313, a value obtained by advancing the crank angle previous value θ (old) by a predetermined angle const2 [deg] is set as the crank angle current value θ (new). Here, the predetermined angle const2 is a value that determines a cycle for calculating the instantaneous work dW described later. If the value of const2 is increased too much, the calculation time of the illustrated positive work estimated value decreases although the calculation time is shortened. On the other hand, if the value of const2 is made too small, the calculation accuracy increases, but the calculation time is prolonged. Therefore, an optimal value is set by conformance. Here, the predetermined angle const2 is considered to be about 1 deg.

In step 314, the crank angle current value θ (new) is compared with the exhaust valve opening timing EVO. If the current crank angle value θ (new) has not reached the exhaust valve opening timing EVO, the routine proceeds to step 315, where the volume Vcyl (new) [m ³ ] at the current crank angle value θ (new) of the combustion chamber 5 is calculated. For this calculation, the piston displacement [m] is first calculated from the current crank angle value θ (new) by the following equation.

Piston displacement = (Crank length + Connecting rod length)-Crank length · cos (θ (new))
-Sqrt {connector length ^ 2-crank length ^ 2 · sin ((θ (new) ^ 2)}
... (81)
The crank length [m] and the connecting rod length [m] in equation (81) are values determined by engine specifications. “Sqrt” in the equation (81) represents taking a square root. However, for the sake of simplicity, an expression in the case where the rotation center of the crankshaft is not offset from the center axis of the cylinder is shown here.

The volume Vcyl (new) [m ³ ] at the crank angle current value θ (new) of the combustion chamber 5 is calculated by the following equation using the piston displacement calculated by the above equation (81).

Vcyl (new) = (π / 4) · D ^ 2 · Piston displacement + Vc (82)
Where D: cylinder bore diameter [m]
Vc: gap volume [m ³ ],
In step 316, the pressure P (new) and temperature T (new) at the crank angle current value θ (new) of the combustion chamber 5 are calculated by the following equations.

P (new) = P (old) × (Vcyl (old) / Vcyl (new)) ^ κ (83)
T (new) = T (old) × (Vcyl (old) / Vcyl (new)) ^ (κ−1)
... (84)
Expressions (83) and (84) are expressions expressing pressure and temperature when the gas is adiabatically compressed in the combustion chamber 5 in accordance with the laws of thermodynamics.

  Here, Vcyl (old) in the equations (83) and (84) is the previous volume value at the previous crank angle value θ (old) of the combustion chamber 5, and P (old) and T (old) are the cranks of the combustion chamber 5. The previous pressure value and the previous temperature value at the previous angle value θ (old) are obtained by step 329 in FIG. 39 described later. As an initial value of P (old), the pressure PINI at the closing timing of the intake valve of the combustion chamber 5 calculated in step 14 of FIG. 5 is calculated as an initial value of T (old) in step 13 of FIG. As the initial value of Vcyl (old), the volume VIVC of the combustion chamber 5 at the intake valve closing timing calculated in step 12 of FIG.

  When the piston 6 rises and the gas is adiabatically compressed in the combustion chamber 5, Vcyl (old) / Vcyl (new)> 1.0. Therefore, from the above equations (83) and (84), the crank angle at the previous value is obtained. The pressure P (new) at the current crank angle value is higher than the previous pressure value P (old), and the temperature T (new) at the current crank angle value is higher than the previous temperature value T (old) at the previous crank angle value. Is higher.

  In step 317, the crank angle current value θ (new) is substituted into the first approximate expression X1 ′ of the combustion mass ratio obtained in step 227 of FIG. 33A, that is, the crank angle current value θ (new). The combustion mass ratio MB (new) is calculated by the following equation.

MB (new) = 1−exp [−k · (θ (new)
-([Theta] sign-IGNEADref1)) ^ (m + 1)]
... (85)
In step 318, the combustion mass ratio MB (new) and the predetermined value const3 are compared, and the crank angle current value θ (new) and the heat generation start crank angle (θign−IGNEADref1) are compared. Here, since combustion ends when the combustion mass ratio reaches 100%, the predetermined value const3 is set to a value sufficiently close to 100%. That is, when the combustion mass ratio MB (new) becomes equal to or greater than the predetermined value const3, it indicates that the combustion is finished.

  When the crank angle current value θ (new) has not reached the heat generation start crank angle (θign−IGNEADref1), no heat is generated yet, and when the combustion mass ratio MB (new) is equal to or greater than the predetermined value const3, the combustion is finished. Therefore, in these cases, the routine proceeds to step 320, where the generated heat amount dQburn from the previous crank angle value to the current crank angle value is set to zero.

  On the other hand, when the crank angle current value θ (new) exceeds the heat generation start crank angle (θsign−IGNEADrefref1) and the combustion mass ratio MB (new) is less than the predetermined value const3, heat is generated and combustion ends. If it is determined that it is not, the routine proceeds from step 318 to step 319, where the amount of heat generated dQburn [kJ] from the previous crank angle value to the current crank angle value is calculated by the following equation.

dQburn = (MB (new) −MB (old)) · QINJ · CF # (86)
Where CF #: lower heating value of fuel [kJ / kg],
Here, MB (old) in equation (86) is the combustion mass ratio previous value at the crank angle previous value θ (old), and is obtained by step 329 of FIG.

  In step 321, a cooling loss dQloss [kJ] from the previous crank angle value to the current crank angle value is calculated by the following equation. This calculation will be described with reference to the flowchart of FIG.

In FIG. 40 (subroutine of step 321 in FIG. 38), in step 341, the heat transfer coefficient h [kcal / (m ² · h · K)] is calculated by the following equation.

h = 110 · D ^ (− 0.2) · P (old) ^ 0.8 · T (old) ^ (− 0.53)
・ [C1 ・ Cm + C2 ・ (Vcyl (old) ・ TINI)
/(PINI·VIVC)·(P(old)-PM(old))]^0.8
... (87)
Where D: cylinder bore diameter [m]
Cm: average piston speed [m / s],
PM (old): previous value of motoring pressure [Pa]
Here, PM (old) in the equation (87) is the previous pressure value at the time of motoring at the previous crank angle value θ (old) of the combustion chamber 5 and is obtained in step 344 described later. As an initial value of PM (old), the collector internal pressure PCOL detected by the pressure sensor 44 is entered. In the above equation (87), P (old), T (old), and Vcyl (old) are the previous pressure value, the previous temperature value, and the previous volume value at the crank angle previous value θ (old) of the combustion chamber 5, respectively. This is obtained by step 329 of FIG. 39 described later.

  The two constants C1 and C2 in the above equation (87) are commonly used constants. For example, one constant C1 is 6.18 in the exhaust stroke and 2.28 in the compression stroke and the expansion stroke. The other constant C2 is 3.24 × 10 ^ (− 3) [m / (s · K)].

  In step 342, a cooling loss dQloss [kJ] from the previous crank angle value to the current crank angle value is calculated by the following equation using the heat transfer coefficient h.

dQloss = h · A · (T (old) −Twall) · Δt (88)
Where A: the surface area of the combustion chamber [m ² ],
Δt: Calculation interval from the previous crank angle value to the current crank angle value [h]
Since the previous temperature value T (old) on the right side of the equation (88) is higher than the combustion chamber wall temperature Twall, the cooling loss dQloss on the left side of the equation (88) is a positive value.

  Δt in equation (88) may be obtained as follows. That is, since the predetermined angle const2 [deg] is the calculation interval from the crank angle previous value to the crank angle current value, this is converted into time [h] using the engine speed NRPM at that time.

  In step 343, the pressure PM (new) at the time of motoring at the crank angle current value θ (new) of the combustion chamber 5 is calculated by the following equation which is the same as the equation (81).

PM (new) = PM (old) × (Vcyl (old) / Vcyl (new)) ^ κ (89)
In step 344, the pressure PM (new) at the time of motoring at the crank angle current value θ (new) of the combustion chamber 5 calculated in this way is prepared for the next time at the crank angle previous value θ (new) of the combustion chamber 5 in preparation. Move to the previous pressure value PM (old) during motoring.

  When the calculation of the cooling loss dQloss from the previous crank angle value to the current crank angle value is thus completed, the process returns to step 322 in FIG. 38 to calculate the constant volume specific heat Cv [J / (kg · K)]. This is done by retrieving a predetermined table from the previous temperature value T (old) at the previous crank angle value of the combustion chamber 5 to obtain the constant pressure specific heat Cp1 and specific heat ratio κ1 of the burned gas, and also the previous crank angle value of the combustion chamber 5. A predetermined table is searched from the previous temperature value T (old) at, and the constant pressure specific heat Cp2 and specific heat ratio κ2 of the unburned gas are obtained.

  Here, by dividing the constant pressure specific heats Cp1 and Cp2 by the corresponding specific heat ratios by κ1 and κ2, the constant volume specific heats Cv1 (= Cp1 / κ1) and Cv2 (= Cp2 / κ2) can be obtained. The burned gas constant volume specific heat Cv1 and the unburned gas constant volume specific heat Cv2 are weighted and averaged by the combustion mass ratio previous value MB (old), and the calculated weighted average value is the gas in the combustion chamber 5. Is set as the constant volume specific heat Cv. That is, the constant volume specific heat Cv is calculated by the following equation.

Cv = (Cp1 / κ1) × MB (old) + (Cp2 / κ2) × (1−MB (old))
... (90)
However, Cp1: constant pressure specific heat of burned gas,
Cp2: constant pressure specific heat of unburned gas,
κ1: Specific heat ratio of burned gas,
κ2: specific heat ratio of unburned gas,
In step 323, the total gas amount GASALL [kg] in the combustion chamber 5 is calculated. This may be calculated from the fuel amount QINJ and the target equivalent ratio TFBYA according to the following equation.

GASALL = QINJ × (1 + 1 / TFBYA) (91)
In step 324 of FIG. 39, the temperature T (new) and the pressure P (new) at the crank angle current value θ (new) of the combustion chamber 5 are calculated again by the following equations.

T (new) = T (old) × (Vcyl (old) / Vcyl (new)) ^ (κ−1)
+ (DQburn-dQloss) / (GASALL · Cv)
... (92)
P (new) = {P (old) × (Vcyl (old)) / Vcyl (new) ^ κ}
T (new) / {T (old) × (Vcyl (old) / Vcyl (new)) ^ (κ−1)}
... (93)
Here, the first term on the right side of the equation (92) is the same as the above equation (84) because the temperature rises when the gas is adiabatically compressed in the combustion chamber 5, and the right side of the equation (92) The second term is the temperature rise due to combustion. Further, {P (old) × (Vcyl (old)) / Vcyl (new) ^ κ} on the right side of the equation (93) is a pressure increase when the gas is adiabatically compressed in the combustion chamber 5. Same as the above equation (83), T (new) / {T (old) × (Vcyl (old) / Vcyl (new)) ^ (κ-1)} on the right side of equation (93) is the temperature by combustion The pressure rises with the rise.

  In step 325, the average work from the previous crank angle value to the current crank angle value is calculated as the instantaneous work dW [kJ] by the following equation.

dW = (P (new) + P (old)) × (Vcyl (new) −Vcyl (old)) / 2
... (94)
Here, when the piston 6 passes the compression top dead center and descends toward the bottom dead center TDC, the volume Vcyl (new) becomes larger than the previous volume value Vcyl (old), and the value on the left side of the equation (94) is positive. On the other hand, when increasing toward the compression top dead center TDC, the volume Vcyl (new) is smaller than the previous volume value Vcyl (old), and the value on the left side of the equation (94) becomes a negative value. The value on the left side of equation 94) can be positive or negative, but here the absolute value is taken so that the instantaneous work dW is always obtained as a positive value.

  In step 326, the crank angle current value θ (new) is compared with the compression top dead center TDC. If the crank angle current value θ (new) does not reach the compression top dead center TDC, it is determined that the work to be performed by the engine is negative work, and the routine proceeds to step 327, where the negative value at the crank angle current value θ (new) is negative. The work Wn (new) is calculated by the following formula.

Wn (new) = Wn (old) + dW (95)
Here, since the instantaneous work dW on the right side of the equation (95) is always calculated as a positive value, the negative work Wn (new) on the left side of the equation (95) is also calculated as a positive value.

  On the other hand, after the crank angle current value θ (new) reaches the compression top dead center TDC, it is determined that the work to be performed by the engine is a positive work, and the process proceeds from step 326 to step 328 to proceed to the crank angle current value θ (new). The positive work Wp (new) at is calculated by the following equation.

Wp (new) = Wp (old) + dW (96)
Wn (old) and Wp (old) on the right side of the above equations (95) and (96) are the negative work previous value at the crank angle previous value θ (old) and the positive value at the crank angle previous value θ (old), respectively. The previous work value is obtained in step 329 described later.

  In step 329, the crank angle current value θ (new), the combustion mass ratio MB (new) at the crank angle current value θ (new), the pressure P (new) at the crank angle current value θ (new) of the combustion chamber 5, and the combustion The temperature T (new) at the crank angle current value θ (new) of the chamber 5, the volume Vcyl (new) at the crank angle current value θ (new) of the combustion chamber 5, and the negative work at the crank angle current value θ (new) Wn (new), positive work wp (new) at the current crank angle value θ (new), the combustion mass at the previous crank angle value θ (old) and the previous crank angle value θ (old) in preparation for the next time The ratio previous value MB (old), the previous pressure value P (old) at the previous crank angle value θ (old) of the combustion chamber 5, the previous temperature value T (old) at the previous crank angle value θ (old) of the combustion chamber 5, The previous volume value Vcyl (old) at the previous crank angle value θ (old) of the combustion chamber 5, the previous negative work value Wn (old) at the previous crank angle value θ (old), and the previous crank angle value θ (old). Of positive work previously transferred to a previous value Wp (old).

  Here, the initial value of the previous crank angle value θ (old) is given in step 312 of FIG. Further, as described above in step 316 of FIG. 38, the initial value of P (old) is the pressure PINI calculated at step 14 of FIG. 5 at the closing timing of the intake valve of the combustion chamber 5 as T (old). As an initial value, the temperature TINI at the intake valve closing timing of the combustion chamber 5 calculated in step 13 of FIG. 5 is used, and as an initial value of Vcyl (old), the temperature of the combustion chamber 5 calculated in step 12 of FIG. The volume VIVC at the intake valve closing timing is set in advance. As the initial values of the remaining previous values MB (old), Wn (old), and Wp (old), zero is set at the start of the processing of the flow of FIG.

  In step 330, as in step 313 in FIG. 38, a value obtained by advancing the crank angle previous value θ (old) by a predetermined angle const2 is set as the crank angle current value θ (new), and then the process returns to step 314 in FIG. The operations in steps 314 to 323 in FIG. 38 and steps 324 to 330 in FIG. 39 are repeated. That is, when the crank angle θ becomes the heat generation start crank angle (θign−IGNEADREFref1), the generated heat amount dQburn and the cooling loss dQloss are calculated for each predetermined crank angle const2, and the calculated generated heat amount dQburn, the cooling loss dQloss and the heat are calculated. The pressure P and temperature T during combustion in the combustion chamber 5 are calculated based on the laws of mechanics, and the instantaneous work dW is calculated for each predetermined crank angle const2 based on the pressure P. In this case, until the compression top dead center TDC, the instantaneous work dW is assumed to be a negative work, and after the compression top dead center TDC, the instantaneous work dW is assumed to be a positive work. Accumulated every const2.

  Eventually, when the crank angle current value θ (new) reaches the exhaust valve opening timing EVO in step 314 of FIG. 38, it is determined that it is the timing for ending work due to the combustion gas in the combustion chamber 5, and from step 314 to step 331. Next, the first illustrated positive work estimated value Wpos1 [kJ] is calculated by the following equation, and the processing in FIGS. 38 and 39 is completed.

Wpos1 = Wp (old) −Wn (old) (97)
When the calculation of the first illustrated positive work estimated value Wpos1 is thus completed, the process returns to step 229 of FIG. 33A to calculate the first shaft torque estimated value Trquest1. This calculation will be described with reference to the flowchart of FIG.

  In FIG. 41 (subroutine of step 229 in FIG. 33A), in step 351, the first illustrated positive work estimated value Wpos1 [kJ] calculated in step 228 of FIG. 33A is detected by the pressure sensor 44. The collector internal pressure PCOL [Pa], engine load, and rotational speed NRPM [rpm] are read.

  In step 352, the pumping loss Wneg1 [kJ] is calculated by the following equation.

Wneg = PCOL × Hx (98)
However, Hx: difference between the maximum value and the minimum value of the distance from the TDC of the piston pin [m]
In step 353, the friction loss Wneg2 [kJ] is calculated. This may be obtained by searching a map having the contents shown in FIG. 42 from the engine load and the rotational speed NRPM. As shown in FIG. 42, the friction loss Wneg2 is a value that increases as the load increases and the rotational speed NRPM increases.

  In step 354, the pumping loss Wneg1 and the friction loss Wneg2 are added to the first illustrated positive work estimated value Wpos1, that is, the first shaft torque estimated value Trquest1 [kJ] is calculated by the following equation.

Trquest1 = Wpos + Wneg1 + Wneg2 (99)
This completes the description of the flow in FIG.

  Next, FIG. 33 (B) (subroutine of step 2 of FIG. 2) is different from FIG. 33 (A) only in the handling of the ignition timing. That is, in FIG. 33B, in step 231, the basic ignition timing MBTCAL [degBTDC] and the ignition timing retard amount previous value RET (old) [deg] calculated in step 43 of FIG. 13 are read. Here, the previous ignition timing retard amount RET (old) is a value obtained in step 6 of FIG. 2 (described later).

  In step 232, the value of the ignition timing θign [degBTDC] is moved to the previous ignition timing value θignz [degBTDC] and the ignition timing retard amount previous value RET (old) is subtracted from the basic ignition timing MBTCAL (that is, from the basic ignition timing MBTCAL). The ignition timing retard amount previous value RET (old) is a value on the retard side) is calculated as the ignition timing θign [degBTDC]. Note that zero is set as the initial value of RET (old).

  In step 234, the second initial combustion period BURN1ref2 [deg], in step 235, the second ignition dead time equivalent crank angle IGNDEADref2 [deg], in step 236, the second main combustion period BURN2ref2 [deg], and in step 237, the combustion mass. In step 238, the second approximate positive work estimate value Wpos2 [kJ] is calculated, and in step 239, the second shaft torque estimated value Trquest2 [kJ] is calculated.

  Here, the calculation of the second initial combustion period BURN1ref2 is performed according to the flow of FIG. 43, the calculation of the second ignition dead time equivalent crank angle IGNDEADref2 is performed according to the flow of FIG. 44, and the calculation of the second main combustion period BURN2ref2 is illustrated in FIG. 46, the calculation of the second approximate expression X2 ′ of the combustion mass ratio is performed according to the flow of FIG. 46, and the calculation of the second illustrated positive work estimated value Wpos2 is calculated according to the flows of FIG. 47 and FIG. The calculation of the value Trquest2 will be described with reference to the flow of FIG. 49 (subroutine of step 239 in FIG. 33B).

  In this case, FIG. 43 (subroutine of step 234 in FIG. 33B), FIG. 44 (FIG. 33B), FIG. 45 (subroutine of step 236 in FIG. 33B), FIG. 46 (FIG. 33B) ) Of FIG. 34, FIG. 35, FIG. 36, FIG. 37, FIG. 38 and FIG. 39 are different from those of FIG. 34, ignition timing previous value θignz in step 241; ignition timing θign in step 261 in FIG. 35; ignition timing θignn in step 281 in FIG. 36; ignition timing previous value θignz; ignition timing θignn in step 301 in FIG. The ignition timing θign in step 311 of 38 is all the basic ignition timing MBTCAL, whereas the step of FIG. The ignition timing previous value θignz at 361, the ignition timing θignn at step 381 in FIG. 44, the ignition timing θignn at step 401 in FIG. 45, the previous ignition timing value θignz, the ignition timing θignn at step 421 in FIG. 46, and the step 431 in FIG. The ignition timing θign is a value obtained by retarding the ignition timing retard amount by the previous value RET (old) from the basic ignition timing MBTCAL, and the others are the same, and thus the description thereof is omitted.

  Thus, in FIG. 33 (B), the initial combustion period, the ignition dead time equivalent crank angle, the main combustion period, the approximate expression of the combustion mass ratio, the indicated positive work estimated value, and the estimated shaft torque value are shown in FIG. 33 (A). Separately calculated are the initial combustion period, ignition dead time equivalent crank angle, main combustion period, combustion mass when ignition is performed by the ignition timing retarded amount RET (old) from the basic ignition timing MBTCAL. The approximate expression of the ratio, the estimated positive work estimated value, and the estimated shaft torque value are the initial combustion period when the ignition is performed at the basic ignition timing MBTCAL, the crank angle equivalent to the ignition dead time, the main combustion period, the approximate expression of the combustion mass ratio, the illustrated positive This is because the work estimated value and the shaft torque estimated value are different values.

  When the calculation of the first axis torque estimated value Trquest1 and the second axis torque estimated value Trquest2 is completed in this way, the requested axis torque Trqtgt is calculated in step 3 returning to FIG. The required shaft torque Trqtgt can be calculated by searching a map having the content shown in FIG. 50 from the accelerator opening (accelerator pedal depression amount) APO detected by the accelerator sensor 42 and the rotational speed NRPM. . As shown in FIG. 50, the required shaft torque Trqtgt is a value that increases as the accelerator opening APO increases and the rotational speed NRPM increases.

  In step 4, the ignition timing retard amount RET (new) is calculated. The calculation of the ignition timing retard amount RET (new) will be described with reference to the flowchart of FIG.

  51 (subroutine of step 4 in FIG. 2), in step 471, the first shaft torque estimated value Trquest1 [kJ] calculated in step 2 of FIG. 2, and the required shaft torque Trqtgt calculated in step 3 of FIG. [KJ] is read.

  In step 472, the ignition timing retard amount RET (new) [deg] is calculated based on the difference between the first shaft torque estimated value Trquest1 and the required shaft torque Trqtgt, that is, the following equation.

RET (new) = (Trquest1-Trqtgt) × const201
... (100)
However, const201: predetermined value,
This ignition timing retard amount RET (new) is calculated when the first shaft torque estimated value Trquest1 that is the estimated shaft torque at the basic ignition timing MBTCAL is larger than the required shaft torque Trqtgt, that is, when the shaft torque is excessive. Is retarded from the basic ignition timing MBTCAL so that the required shaft torque can be obtained. Therefore, according to equation (100), the first shaft torque estimated value Trquest1 may be smaller than the required shaft torque Trqtgt, but this is not considered.

  The predetermined value const 201 is a value for converting the shaft torque into the ignition timing retard amount, and is adapted in advance.

  In step 473, a shaft torque estimated value is calculated as a third shaft torque estimated value Trquest3 when it is assumed that the ignition timing retard amount RET (new) is retarded from the basic ignition timing MBTCAL.

  The calculation of the third axis torque estimated value Trquest3 will be described with reference to the flowchart of FIG. In FIG. 52, the same steps as those in FIG.

  Here, the calculation method itself of the third axis torque estimation value Trquest3 is the same as the calculation method of the second axis torque estimation value Trquest2 described with reference to FIG. In FIG. 52, the difference from the case of FIG. 33B is that the ignition timing retard amount previous value RET (old) was used in FIG. 33B, whereas in FIG. 52 it was calculated in step 472 of FIG. The only point is to use the ignition timing retard amount RET (new). That is, in step 481, the basic ignition timing MBTCAL [degBTDC] calculated in step 43 of FIG. 13 and the ignition timing retard amount RET (new) [deg] calculated in step 472 of FIG. The value of the ignition timing θign [degBTDC] is shifted to the previous ignition timing value θignz [degBTDC], and the value obtained by subtracting the ignition timing retard amount RET (new) from the basic ignition timing MBTCAL (that is, the ignition timing retard amount RET from the basic ignition timing MBTCAL). (New value on the retard side) is calculated as the ignition timing θign [degBTDC].

  In Steps 233 to 239 of FIG. 52, the second axis torque estimated value Trquest2 is calculated in the same manner as Steps 233 to 239 of FIG. 33B. In Step 483 of FIG. 52, the calculated second axis torque estimated value Trquest2 is calculated. The process shifts to the triaxial torque estimated value Trquest3.

  When the calculation of the third axis torque estimated value Trqest3 is thus completed, the process returns to step 474 in FIG. 51 to subtract the required axis torque Trqtgt from the third axis torque estimated value Trquest3 to obtain the axis torque difference dTrq (= Trqest3-Trqtgt). In step 475, the shaft torque difference dTrq and the predetermined value const202 are compared.

  Here, the predetermined value const 202 is a value that defines an allowable range of the required shaft torque. If the shaft torque difference dTrq is equal to or greater than the predetermined value const202, the third shaft torque estimated value Trquest3 is not within the allowable range of the required shaft torque Trqtgt. In other words, according to the third shaft torque estimated value Trquest3 at that time, the shaft torque is excessive. In step 476, the ignition timing retard amount RET (new) is moved to the ignition timing retard amount previous value RET (old). In step 477, the ignition timing retard amount RET (new) [deg] is set to a predetermined value. Just make it bigger. That is, the ignition timing retard amount RET (new) is updated by the following equation.

RET (new) = RET (old) + (Trquest3-Trqtgt)
× const201
... (101)
However, RET (new): ignition timing retard amount after update,
RET (old): ignition timing retard amount before update,
const201: the degree of update to the larger side,
The ignition timing retard amount RET (new) is updated to the larger ignition timing when the basic ignition timing MBTCAL is retarded by a predetermined ignition timing retard amount (retard amount in the above equation (100)). The third shaft torque estimated value Trquest1 that is the estimated shaft torque is also larger than the required shaft torque Trqtgt, that is, when the shaft torque is excessive, the ignition timing retard amount is further increased to obtain the required shaft torque. Is for. Therefore, according to the equation (101), the third shaft torque estimated value Trquest3 may be smaller than the required shaft torque Trqtgt, but this case is not considered.

  The predetermined value const201 in the equation (101) is the same as the equation (100).

  Thereafter, the process returns to step 473 and the operations of steps 473 to 477 are repeated. By this operation, when the third shaft torque estimated value Trquest3 decreases and the shaft torque difference dTrq becomes less than the predetermined value const202, the process proceeds from step 475 to “END”, and the process of FIG. In this case, the operations in steps 473 to 477 are performed without actually igniting, and are completed almost instantaneously.

  On the other hand, when the shaft torque difference dTrq is less than the predetermined value const202 in Step 475, that is, when the third shaft torque estimated value Trquest3 is within the allowable range of the requested shaft torque Trqtgt from the beginning, “END” is also obtained from Step 475. Then, the process of FIG. 51 is terminated.

  When the calculation of the ignition timing retard amount RET (new) is completed in this way, the process returns to step 5 in FIG. 2, and the ignition timing retard amount RET (new) is retarded from the basic ignition timing MBTCAL. Command value QADV [degBTDC] is set.

QADV = MBTCAL-RET (new) (102)
In step 6, the ignition timing retard amount RET (new) is moved to the ignition timing retard amount previous value RET (old). The value of the ignition timing retard amount previous value RET (old) is used in steps 231 and 232 of FIG.

  The ignition timing command value QADV set by the above equation (102) is transferred to the ignition register in step 7, and ignition is performed to cut off the primary current from the engine controller 31 at the timing when the actual crank angle coincides with the ignition timing command value QADV. A signal is output to the ignition coil 13.

  Next, FIG. 53 is for calculating the required air increase rate Rqai as the engine torque correction amount, and is executed for each fuel injection.

  Here, the reason why the required air increase rate Rqai is introduced is as follows. That is, as described above with reference to FIG. 51, the ignition timing is retarded when the first shaft torque estimated value Trquest1 is greater than the required shaft torque Trqtgt (see step 472 in FIG. 51), so the ignition retarded from the basic ignition timing MBTCAL. The shaft torque generated when ignition is performed at the timing (that is, the second shaft torque estimated value Trquest2) is lower than the shaft torque generated when ignition is performed at the basic ignition timing MBTCAL (that is, the first shaft torque estimated value Trquest1). In order to avoid a decrease in shaft torque accompanying this ignition timing retard, torque compensation is performed.

  First, in step 491, the first axis torque estimated value Trquest1 [kJ] and the second axis torque estimated value Trquest2 [kJ] calculated in step 2 of FIG. 2 are read.

  In step 492, the shaft torque shortage amount Ts [kJ] from the first shaft torque estimated value Trquest1 is calculated by the following equation.

Ts = Trquest1-Trquest2 (103)
In step 493, the shaft torque shortage amount Ts is compared with a predetermined value const203. Here, the predetermined value const 203 defines an allowable range of the shaft torque shortage and is set in advance. When the shaft torque shortage amount Ts is equal to or greater than the predetermined value const203, the process proceeds from step 493 to step 494, where the required air increase rate Rqai (new) is increased by the predetermined value const204. That is, the required air increase rate Rqai (new) is updated by the following equation.

Rqai (new) = Rqai (old) + const204 (104)
However, Rqai (new): Requested air increase rate after update,
Rqai (old): required air increase rate before update,
const204: Amount of update to the side to be increased (positive value),
Here, 1.0 is set as the initial value of Rqai (old). When the shaft torque shortage amount Ts becomes equal to or greater than the predetermined value const203 for the first time after the engine is operated, the required air increase rate Rqai (new) becomes a value exceeding 1.0 because Rqai (old) = 1.0.

  In step 495, the value of Rqai (new) is moved to Rqai (old) representing the previous value in preparation for the next time, and the current process is terminated.

  In the present embodiment, the same control as described in Japanese Patent Laid-Open No. 9-287513 is performed for the fuel amount and the air amount. The control of the fuel amount and the air amount will be described. The reference target intake air amount tTP for obtaining the theoretical air-fuel ratio based on the accelerator opening APS detected by the accelerator sensor 42 and the engine speed NRPM, the stoichiometric air-fuel ratio. Target equivalent ratio tDML for enabling operation at a lean air-fuel ratio or rich air-fuel ratio, and target intake air for obtaining the target air-fuel ratio from these reference target intake air amount tTP and target equivalent ratio tDML An amount tTP ′ (= tTP / tDML) is calculated, a target throttle valve opening tTPS is calculated based on the target intake air amount tTP ′, and the throttle motor 24 is configured to realize the target throttle valve opening tTPS. The actual intake air detected by the engine speed and the air flow meter 32 while controlling the air amount by controlling the opening of the throttle valve 23 By giving the fuel injector 21 the fuel injection pulse width TI [ms] obtained by correcting the basic fuel injection pulse width TP [ms] calculated based on the above with the target equivalent ratio tDML and the invalid pulse width TS. The amount is controlled.

  In the present embodiment, on the assumption that the fuel amount and the air amount are controlled, the required air increase rate Rqai updated in step 494 (or step 496) in FIG. 53 to the target intake air amount (tTP or tTP ′). The value obtained by multiplying (new) is used as the target intake air amount (tTP or tTP ′) to control the air amount. Since the required air increase rate Rqai (new) updated at step 494 in FIG. 53 is a value exceeding 1.0, the required air increase rate Rqai (new) exceeding 1.0 is the target intake air amount. Is multiplied, the target air amount after multiplication becomes larger than before multiplication, and the opening of the throttle valve 23 is made larger than before multiplication so that the target air amount after multiplication is obtained. As a result, the amount of intake air is increased, the amount of fuel is increased through the increase in the amount of air, and the shaft torque generated by the engine is increased even at the same ignition timing. Then, in parallel with the increase in the target air amount based on the required air increase rate Rqai (new), the second shaft torque estimated value Trquest2 is calculated in step 2 of FIG. 2, and the target air based on the required air increase rate Rqai (new) is calculated. If the shaft torque increases as the amount increases, the calculated second shaft torque estimated value Trquest2 also increases, and the shaft torque shortage amount Ts calculated in step 492 of FIG. 53 decreases. If the shaft torque shortage amount Ts is still greater than or equal to the predetermined value const203, the required air increase rate Rqai (new) is updated to a larger side in step 494 of FIG. As a result, the shaft torque generated by the engine increases and the shaft torque shortage amount Ts decreases, so that the shaft torque shortage amount Ts eventually falls below the predetermined value const203 by repeating step 494 of FIG. At this time, in order to maintain the value of the required air increase rate Rqai (new) at that time, the routine proceeds from step 493 to step 496 to move the previous required air increase rate value Rqai (old) to the required increase rate current value Rqai (new). This process is terminated.

  Here, the effect of this embodiment is demonstrated.

  According to the present embodiment (the invention described in claim 7), the operating state of the engine is detected, and based on the detected operating state, the combustion mass ratio (X1 ′, X2 ′), the shaft torque generated by the engine is estimated based on the estimated combustion mass ratio, and torque control (intake air amount increase correction and ignition is performed based on the estimated shaft torque (Trquest1, Trquest2). (Time retard correction), the shaft torque generated by the engine can be estimated without directly detecting the pressure in the combustion chamber by the pressure sensor.

  According to the present embodiment (the invention described in claims 2 and 8), the axial torque estimating means calculates the approximate expression (X1 ′, X2 ′) of the combustion mass ratio at the crank angle during combustion. (See Step 227 in FIG. 33A and Step 237 in FIG. 33B) and approximate positive work estimated values (Wpos1, Wpos2) based on approximate equations (X1 ′, X2 ′) of the combustion mass ratio. The figured positive work estimated value calculation means for calculating (see step 228 in FIG. 33A, step 238 in FIG. 33B), and the pumping loss Wneg1 and friction loss Wneg2 to the figured positive work estimated values (Wpos1, Wpos2). Is calculated as a shaft torque estimated value (Trquest1, Trquest2) (step 229 in FIG. 33A, step 229 in FIG. 33B). Because it contains 239 references) and can be calculated shown positive work estimate the (Wpos1, Wpos2) accurately.

  According to the present embodiment (the invention according to claims 3 and 9), the engine operating state is the heat generation start time θs, and the approximate expression (X1 ′, X2 ′) of the combustion mass ratio is the heat generation start time. Since it is an expression of an exponential function using θs and two constants (m, k) (see step 305 in FIG. 37 and step 425 in FIG. 46), the change in the combustion mass ratio from the start of combustion is given almost strictly. be able to.

  According to the present embodiment (inventions according to claims 4 and 10), when ignition is performed at the basic ignition timing MBTCAL (ignition timing at which MBT is obtained), the ignition timing retard amount previous value RET (old) from the basic ignition timing MBTCAL. Therefore, the approximate expression of the combustion mass ratio is different from that when ignition is performed at a retard side value (ignition timing retarded from the ignition timing at which MBT is obtained) (steps 222, 227, and 37 in FIG. 33A). Step 305, Steps 232 and 237 in FIG. 33 (B), Step 425 in FIG. 46), the first shaft torque estimated value Trquest1 (the shaft torque estimated value when ignition is performed at the ignition timing at which MBT is obtained), the second The accuracy of the shaft torque estimated value Trquest2 (the shaft torque estimated value when ignition is performed at an ignition timing retarded from the ignition timing at which MBT is obtained) is accurate. Calculation can be.

  According to the present embodiment (the invention described in claim 13), the shaft torque estimated value calculating means calculates the shaft torque estimated value when ignition is performed at the ignition timing at which MBT is obtained as the first shaft torque estimated value Trquest1. The first shaft torque estimated value calculation means (see step 229 in FIG. 33A) and the shaft torque estimated value when ignition is performed at an ignition timing retarded from the ignition timing obtained by MBT, are calculated as the second shaft torque estimated value Trquest2. Second torque estimated value calculation means (see step 239 in FIG. 33B), and the torque control means matches the second axis torque estimated value Trquest2 with the first axis torque estimated value Trquest1. Thus, since the air amount increasing means for increasing the intake air amount (see FIG. 52), the ignition timing is retarded from the ignition timing at which MBT is obtained (step in FIG. 2). 4 and 5), even if a decrease in the shaft torque not intended by the driver occurs, the same shaft torque as the first shaft torque estimated value Trquest1 (the shaft torque when ignition is performed at the ignition timing at which MBT is obtained) Is obtained.

  When the first shaft torque estimated value Trquest1 (the shaft torque estimated value when ignition is performed at the ignition timing at which MBT is obtained) exceeds the required shaft torque Trqtgt, the shaft torque is excessively generated. According to the fourteenth aspect of the present invention, when the first shaft torque estimated value Trquest1 exceeds the required shaft torque Trqtgt, the torque control means makes the first shaft torque estimated value Trquest1 coincide with the required shaft torque Trqtgt. Includes an ignition timing correction means (see steps 471 and 472 in FIG. 51) for correcting the ignition timing to the retard side by a predetermined ignition timing retard amount RET (new), so that the first shaft torque estimated value Trquest1 is the required shaft torque Trqtgt. The shaft torque as intended by the driver can be generated even when

  When the shaft torque when the ignition timing is retarded by the predetermined ignition timing retard amount RET (new) exceeds the required shaft torque Trqtgt, the shaft torque is still excessively generated. 15), it is assumed that the shaft torque estimated value calculation means retards the ignition timing by the predetermined ignition timing retard amount RET (new) from the basic ignition timing MBTCAL (ignition timing at which MBT is obtained). Including a third axis torque estimated value calculating means (see step 473 of FIG. 51, FIG. 52) for calculating the estimated shaft torque value as the third axis torque estimated value Trquest3. When the estimated value Trqest3 exceeds the required shaft torque Trqtgt, the third shaft torque estimated value Trqest3 becomes the required shaft torque Trqtgt. Since ignition timing correction means (see steps 473 to 477 in FIG. 51) for correcting the ignition timing to the retard side so as to match is included, even when the third shaft torque estimated value Trquest3 still exceeds the required shaft torque Trqtgt. It is possible to generate shaft torque as intended by the driver.

  In the embodiment, as shown in FIG. 51, a shaft torque estimated value is obtained when it is assumed that the first shaft torque estimated value Trquest1 is retarded by an ignition timing retard amount for reducing the first shaft torque estimated value Trquest1 to the required shaft torque Trqtgt. In the case where the estimated torque value is still larger than the required shaft torque Trqtgt, the ignition timing retard amount is updated to the larger side, but a mode in which the operations in steps 473 to 477 in FIG. 51 are omitted can be considered.

  In the embodiment, the engine operating state is the heat generation start time θs, and the approximate expression of the combustion mass ratio is an exponential function expression using the heat generation start time θs and two constants (m, k). As described above, the engine operating state is a combustion period from the start of combustion to a predetermined crank angle, and is provided with combustion period calculation means for dividing the combustion period into a plurality of calculation periods, and based on the plurality of combustion periods. An approximate expression of the combustion mass ratio can also be set (second embodiment).

  The second embodiment will now be described. As a crank angle for obtaining the combustion mass ratio, a crank having a compression top dead center TDC of zero as a reference, a retarded side is added to this, and an advanced side is set to a minus. Consider the angle Θ [degATDC]. At this time, the combustion mass ratio BR when ignited at the basic ignition timing MBTCAL is as follows using the basic ignition timing MBTCAL, the first ignition dead time equivalent crank angle IGNDEADref1, the first initial combustion period BURN1ref1, and the first main combustion period BURN2ref1. It can be approximated by such a linear expression.

First combustion delay period;
BR = 0 (105)
First initial combustion period;
BR = SS1 × (Θ + MBTCAL-IGNEADREFref1)
... (106)
First main combustion period;
BR = 0.02 + SS2 × (Θ + MBTCAL-IGNEADREFref1
-BURN1ref1)
... (107)
However, SS1: 0.02 / BURN1ref1,
SS2: 0.58 / BURN2ref1,
Similarly, the combustion mass ratio BR when igniting at the ignition timing retarded by the ignition timing retard amount previous value RET (old) from the basic ignition timing MBTCAL is the basic ignition timing MBTCAL, the second ignition dead time equivalent crank angle IGNDEADref2, Using the second initial combustion period BURN1ref2 and the second main combustion period BURN2ref2, it can be approximated by the following linear expression.

Second combustion delay period;
BR = 0 (108)
Second initial combustion period;
BR = SS1 × (Θ + MBTCAL-IGNDEADref2)
... (109)
Second main combustion period;
BR = 0.02 + SS2 × (Θ + MBTCAL−IGNEADREFref2
-BURN1ref2)
... (110)
However, SS1: 0.02 / BURN1ref2,
SS2: 0.58 / BURN2ref2,
FIG. 54 is a linear expression of the combustion mass ratio when ignition is performed at the basic ignition timing MBTCAL and the combustion mass ratio when ignition is performed at an ignition timing retarded by the ignition timing retard amount previous value RET (old) from the basic ignition timing MBTCAL. The approximate waveform is shown as a model.

  Thus, according to the second embodiment (the inventions described in claims 5 and 11), the engine operating state is a combustion period from the start of combustion to a predetermined crank angle, and this combustion period is defined as an initial combustion period and a main combustion period. Combustion period calculation means (see Steps 224 to 226 in FIG. 33A and Steps 234 to 236 in FIG. 33B) for calculating by dividing into two (plural) periods, Since the combustion mass ratio approximation formula (see formulas (105) to (110) above) is set based on BURN1ref1, BURN1ref2, BURN2ref1, BURN2ref2, IGNDEADref1, IGNDEADref2), the change of the combustion mass ratio is approximated by a simple linear formula This makes it easy to set an approximate expression for the combustion mass ratio.

  In the embodiment, the combustion mass ratio (X1 ′, X2 ′) at the crank angle during combustion of the gas in the combustion chamber has been described, but instead, combustion at the crank angle during combustion of the gas in the combustion chamber is described. You may make it estimate the heat generation pattern which is a differential value of mass, this combustion mass, or the said combustion mass ratio.

  In the embodiment, the operating state detecting means for detecting the operating state of the engine, the combustion mass at the crank angle during combustion of the gas in the combustion chamber based on the detected operating state, and at the crank angle during combustion of the gas within the combustion chamber Estimating means for estimating any one of the combustion mass ratio, the heat generation pattern that is a differential value of the combustion mass or the combustion mass ratio, and estimating the shaft torque generated by the engine based on any one of the estimated values In the above description, the description has been given of the case of including the torque control means for performing the torque control based on the estimated shaft torque (the invention according to claim 7). And the combustion mass at the crank angle during combustion of the gas in the combustion chamber, the crank angle during combustion of the gas within the combustion chamber based on the detected operating state Shaft torque by estimating means for estimating any one of the combustion mass ratio and heat generation pattern and shaft torque estimating means for estimating the shaft torque generated by the engine based on any one of the estimated values An arithmetic device can be configured.

  In the first and seventh aspects of the invention, the function of the operating state detecting means is step 302 in FIG. 37, step 422 in FIG. 46 or steps 224 to 226 in FIG. 33 (A), and steps 234 to 226 in FIG. 236, the function of the estimation means is step 227 of FIG. 33A and step 237 of FIG. 33B, and the function of the shaft torque estimation means is steps 227 to 229 of FIG. 33A, FIG. 33B. By the steps 237 to 239, the function of the torque control means is fulfilled by FIG. 51 or FIG. 53, respectively.

The engine control system figure of one Embodiment. The flowchart for demonstrating ignition timing control. The pressure change figure of a combustion chamber. The characteristic view explaining the change of a combustion mass ratio. The flowchart for demonstrating calculation of a physical quantity. Diagram explaining the positional relationship between the crankshaft of the engine and the connecting rod. The characteristic diagram of a water temperature correction coefficient. The characteristic view of an equivalence ratio correction coefficient. The characteristic figure of a standard crank angle. The flowchart for demonstrating calculation of an initial combustion period. The characteristic figure of a temperature rise rate. The flowchart for demonstrating calculation of the main combustion period. The flowchart for demonstrating calculation of basic ignition timing. The flowchart for demonstrating calculation of an internal inert gas rate. The flowchart for demonstrating calculation of an internal inert gas amount. The flowchart for demonstrating calculation of the amount of inert gas at the time of EVC. The flowchart for demonstrating calculation of the amount of inert gas blown back during overlap. The flowchart for demonstrating the setting of a supercharging determination flag and a choke determination flag. The flowchart for demonstrating calculation of the inactive gas flow rate during the overlap when there is no supercharging and there is no choke. The flowchart for demonstrating the calculation of the inactive gas flow rate during the overlap at the time of supercharging without a choke. The flowchart for demonstrating the calculation of the inactive gas flow rate during the overlap at the time of supercharging and no choke. The flowchart for demonstrating calculation of the inactive gas flow rate during the overlap at the time of supercharging and a choke. The characteristic figure of the combustion chamber volume in the exhaust valve closing timing. The characteristic figure of the gas constant of an inert gas. The characteristic figure of the integrated effective area during overlap. Explanatory drawing of the integrated effective area during overlap. The characteristic figure of the specific heat ratio of an inert gas. The characteristic figure of the specific heat ratio of air-fuel | gaseous mixture. The flowchart for demonstrating calculation of a crank angle equivalent to ignition dead time. The characteristic figure of ignition dead time basic value. The characteristic figure of a temperature correction coefficient. The characteristic diagram of an excess air ratio correction coefficient. The flowchart for demonstrating calculation of a 1st axis | shaft torque estimated value. The flowchart for demonstrating calculation of a 2nd axis | shaft torque estimated value. The flowchart for demonstrating calculation of a 1st initial combustion period. The flowchart for demonstrating calculation of the 1st ignition dead time equivalent crank angle. The flowchart for demonstrating calculation of a 1st main combustion period. The flowchart for demonstrating calculation of the 1st approximate expression of a combustion mass ratio. The flowchart for demonstrating calculation of a 1st illustration positive work estimated value. The flowchart for demonstrating calculation of a 1st illustration positive work estimated value. The flowchart for demonstrating calculation of a cooling loss. The flowchart for demonstrating calculation of a 1st axis | shaft torque estimated value. The characteristic diagram of friction loss. The flowchart for demonstrating calculation of a 2nd initial stage combustion period. The flowchart for demonstrating calculation of a 2nd ignition dead time equivalent crank angle. The flowchart for demonstrating calculation of a 2nd main combustion period. The flowchart for demonstrating calculation of the 2nd approximate expression of a combustion mass ratio. The flowchart for demonstrating calculation of the 2nd illustration positive work estimated value. The flowchart for demonstrating calculation of the 2nd illustration positive work estimated value. The flowchart for demonstrating calculation of a 2nd initial stage combustion period. Fig. 4 is a characteristic diagram of required shaft torque. The flowchart for demonstrating calculation of the ignition timing retard amount. The flowchart for demonstrating calculation of a 3rd axis torque estimated value. The flowchart for demonstrating calculation of a request | requirement air increase rate. The wave form diagram which shows the change of a combustion mass ratio.

Explanation of symbols

1 Engine 5 Combustion chamber 11 Ignition device (spark ignition means)
21 Fuel injector 31 Engine controller

Claims (15)

An operating state detecting means for detecting the operating state of the engine;
Estimating means for estimating one of the combustion mass at the crank angle during combustion of the gas in the combustion chamber, the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber, and the heat generation pattern based on the detected operating state When,
A shaft torque calculating device comprising: shaft torque estimating means for estimating a shaft torque generated by the engine based on any one of the estimated values. When the estimation means estimates the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber,
The shaft torque estimating means includes
An approximate expression calculating means for calculating an approximate expression of the combustion mass ratio at the crank angle during combustion;
Illustrated positive work estimated value calculating means for calculating the indicated positive work estimated value based on the approximate expression of the combustion mass ratio,
The shaft torque calculation device according to claim 1, further comprising: a shaft torque estimated value calculation unit that calculates a value obtained by adding a pumping loss and a friction loss to the illustrated positive work estimated value as a shaft torque estimated value.   The operating state of the engine is a heat generation start time, and the approximate expression of the combustion mass ratio is an exponential function expression using the heat generation start time and two constants. Shaft torque calculation device.   The approximate expression of the combustion mass ratio is different between when ignition is performed at an ignition timing at which MBT is obtained and when ignition is performed at an ignition timing retarded from the ignition timing at which MBT is obtained. Shaft torque calculation device.   The operating state of the engine is a combustion period from the start of combustion to a predetermined crank angle, and includes combustion period calculation means for dividing the combustion period into a plurality of calculations, and based on the plurality of combustion periods, the combustion mass ratio is calculated. The shaft torque calculation device according to claim 2, wherein an approximate expression is set.   6. The approximate expression of the combustion mass ratio is set separately for when ignition is performed at an ignition timing at which MBT is obtained and when ignition is performed at an ignition timing retarded from the ignition timing at which MBT is obtained. The described shaft torque calculation device. An operating state detecting means for detecting the operating state of the engine;
Estimating means for estimating one of the combustion mass at the crank angle during combustion of the gas in the combustion chamber, the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber, and the heat generation pattern based on the detected operating state When,
Shaft torque estimating means for estimating shaft torque generated by the engine based on any one of the estimated values;
An engine control device comprising: torque control means for performing torque control based on the estimated shaft torque. When the estimation means estimates the combustion mass ratio at the crank angle during combustion of the gas in the combustion chamber,
The shaft torque estimating means includes
An approximate expression calculating means for calculating an approximate expression of the combustion mass ratio at the crank angle during combustion;
Illustrated positive work estimated value calculating means for calculating the indicated positive work estimated value based on the approximate expression of the combustion mass ratio,
The engine control apparatus according to claim 7, further comprising: a shaft torque estimated value calculating unit that calculates a value obtained by adding a pumping loss and a friction loss to the illustrated positive work estimated value as a shaft torque estimated value.   The operation state of the engine is a heat generation start time, and the approximate expression of the combustion mass ratio is an exponential function expression using the heat generation start time and two constants. Engine control device.   The approximate expression of the combustion mass ratio is different between when ignition is performed at an ignition timing at which MBT is obtained and when ignition is performed at an ignition timing retarded from the ignition timing at which MBT is obtained. Engine control device.   The operating state of the engine is a combustion period from the start of combustion to a predetermined crank angle, and includes combustion period calculation means for dividing the combustion period into a plurality of calculations, and based on the plurality of combustion periods, the combustion mass ratio is calculated. The engine control apparatus according to claim 8, wherein an approximate expression is set.   12. The approximate expression of the combustion mass ratio is set separately for when ignition is performed at an ignition timing at which MBT is obtained and when ignition is performed at an ignition timing retarded from the ignition timing at which MBT is obtained. The engine control apparatus described in 1. The shaft torque estimated value calculation means includes
First shaft torque estimated value calculation means for calculating a shaft torque estimated value when ignition is performed at an ignition timing at which MBT is obtained as a first shaft torque estimated value;
A second shaft torque estimated value calculating means for calculating a shaft torque estimated value when ignition is performed at an ignition timing retarded from an ignition timing obtained by MBT as a second shaft torque estimated value;
9. The engine control according to claim 8, wherein the torque control means is an air amount increasing means for increasing an intake air amount so that the second shaft torque estimated value coincides with the required request shaft torque. apparatus.   The torque control means sets the ignition timing to the retard side by a predetermined ignition timing retard amount so that the first axis torque estimated value matches the required axis torque when the estimated value of the first axis torque exceeds the required axis torque. 14. The engine control device according to claim 13, further comprising ignition timing correcting means for correcting the engine timing.   The shaft torque estimated value calculating means calculates a shaft torque estimated value when the ignition timing is retarded by the predetermined ignition timing retard amount from the ignition timing at which MBT is obtained as a third shaft torque estimated value. Shaft torque estimated value calculating means, wherein the torque control means sets the ignition timing so that when the third axis torque estimated value exceeds the required axis torque, the third axis torque estimated value coincides with the required axis torque. 15. The engine control device according to claim 14, further comprising ignition timing correction means for correcting to the retard side.

Images