F2000A177 Seoul 2000 FISITA World Automotive Congress June 12-15, 2000, Seoul, Korea A Nonlinear Dynamic Model of SI Engines for Designing Controller Paljoo Yoon1) Seungbum Park2) Myoungho Sunwoo3)* 1)2) 3) Grdauate Student, Hanyang University, Seoul, Korea Associate Professor, Hanyang University, Seoul, Korea In this paper, a nonlinear dynamic engine model is introduced, which is developed to represent a SI engine over a wide range of operating conditions. The model includes intake manifold dynamics, fuel film dynamics, and engine rotational dynamics with transport delays inherent in the four stroke engine cycles, and can be used for designing engine controllers. The model is validated with engine-dynamometer experimental data. The accuracy of the model is evaluated by the comparison of the simulated and the measured data obtained from a 2.0L inline four cylinder engine over wide operating ranges. The test data are obtained from 42 operating conditions of the engine. The speed range is from 1500[rpm] to 4000[rpm], and the load range is from 0.4[bar] to WOT. The results show that the simulation data from the model and the measured data during the engine test are in good agreement. Furthermore, this nonlinear engine model is mathematically compact enough to run in real time, and can be used as an embedded model within a control algorithm or an observer when a powertrain controller is designed and developed. Keywords: Engine Control System, Intake Manifold Pressure Dynamics, Fuel Film Dynamics, Engine Rotational Dynamics previously were developed under the assumption that the engine is operated near the stoichiometric AFR only, and they did not consider the change of MBT with respect to AFR. INTRODUCTION Due to more stringent fuel economy and emission legislations, improvements in transient engine control performance have arisen as important issues. To reduce tail-pipe emissions using a three-way catalytic converter, the air-to-fuel ratio (AFR) should be controlled very precisely in both the steady state and transient engine operations. (Mean AFR variation : ±0.2%) (2). In order to regulate the AFR precisely, most of the conventional engine control units currently used adopt the feedforward control schemes parallel to oxygen sensor feedback. However, most of the current engine control algorithms are designed on the basis of steady state engine maps, and open-loop correction terms are added to compensate for transient engine dynamics. Therefore, this type of controller cannot compensate control gain from manufacturing tolerance or engine aging, and it does not provide robustness to external disturbances. Many studies have conducted to alleviate these problems regarding the development of the control-oriented dynamic engine model (5)(7)(10)(12) and model-based control since the 1980s. Moskwa and Hedrick (1992) developed the mean torque predictive engine model based on the research results of Dobner (1983). Both presented very nice model descriptions in the intake manifold and fuel delivery model. However, the torque function they used to estimate the maximum possible torque comes from general simulation results, and which may cause errors for the application of specific engines. Additionally, the model validation results they presented showed relatively large errors. Hendricks and Sorenson (1990) developed a very accurate mean value engine model based on engine power instead of engine torque. However, engine torque is more adequate in describing the dynamic behavior of an engine than engine power. Most of the controllers dedicated to vehicle motion control require not engine power, but engine torque information for the proper engine control action. Furthermore, most engine models published This study introduces a nonlinear dynamic engine model that is applicable to a wide operating range of SI engines. The dynamic model uses the test data of both steady and transient operations and nonlinear estimation techniques. As shown in Figure 1, the engine model consists of three input variables (throttle angle, fuel flow rate, and spark timing), one disturbance (load torque), and three state variables (intake manifold pressure, engine speed, and fuel mass in the fuel film). The time constants of mass airflow through the throttle body and torque production model are relatively smaller than three state variables, so they are expressed as algebraic relations. Figure 1. Block Diagram of 3-State Engine Model EXPERIMENTAL SETUP For the development of the engine model, 2.0L, an inline 4 cylinder DOHC engine and an eddy current type dynamometer are used. The ignition timing, fuel injection timing and duration are controlled by a PC-based Engine 1 * Corresponding author. e-mail: [email protected] The discharge (or flow) coefficient, C D is a function of the throttle angle (i.e., throttle geometry) and the pressure ratio across the throttle body. It is obtained from the ratio of the ideal flow rate to the actual flow rate from the steady state flow test data from various operating conditions. Control System (ECS) developed at Hanyang University (17). The ignition (dwell time and spark timing) and fuel injection (timing and duration) of each cylinder can be precisely controlled by the ECS with the use of a crankshaft encoder and a cam sensor. An eight channel analog-to-digital converter, which the ECS employs, is used for the purposes of data acquisition and closed-loop control. The analog signals from the engine-dynamometer system, such as intake manifold pressure, mass airflow rate, AFR, throttle angle, and brake torque are measured at each engine event. The AFR is measured by a wideband oxygen sensor. A hot film type airflow sensor is used for the measurement of mass airflow through the throttle body. Furthermore, IMEP is calculated with a PC-based combustion analysis system with a cylinder pressure sensor. INTAKE MANIFOLD PRESSURE DYNAMICS For the derivation of the manifold pressure state equation, the conservation of air mass in the intake manifold and the differential form of the ideal gas law are used. With simple mathematical manipulations, the manifold pressure state equation, Eq (5) is obtained. RTman Pman = m − m ap Vman at ENGINE MODELING The mass airflow rate out of the intake manifold, m ap is represented by a speed-density algorithm shown in Eq (6). MASS AIRFLOW RATE THROUGH THROTTLE BODY m ap = (1) man = ηvol Pman = f1(n) + f 2 (n) Pman MA is the maximum possible flow rate through the specific throttle, and TC is a normalized flow as a function of the cross-sectional area, and PRI is a normalized flow as a function of pressure ratio. These are as follows: ( P A α MA = amb th max RTamb ), TC = PRI = γ +1 γ −1 2 2 2(1−γ ) 2 1 − PR γ γ +1 γ − 1 f1(n) = f10 + f11n + f12n 2 + f13n3 + f14n 4 , Ath (α ) Ath (α max ) PRI = 1.0 where PR = Pman , Pamb if PR ≥ PCR (6) (7) where f 2 (n) = f 20 + f 21n + f 22 n 2 Figure 2 shows the normalized air charge with various engine operating conditions. It can be easily understood from the figure that the man is approximately proportional to the manifold pressure at a constant engine speed, and shows very small variations with the engine speed. 1 1 PRγ VD η P N 120 RTman vol man Apart from constants, the speed-density relationship is represented by (ηvol Pman )× N , and the value in the parenthesis, η vol Pman , is proportional to the actual air charge per stroke (9). This quantity, η vol Pman , can be called the normalized air charge, m an . The normalized air charge is obtained by the steady state engine test, and is approximated with a polynomial equation. The flow of air through a throttle valve can be considered as a special case of the isentropic flow of compressible fluid (7). The mass air flow rate through the throttle, m at , including the discharge coefficient, C D can be rewritten by Eq (1). m at = C D ⋅ MA ⋅ TC ⋅ PRI (5) (2) otherwise γ γ −1 2 PCR = γ +1 If the shape of the throttle bore is assumed to be circular and the projection area of the throttle plate about the vertical plane is ellipsoidal, the cross-sectional area of throttle opening, Ath (α ) is represented Eq (3). cosα D 2 π d d sin − 1 − + D2 − d 2 Ath (α ) = D 2 1 − 4 cosαo 2 D 2 d cosα o D 2 cosα + d sin − 1 + D 2 cos 2 α − d 2 cos 2 α o D cosα 2 cosα 2 cosα o (3) Figure 2 Normalized air charge with Pman and N Since the SI engine operation is based on the engine events, an air charge per stroke has more important meanings than the normalized air charge during the process of the development of the engine model. As shown in Eq (8.1), the air charge per stroke is obtained from the integration of the mass air flow into the intake port between intake events, but the value obtained during At the small throttle opening, there is an error between the ideal cross-sectional area and the actual area due to machining variation (6). A correction for this error is: α o = 0.91α o* − 2.59 (4) 2 the steady state test is widely used because of its mathematically simple form. m = ∫ t IVC m (t )dt ac t IVO ap − m ff (t ) = e (8.1) kg m ap 30m ap VD s = = mac = 4 RT N rev min 2 stroke man N min 60 s rev where m an (8.2) ( 120V man , V Dη vol N 120 RT man m Pman, ss = V Dη vol N at t − τ τ f X m fi (0) + τ f X m fi (t )1 − e f (12.1) = constant at step test, and at steady state Therefore, fuel flow rate into the cylinder is m fc (t ) = (1 − X )m fi (t ) + ) where (9.1) 1 τf − m ff (t ) = m fi (t ) − X ∆m fi (t ) e t τ f (12.2) ∆m fi (t ) = m fi (t ) − m fi (0) If the throttle angle is assumed to be constant and only step fuel perturbation is applied, the equivalence ratio in the cylinder is represented by Eq (12.3). where τ man corresponds to the time constant of the intake manifold pressure state equation represented by the first order differential equation, and Pman,ss is the manifold pressure determined from the throttle angle and engine speed at a steady state. τ man = f m ff (0) = τ f X m fi (0) An intake manifold pressure state equation, Eq (9.1), can be obtained by substituting Eq (6) into Eq (5). dPman V η N RTman = − D vol Pman + m at α , Pman 120V man dt Vman 1 =− − Pman, ss P τ man man m fi (t ) t τ φ c (t ) = (9.2) t − (A/ F) s m (t ) − X ∆m (t ) e τ f fi m ac (t ) fi (12.3) Figure 3 and Figure 4 show the AFR response to the perturbations of fuel supply rate according to Eq (12.3). FUEL FILM DYNAMICS The characteristics of the fuel injection process is very complex and is largely influenced by various factors, such as injection schemes, spray patterns, injection timing, port geometry, intake manifold wall temperature, fuel rail pressure, and so on. However it is actually not easy to implement such complicated fuel delivery models which consider all aspects mentioned above to the real time engine control applications. Thus, the first order lumped parameter model expressed in Eq (10) is widely used to represent the fuel delivery process (1)(7). m ff = − 1 Figure 3 Sensitivity of A/F response to X produced by fuel perturbations with τ f =0.6 m ff + Xm fi τf m fv = (1 − X ) m fi m fc = m fv + (10) 1 m τ f ff In Eq (10), m fi is the injected fuel flow rate from the injector, X is a fraction of fuel that forms a fuel puddle on the walls, m ff is mass of fuel in the fuel puddle, and the actual fuel flow rate, m fc , which enters the cylinder, is the sum of m fv and the evaporated fuel from the film with the time constant τ . Therefore, the actual AFR in the cylinder is given by m ac / m fc , and it can also be represented by an equivalence ratio ( φ c = 1 / λ ). f φ c (t ) = ( A / F ) s m fc (t ) m ac (t ) Figure 4 Sensitivity of A/F response to τ f produced by fuel perturbations with X = 0.6 (11) Figure 3 shows how the AFR response is influenced by the value of X , Figure 4 shows the influence of the value of τ f . The AFR at the initial transient period goes up rapidly due to the influence of fast fuel m fv , and is changed exponentially according to the evaporation behavior of the fuel film. As shown in these figures, steady state fuel film dynamics is independent of fuel delivery model parameters, X and τ f , and these parameters affect the In order to understand the dynamic characteristics of AFR with the variation of important fuel delivery model parameters ( X and τ f ), the transient response of AFR with step fuel perturbation is predicted. If Eq (10) is regarded as a kind of an initial value problem, the solution is given by 3 Figure 5 and Figure 6 represent the (MBT )λs and test engine, respectively. transient AFR only. Furthermore, the changes of AFR at the initial stage of the transient period are determined by X , and the AFR response of the successive transient period is affected by the time constant, τ f . Actually, the fuel film model parameters, X and τ f , depend on various engine operating conditions (engine speed, load, AFR and coolant temperature, and so on), and the engine speed and coolant temperature give the dominant effects on model parameters (14). Eq. (13) represents the typical relationship between engine operating conditions and the parameters of fuel film model. X = a 0 + b0 ( N ) + [a1 + b1 ( N )]T EC (13.1) 2 + c T 3 ]/ N τ f = [c 0 + c1TEC + c 2TEC 3 EC (13.2) In this study, the AFR under various engine operating conditions are measured from the step fuel perturbation test. The wideband oxygen sensor is used to measure the AFR at the exhaust manifold of the test engine. The Levenberg-Marquardt Method is employed to estimate the nonlinear fuel film model parameters (15). ∆MBT (λ ) of Figure 5 MBT at stoichiometry TORQUE PRODUCTION MODEL The combustion process in the cylinder produces engine torque, and the quantity of the torque produced is influenced by the AFR of the mixture in the cylinder, spark timing, and combustion efficiency. The transport delay with respect to each engine event gives dynamic characteristics to the torque production model. In this study, it is assumed that the dynamic torque of the engine is a function of engine speed, air charge per stroke, spark timing, and the AFR. Under this assumption, the predictive torque production model is derived from the steady state engine experiments. One of the most important things in the development of the torque production model is the identification of the optimal spark timing (MBT) at each of the operating conditions. Most engine models published previously(5)(7)(10) were developed under the assumption that the engine is operated near the stoichiometric AFR only, and they did not consider the change of the MBT with respect to the change of the AFR. The recent advanced engine technologies, such as lean burn, and gasoline direct injection engine, motivate the development of the engine model, which is applicable to a wide range of AFR. The MBT at the stoichiometric AFR is identified under various engine speed and load conditions. After that, the identification of the MBT at various AFR conditions is proceeded. The difference between the MBT at the stoichiometric AFR and the MBT at an arbitrary AFR at various engine operating conditions is denoted by ∆MBT (λ ) . Based upon this approach, the MBT at an arbitrary AFR can be represented by Eq (14). (MBT )λ s ( (MBT )λ = (MBT )λ where ) ( ) = f m ac , n = m0 + m1n + m 2 n 2 + m3 + m 4 n m ac s + ∆MBT (λ ) Figure 6 MBT variation with respect to λ Under the assumption of negligible cylinder-bycylinder engine torque variations, the indicated torque is obtained from the measurement of IMEP of cylinder #1 using a cylinder pressure sensor (3). If AFR and spark timing are fixed, the indicated torque of an engine is a function of engine speed and air charge. Therefore the indicated torque at the MBT and the stoichiometric AFR can be represented by Eq (15). Figure 7 shows the test data and model data generated from Eq (15). (Tind )MBT , λ s ( ) = t 0 + t1n + t 2 n 2 + t 3 + t 4 n m ac (15) (14.1) (14.2) Figure 7 Indicated torque at MBT and ∆MBT (λ ) = d 0 + d1λ + d 2 λ 2 λ =1 The indicated torque at an arbitrary engine condition is obtained from Eq (15) multiplied by the efficiencies of the spark timing and the AFR which reflect the changed operating conditions. (MBT )λ : MBT at Stoichiometric AFR s ( (MBT )λ : MBT at arbitrary AFR ) ( ) ⋅ η ∆SA ⋅ η λ Tind n, m ac , η ∆SA , η λ = Tind MBT , λ s 4 (16) where η ∆SA = s 0 + s1 (∆SA) + s 2 (∆SA)2 production model, three different transport delays are introduced. This approach is reasonable for real time control applications because the combustion process of the engine occurs much faster than the transport dynamics of air and fuel. However, these transport delays are subjected to provide stability problems for feedback systems, which must be considered in the control synthesis (10). ∆SA = (MBT )λ − SA η λ = l 0 + l1λ + l 2 λ 2 The η∆SA , shown in Figure 8, denotes the spark timing efficiency, which is a function of the difference of spark timing between the MBT at an arbitrary AFR condition and spark timing at the current operation. The value of η ∆SA at the MBT is 1. The η λ , depicted in Figure 9, represents the influence of the AFR to the indicated torque, and is normalized by 1 at a stoichiometric AFR. Figure 10 Schematic diagram of torque production model Figure 8 Spark timing efficiency TRANSPORT DELAY BETWEEN ENGINE EVENTS In the instant of closing the intake valve, the amount and AFR of the mixture in the cylinder that are available for torque production are determined. Therefore, it seems natural to synchronize the sample times with the intake events. If each event of the engine is concentrated at one instant of time (5), the useful transport delays used in this engine model are: Intake to exhaust delay, Figure 9 A/F efficiency Intake to torque delay, Actually, the indicated torque generated by the combustion of an in-cylinder mixture is decreased by both engine friction and pumping loss. The difference between the indicated torque and the friction and pumping loss becomes the brake torque. The friction and pumping torque, T , is estimated from the difference between the indicated torque measured from the steady state cylinder pressure and the brake torque obtained from an engine dynamometer. It is well known that the friction and pumping loss is a function of the engine speed and intake manifold pressure (16). As a result of that, the friction and pumping torque can be approximated by a polynomial of engine speed and air charge. ( ) ∆ϕ IT1 Injection to torque delay, Spark to torque delay, = 489° = 279° ∆ϕ IT 2 ∆ϕ ST (18) = 585° = 75° ROTATIONAL DYNAMICS OF ENGINE f /p T f / p = p 0 + p1n + p 2 n 2 + p 3 + p 4 n m ac ∆ϕ IE The rotational dynamics of an engine is modeled under the assumption of a lumped parameter system with constant inertia. Using Newton's Second Law, the state equation for engine speed is given by: J eff (17) where Summarizing the torque production model proposed in this study, this model has four inputs (engine speed, air charge, fuel flow rate, and spark timing), six submodels approximated by polynomials, and brake torque as an output. The overall structure of the torque production model developed in this study is shown in Figure 10. In order to provide the dynamic characteristics of the torque with dN (t ) 30 (t ) − T / (t ) − T (t ) = T f p L π ind dt ( ) (t − ∆t IT1 ) ⋅ η ∆SA (t − ∆t SA ) ⋅ η λ (t − ∆t IT 2 ) Tind (t ) = Tind MBT , λ s ∆t IT1 = ∆ϕ IT1 6 N (t ) , ∆t IT 2 = ∆ϕ IT 2 6 N (t ) , ∆t ST = ENGINE MODEL VALIDATION 5 (19) ∆ϕ ST 6 N (t ) To validate the developed engine model, steady and transient tests in a wide operating range must be carried out. The algebraic equation derived in this study is an approximation of experimental data using the LevenbergMarquardt Method (18), which is a kind of nonlinear least square estimation technique. The test data and approximation results are depicted in each figure in a comparative manner. The sampling interval for the simulation is selected to be synchronized with the intake event because of the event-based operation characteristics of the engine. ts = 120 N cyl N Case I : Constant Engine Speed To prove the effectiveness of the torque production model proposed in this study, the dynamometer is operated at constant speed mode. Figure 11 shows the comparison between the model and the measurement when step throttle input is applied. The model outputs of the manifold pressure and the brake torque show in good agreement with experimental data. However, for the brake torque, a slight time delay is observed during the transient period between the predicted and the measured data. The time delay may be caused by the digital-to-analog conversion and the signal output process in the dynamometer controller. (20) In order to evaluate the accuracy of the model, test data from various operating conditions are collected, and a comparative study is done with the predicted values of various variables of the model using the same input profiles (throttle angle, fuel flow rate, and spark timing). MODEL VALIDATION AT STEADY STATE In order to confirm the accuracy of the model, the mean values of the errors in the predictions of the model and the test data are calculated from 42 operating conditions of the engine. The speed range is from 1500[rpm] to 4000[rpm], and the load range is from 0.4[bar] to WOT. The results given in Table 1 show that the relatively large error (about 10%) is observed at low load conditions but the slight error is observed in the remaining operating regions. The larger error observed in the low load region is due to the characteristics of the eddy current dynamometer, because it is well known that the eddy current type dynamometer causes poor torque control performance at low load conditions. Furthermore, the shaft that connects the engine and dynamometer is not rigid enough, and thus the influence of the shaft compliance on the test data is much larger at low load conditions, further studies are needed to accommodate the effects of the dynamometer and shaft compliance on the model performance. Table 1. Mean errors in the steady state predictions of the engine model Variable % error N Pman m at -5.01 Figure 11 Validation of engine model at constant engine speed (N=2000rpm) 5.68 Case II : Constant Load Torque -3.41 In this case, the dynamometer is set at a constant torque operation, and the dynamic behaviors of three state variables are examined. Figure 12 shows the comparison between the model and the measurement when step throttle transient is applied. The predicted manifold pressure, engine speed, and the AFR at the exhaust pipe show in good agreement with the experimental data. MODEL VALIDATION AT TRANSIENT STATE For validation of the engine model during transients, the dynamometer is operated at a constant speed and a constant torque mode, respectively. A set of time responses was recorded using throttle transients. Among three state variables, the engine speed and the manifold pressure are measured directly, and the AFR at the exhaust pipe is measured instead of fuel mass in the film because the fuel film dynamics is not easy to measure. The torque production model is verified with the brake torque from the output of the dynamometer controller. 6 (4) The engine model is validated to an accuracy of ±5% by the comparison of the predicted and the measured data obtained from wide operating ranges. (5) This model is mathematically compact enough to run in real time, and can be used as an embedded model within a control algorithm or an observer. (6) For the improvement of the reliability of this model in low load regions, the effect of the eddy current dynamometer's characteristics and shaft compliance to the model accuracy should be considered. REFERENCES [1] Aquino, C.F. 1981. Transient A/F Control Characteristics of the 5 Liter Central Fuel Injection Engine. SAE Paper 810494 [2] Benninger, N.F. and Plapp, G. 1991. Requirements and Performance of Engine Management Systems under Transient Conditions. SAE Paper 910083. [3] Brunt, M.F. and Emtage, A.L. 1996. Evaluation of IMEP Routines and Analysis Errors. SAE Paper 960609 [4] Chang, C.-F., Fekete, N.P. and Powell, J.D. 1993. Engine Air-Fuel Ratio Control Using an Event-Based Observer. SAE Paper 930766 [5] Dobner, D.J. 1983. Dynamic Engine Models for Control Development - Part I : Nonlinear and Linear Model Formulation. International Journal of Vehicle Design, Technological Advances in Vehicle Design Series [6] Harrington, D.L. and Bolt, J.A. 1970. Analysis and Digital Simulation of Carburetor Metering. SAE Paper 700082 [7] Hendricks, E. and Sorenson, S.C. 1990. Mean Value Modeling of Spark Ignition Engines. SAE Paper 900616 [8] Hendricks, E. and Sorenson, S.C. 1991. SI Engine Controls and Mean Value Engine Modeling. SAE Paper 910258 [9] Hendricks, E., Chevalier, A, Jensen, A. and Sorenson, S.C. 1996. Modeling of the Intake Manifold Filling Dynamics. SAE Paper 960037 [10] Moskwa, J.J. and Hedrick, J.K. 1992. Modeling and Validation of Automotive Engines for Control Algorithm Development. Trans. ASME, J. Dynamic Systems, Measurement, and Control, Vol. 114 [11] Nasu, M., Obata, A. and Abe, S. 1996. Model-based Fuel Injection System for SI Engines. SAE Paper 961188 [12] Powell, B.K. and Cook, J.A. 1987. Nonlinear Low Frequency Phenomenological Engine Modeling and Analysis. Proc. American Control Conference [13] Powell, J.D., Fekete, N.P. and Chang, C.-F. 1998. Observer-based Air-Fuel Ratio Control. IEEE Control Systems, Vol. 18, No. 5 [14] Shayler, P.J., Teo, Y.C. and Scarisbrick, A. 1995. Fuel Transport Characteristics of SI engines for Transient Fuel Compensation. SAE Paper 950067 [15] Sunwoo, M., Yoon, P. and Ju, J. 1998. A Study on the Design of Fuel Film Compensator and A/F Controller in SI Figure 12 Validation of engine model at constant load torque CONCLUSION A nonlinear dynamic engine model which is applicable to a wide operating range of SI engines is developed by using the steady and transient test data and a nonlinear estimation technique. The following conclusions are derived in this study. (1) The engine model is comprised of three input variables (throttle angle, fuel flow rate, and spark timing) including one disturbance (load torque). The model is described with three state variables (intake manifold pressure, engine speed, and fuel mass in the fuel film), and major subsystems (mass airflow through throttle, torque production model, etc.) of the engine expressed in algebraic equations are included. (2) The applicable operating range of the engine model is expanded by the inclusion of the changes of the MBT with respect to AFRs, and it also can be extended to the lean burn operations. (3) The influence of fuel delivery model parameters to transient AFR characteristics is examined through the analysis of the sensitivity of AFR response to model parameters. The AFR change at the initial stage during the transient period is determined by X , and the response of the successive transient period is mainly affected by time constant, τ f . 7 Engines. Autumn Conference Proceeding of KSAE, pp. 154~159 (in Korean) [16] Taylor, C.F. 1985. The Internal Combustion Engine in Theory and Practice. Vol. 1, 2nd ed. MIT Press [17] Yoon, P., Kim, M. and Sunwoo, M. 1998. A Study on Design and Development of an Engine Control System based on Crank Angle. Transactions of KSAE, Vol. 6, No .4 (in Korean) [18] Zielinski, T.J. and Allendoerfer, R.D. (1997) Least Squares Fitting of Nonlinear Data in the Undergraduate Laboratory. J. of Chemical Education, Vol. 74, No. 8 8
© Copyright 2024