【溫馨提示】 dwg后綴的文件為CAD圖，可編輯，無水印，高清圖，壓縮包內文檔可直接點開預覽，需要原稿請自助充值下載，請見壓縮包內的文件及預覽，所見才能所得，請細心查看有疑問可以咨詢QQ：414951605或1304139763
Journal of Engineering Mathematics 44: 57–82, 2002. ? 2002 Kluwer Academic Publishers. Printed in the Netherlands. Modelling gas motion in a rapid-compression machine M.G. MEERE 1 , B. GLEESON 1 and J.M. SIMMIE 2 Department of Mathematical Physics, NUI, Galway, Ireland 2 Department of Chemistry, NUI, Galway, Ireland Received 25 July 2001; accepted in revised form 8 May 2002 Abstract. In this paper, a model which describes the behaviour of the pressure, density and temperature of a gas mixture in a rapid compression machine is developed and analyzed. The model consists of a coupled system of nonlinear partial differential equations, and both formal asymptotic and numerical solutions are presented. Using asymptotic techniques, a simple discrete algorithm which tracks the time evolution of the pressure, temperature and density of the gas in the chamber core is derived. The results which this algorithm predict are in good agreement with experimental data. Key words: gasdynamics, rapid-compression machines, shock-waves, singular perturbation theory 1. Introduction 1.1. RAPID-COMPRESSION MACHINES A rapid-compression machine is a device used to study the auto-ignition of gas mixtures at high pressures and temperatures, with particular reference to auto-ignition in internal combus- tion engines; see [1–3]. A typical combustion engine is a very dirty and complex environment, and this has prompted the development of rapid-compression machines which enable the scientific study of compression and ignition in engines in a cleaner and simpler setting. In Figure 1 we schematically illustrate a two-piston rapid-compression machine, such as the one in the department of Chemistry at NUI, Galway. However, single-piston machines, with a piston at one end and a stationary solid wall at the other, are more typical. The analysis developed in this paper is appropriate to both single- and two-piston machines. The operation of a rapid-compression machine is very simple - the pistons are simul- taneously driven in pneumatically, compressing the enclosed gas mixture, thereby causing the gas pressure, temperature and density to rise quickly. In Figures 1(a), 1(b) and 1(c) we schematically represent a rapid-compression machine prior to, during, and after compression, respectively. The ratio of the final volume to the initial volume of the compression chamber for the machine at NUI, Galway is about 1:12, this value being typical of other machines. At the end of the compression, the gas mixture will typically have been pushed into a pressure and temperature regime where auto-ignition can occur. In Figure 2, we display an experimental pressure profile for a H 2 /O 2 /N 2 /Ar mixture which has been taken from Brett et al. [4], with the kind permission of the authors. In this graph, the time t = 0 corresponds to the end of the compression time. We note that, for the greater part of the compression, the pressure in the chamber is rising quite gently, but that towards the end of the compression (that is, just before t = 0), there is a steep rise in the pressure. After compression, the pressure profile levels off as expected; the extremely steep rise at the end of 58 M.G. Meere et al. Figure 1. Schematic illustrating the operation of a rapid-compression machine; we have shown the configuration (a) prior to compression, (b) during compression and (c) after compression. Figure 2. An experimental pressure profile for a gas mixture H 2 /O 2 /N 2 /Ar = 2/1/2/3, as measured in the rapid-compression machine at NUI, Galway. It is taken from [4], and has an initial pressure of 0·05 MPa and an initial temperature of 344 K. Modelling gas motion in a rapid-compression machine 59 the profile corresponds to the ignition of the mixture. We note that the compression time and the time delay to ignition after compression are both O(10) ms. The pressure history is the only quantity which is measured in experiments. However, the temperature in the core after compression is the quantity which is of primary interest to chemists since reaction rates depend mainly on temperature for almost all systems, although there may also be some weaker pressure dependence. Measuring temperature accurately in the core can be problematic because of the presence of a thermal boundary layer; see the comments below on roll-up vortices. However, with the experimental pressure data in hand, the corresponding temperatures can be estimated using the isentropic relation ln(p/p i ) = integraldisplay T T i γ(s) s(γ(s)? 1) ds, (1) where (T i ,p i ) are the initial values for the core temperature and pressure, (T,p) are these quantities at some later time, and γ(s)is the specific heat ratio at temperature s. In exper- iments, the initial core temperature is typically O(300 K), while the core temperature after compression is usually O(1000 K). In this paper, we shall consider only the behaviour of the gas mixture during compression; the post-compression behaviour is not considered here, but this will form the subject for future work. Nevertheless, the model presented here does provide a reasonable description of the post-compression behaviour of a single species pure gas, or an inert gas mixture; see Section 3.5. 1.2. THE MODEL We suppose that the compression chamber is located along 0 x 0. This assumption is actually quite a strong one in this context since higher-dimensional effects are frequently observed in experiments, roll-up vortices near the corner regions defined by the piston heads and the chamber wall being particularly noteworthy; see, for example, [5]. These vortices arise due to the scraping by the pistons of the thermal boundary layer at the chamber wall, and they can, and frequently do, disturb the gas motion in the core of the compression chamber. However, the justification for the one-dimensional model studied here is twofold: (i) the corner vortices can be successfully suppressed by introducing crevices at the piston heads which swallow the thermal boundary layer as the pistons move in (see Lee [6]), rendering the gas motion away from the chamber walls one-dimensional to a good approximation, and, (ii) the study of the one-dimensional model provides a useful preliminary to the study of higher-dimensional models. We now give the governing equations for our one-dimensional model. A reasonably com- plete derivation of the governing equations for a multi-component reacting gas can be found in the appendices of [7]; these standard derivations are not reproduced here. The model which we shall study includes a number of simplifying assumptions and these will be clearly stated as they arise. The equation expressing conservation of mass is given by ?ρ ?t + ? ?x (ρv) = 0, 60 M.G. Meere et al. where ρ = ρ(x,t)and v = v(x,t) are the density and the velocity of the gas, respectively, at location x and time t. It should be emphasized that these quantities refer to a gas mixture, so that if there are N different species in the mixture then ρ = N summationdisplay i=1 ρ i , where ρ i = ρ i (x,t) is the density of species i. Also, the velocity v above refers to the mass- averaged velocity of the mixture, that is, v = N summationdisplay i=1 Y i v i , where Y i = ρ i /ρ and v i = v i (x,t) are the mass fraction and velocity, respectively, of species i; see [7]. Neglecting body forces and viscous effects, the equation expressing conservation of mo- mentum is given by ?v ?t +v ?v ?x =? 1 ρ ?p ?x , where p = p(x,t)is the pressure. We assume that the gas mixture is ideal, so that the equation of state is given by p = R M ρT, (2) where T = T(x,t)is the temperature, R is the universal gas constant (8·314 JK ?1 mol ?1 ), and M is the average molecular mass of the mixture. This last quantity is given by M = N summationdisplay i=1 n i W i (mA), where n i and W i give the number fraction and molecular weight, respectively, of species i, m is the atomic mass unit (1·661 × 10 ?27 kg) and A is Avogadro’s number (6·022 × 10 23 molecules mol ?1 ). The equation expressing conservation of energy is given by (see [7] or [8]) ρ parenleftbigg ?u ?t +v ?u ?x parenrightbigg =?M parenleftbigg ?q ?x +p ?v ?x parenrightbigg , (3) where u = u(x,t) is the internal energy of the gas mixture and q = q(x,t) is the heat flux. We also have the thermodynamic identity u= N summationdisplay i=1 h i Y i ?Mp/ρ, (4) where the enthalpies h i = h i (T) are given by h i (T) = h i (T 0 )+ integraldisplay T T 0 c p,i (s)ds, i = 1,2,.,N, (5) Modelling gas motion in a rapid-compression machine 61 where T 0 is some reference temperature and the c p,i (T) are the specific heats at constant pressure for the N species. When diffusion velocities and the radiant heat (again, see [7] for more details) are neglected, the expression for the heat flux is given by q =?λ(T) ?T ?x , (6) where λ(T) is the thermal conductivity. The mass fractions Y i = ρ i /ρ are not necessarily constant since chemical reactions can change the composition of the mixture. However, for many systems such chemical effects can be neglected in the analysis of the compression because the gas mixture is ‘cold’ for most of the compression time. The core temperature will only rise to a level where chemical reactions can have a significant effect near the end of the compression, and the duration of this period is typically very short (a couple of milliseconds usually). Nevertheless, it is possible for some chemical reactions to proceed sufficiently rapidly for them to significantly influence the compression behaviour. However, we do not attempt to model systems which exhibit this behaviour here and take the Y i to be constant during compression. Substituting (4) and (6) in (3), and using (5), we have the final form of the equation expressing conservation of energy: ?T ?t +v ?T ?x = M ρ(c p (T)?R) parenleftbigg ? ?x parenleftbigg λ(T) ?T ?x parenrightbigg ?p ?v ?x parenrightbigg , where c p = N summationdisplay i=1 Y i c p,i is the mass-averaged specific heat. 1.3. BOUNDARY AND INITIAL CONDITIONS We suppose that the left and right pistons move with constant velocities V 0 and ?V 0 ,re- spectively, so that their motions are given by x = V 0 t and x = 2L ? V 0 t. In reality, the pistons in a rapid-compression machine will spend some of the compression time accelerating from rest and decelerating to rest, and this is not difficult to incorporate into the analysis given below. However, rather than complicate the analysis unnecessarily at the outset by considering variable piston velocity, we shall simply quote the results for general piston motion once the constant velocity case has been completed; see Section 3.4. Throughout the compression, we assume that the temperature of the walls of the chamber remain at their initial constant value, which we denote by T 0 . Hence, at the left piston, we impose v = V 0 ,T= T 0 on x = V 0 t, while at the right piston we set v =?V 0 ,T= T 0 on x = 2L?V 0 t. The gas in the chamber is initially at rest and we suppose that v = 0,T= T 0 ,p= p 0 ,ρ= ρ 0 at t = 0, 62 M.G. Meere et al. where p 0 and ρ 0 are constants. Clearly, in view of (2), we have p 0 = R M ρ 0 T 0 . However, the above are not quite the boundary and initial conditions that are considered in this paper. For the conditions described above above we have the symmetry v(x,t)=?v(2L?x,t), T(x,t)= T(2L?x,t), p(x,t)= p(2L?x,t), ρ(x,t)= ρ(2L?x,t). We exploit this behaviour by halving the spatial domain, considering the gas motion in V 0 t xq + ,wehavev = 0, p = ρ = T = 1. For x 0, but this amounts to nothing more than requiring that the pistons travel at velocities which do not exceed the speed of sound in the gas. Recall that the maximum speed of the pistons is O(10 ms ?1 ), while the speed of sound in gases under typical conditions is frequently O(300 ms ?1 ). Substituting (17) in (16) 3 , and integrating subject to the conditions v ? 0 → 0,p ? 0 → 1as z ? →+∞,wehave Modelling gas motion in a rapid-compression machine 71 p ? 0 = 1 +˙q + 0 v ? 0 /θ. (19) Letting z ? →?∞in (19) we obtain P s (t) = 1 +˙q + 0 /θ, (20) so that, T c (t) = (1 +˙q + 0 /θ)(1 ? 1/˙q + 0 ), (21) both of which are constant since, as we shall now show, ˙q + 0 is constant. The prediction that (p,ρ,T)are constant to leading order in the outer region behind the wave-front is clearly consistent with the numerical solution displayed in Figure 3. Substituting (17) and (19) in (16) 4 and integrating subject tov ? 0 → 0,T ? 0 → 1,?T ? 0 /?z ? → 0asz ? →+∞,weget ?˙q + 0 μ(T ? 0 ) = λ(T ? 0 ) ?T ? 0 ?z ? ?(v ? 0 +˙q + 0 v ?2 0 /2θ), (22) where μ(T ? 0 ) = integraldisplay T ? 0 1 ds γ(s)? 1 . Letting z ? →?∞in (22), we obtain ˙q + 0 μ parenleftbig (1 +˙q + 0 /θ)(1 ? 1/˙q + 0 ) parenrightbig = 1 +˙q + 0 /2θ, (23) which determines ˙q + 0 , completing the specification of the leading order outer problem. It is clear that the solution for ˙q + 0 to (23) does not depend on t,sothatq + 0 has the form αt where α is a constant. Using (9), we have μ(T ? 0 ) = 1 γ 1 ln parenleftbigg γ 0 +γ 1 T ? 0 ? 1 γ 0 +γ 1 ? 1 parenrightbigg , so that (23) becomes ˙q + 0 γ 1 ln parenleftbigg γ 0 ? 1 +γ 1 (1 +˙q + 0 /θ)(1 ? 1/˙q + 0 ) γ 0 +γ 1 ? 1 parenrightbigg = 1 +˙q + 0 /2θ, (24) which is an equation that is easily solved numerically for ˙q + 0 for given values of γ 0 , γ 1 and θ. In the limit γ 1 → 0 (so that γ(T)≡ γ 0 ) this expression reduces to a quadratic in ˙q + 0 which can be solved to give ˙q + 0 = 1 4 parenleftBig γ 0 + 1 ± radicalbig (γ 0 + 1) 2 + 16γ 0 θ parenrightBig , with the positive solution being clearly the relevant one here. The numerical solution of (24) for γ 1 negationslash= 0 is usually unnecessary. Recalling that θ = O(10 3 ) typically, and considering the behaviour of (24) for θ greatermuch 1, we can easily show that ˙q + 0 ～ radicalbig (γ 0 +γ 1 )θ for θ greatermuch 1. (25) In dimensional terms, this expression for ˙q + 0 is 72 M.G. Meere et al. radicalBigg γ(T 0 )p 0 ρ 0 , which is the familiar expression for the speed of sound in an ideal gas. We favour the simpler expression (25) over (24) for the algorithm described in Section 3.4. It is worth noting here that the relations (18), (20) and (23) could also have been obtained using integral forms for the conservation laws, and it is not necessary (although it is preferable) to consider the detail of the transition layer. For example, conservation of mass implies that d dt parenleftbiggintegraldisplay 1 t ρ(x,t)dx parenrightbigg = 0, which at leading order gives d dt parenleftBigg integraldisplay q + 0 t ρ ? 0 (x,t)dx + integraldisplay 1 q + 0 1dx parenrightBigg = 0, and this leads to (18). 3.1.4. Summary The motion of the wave-front, x = q + (t;ε), is such that as ε → 0, q + (t;ε) ～ q + 0 (t),where q + 0 (t) is determined by solving (24) subject to q + 0 (0) = 0. For xq + , p = ρ = T = 1andv = 0. 3.2. THE FIRST REFLECTION OF THE WAVE FROM THE CENTRE-LINE When the wave-front reaches the centre-line, it reflects off the identical opposing wave, and then moves from right to left towards the incoming piston. The leading-order behaviour ahead of the wave is now known from the calculations of the previous subsection. A numerical solution illustrating this case is given in Figure 4. The leading-order behaviour in the boundary layer near the piston is clearly unchanged from that considered in Section 3.1.1 and requires no further discussion. 3.2.1. Outer region We denote the motion of the reflected wave by x = q ? (t;ε).Forxq ? ,wehavev = o(1) and we pose p ～ p + 0 (x,t), ρ ～ ρ + 0 (x,t), T ～ T + 0 (x,t) to obtain p + 0 = ρ + 0 T + 0 , ?ρ + 0 ?t = 0, ?p + 0 ?x = 0, ?T + 0 ?t = 0, so that p + 0 = P sr ,ρ + 0 = g(x), T + 0 = P sr /g(x), Modelling gas motion in a rapid-compression machine 73 where P sr , which is constant, and g(x) are determined below by matching. 3.2.2. Transition region This is located at z ? = O(1) where x = q ? (t;ε)+εz ? . It gives the location of the narrow region over which v drops from v ～ 1tov = o(1); the transition region is also clearly identifiable in the solutions for p, ρ and T; see Figure 4. In z ? = O(1) we pose q ? ～ q ? 0 (t), p ～ p ? 0 (z ? ,t),ρ～ ρ ? 0 (z ? ,t),v～ v ? 0 (z ? ,t),T ～ T ? 0 (z ? ,t), to obtain leading-order equations which have precisely the same form as (16). Integrating and matching in a manner similar to that described in Section 3.1.3, we obtain ρ ? 0 = ˙q + 0 (˙q ? 0 ? 1) (˙q + 0 ? 1)(˙q ? 0 ?v ? 0 ) , v ? 0 = 1 + θ(˙q + 0 ? 1) ˙q + 0 (˙q ? 0 ? 1) (p ? 0 ? 1 ?˙q + 0 /θ), ˙q + 0 (˙q ? 0 ? 1) ˙q + 0 ? 1 (μ(T c )?μ(T ? 0 )) = λ(T ? 0 ) ?T ? 0 ?z ? ? parenleftbigg (1 +˙q + 0 /θ)(v ? 0 ? 1)+ ˙q + 0 (˙q ? 0 ? 1) 2θ(˙q + 0 ? 1) (v ?2 0 ? 1) parenrightbigg . (27) Letting z ? →+∞in (27) gives g(x)≡ ˙q + 0 (˙q ? 0 ? 1) ˙q ? 0 (˙q + 0 ? 1) ,P sr = 1 + ˙q + 0 (˙q + 0 ?˙q ? 0 ) θ(˙q + 0 ? 1) , T + 0 = ˙q ? 0 ˙q + 0 (˙q ? 0 ? 1) parenleftbigg ˙q + 0 ? 1 + ˙q + 0 θ (˙q + 0 ?˙q ? 0 ) parenrightbigg , (28) where the constant reflected wave speed ˙q ? 0 is determined as the negative solution to ˙q + 0 (˙q ? 0 ? 1) γ 1 (˙q + 0 ? 1) ln parenleftbigg γ 0 ? 1 +γ 1 T c γ 0 ? 1 +γ 1 T + 0 parenrightbigg = 1 + ˙q + 0 θ ? ˙q + 0 (˙q ? 0 ? 1) 2θ(˙q + 0 ? 1) , (29) where T c is given by (21) and T + 0 is given in (28). Considering the behaviour of this last expression for θ greatermuch 1, we find that ˙q ? 0 ～± √ (γ 0 +γ 1 )θ, the negative solution being the relevant one now. It is this simpler form which we shall use for the algorithm described in Section 3.4. 3.2.3. Summary The location of the reflected wave-front, x = q ? (t;ε), is such that as ε → 0, q ? (t;ε) ～ q ? 0 (t),where˙q ? 0 (t) is determined as the solution to (29). For xq ? ,wehavev = o(1) and p ～ 1 + ˙q + 0 (˙q + 0 ?˙q ? 0 ) θ(˙q + 0 ? 1) ,ρ～ ˙q + 0 (˙q ? 0 ? 1) ˙q ? 0 (˙q + 0 ? 1) , T ～ ˙q ? 0 ˙q + 0 (˙q ? 0 ? 1) parenleftbigg ˙q + 0 ? 1 + ˙q + 0 θ (˙q + 0 ?˙q ? 0 ) parenrightbigg . 74 M.G. Meere et al. 3.3. THE WAVE TRAVELS OVER AND BACK IN THE CHAMBER FOR THE N th TIME Most of the notation required here has previously been introduced in Section 2.3. 3.3.1. The wave travels down the chamber for the N th time Denoting the location of the wave-front by q + N ,wehaveforxq + N ,wehavev = o(1) and p ～ p 2N?2 ,ρ～ ρ 2N?2 ,T～ T 2N?2 , where (p 2N?2 ,ρ 2N?2 ,T 2N?2 ) are constants. If we now consider the transition layer of thick- ness O(ε) about q + N , and perform calculations which are almost identical to those of Sec- tion 3.1.3, we find that ρ 2N?1 = ˙q + N0 ˙q + N0 ? 1 ρ 2N?2 , p 2N?1 = p 2N?2 + ˙q + N0 θ ρ 2N?2 , T 2N?1 = p 2N?1 ρ 2N?1 , ˙q + N0 ρ 2N?2 (μ(T 2N?1 )?μ(T 2N?2 )) = p 2N?2 + ˙q + N0 2θ ρ 2N?2 . (30) Using (9), and considering the behaviour of (30) 4 for θ greatermuch 1, we can easily show that ˙q + N0 ～ radicalbig γ(T 2N?2 )T 2N?2 θ = radicalBigg γ(T 2N?2 )p 2N?2 θ ρ 2N?2 for θ greatermuch 1. We note from this expression that as the pistons compress the gas in the chamber core, the speed of the wave increases in proportion to the square root of the rising temperature (for γ(T)constant). If γ(T)≡ γ 0 , then we can solve exactly for ˙q + N0 to obtain ˙q + N0 = 1 4 parenleftBig γ 0 + 1 + radicalbig (γ 0 + 1) 2 + 16γ 0 θT 2N?2 parenrightBig . (31) 3.3.2. The wave reflects off the centre-line for the N th time We denote the location of the wave-front by x = q ? N .Asε → 0, we have for xq ? N ,wehavev = o(1) and p ～ p 2N ,ρ～ ρ 2N ,T～ T 2N , with (p 2N ,ρ 2N ,T 2N ) constant. Considering the transition layer about q ? N , it is readily shown that Modelling gas motion in a rapid-compression machine 75 ρ 2N = ˙q ? N0 ? 1 ˙q ? N0 ρ 2N?1 , p 2N = p 2N?1 ? ˙q ? N0 θ ρ 2N , T 2N = p 2N ρ 2N , ˙q ? N0 ρ 2N (μ(T 2N?1 )?μ(T 2N )) = p 2N + ˙q ? N0 2θ ρ 2N . (32) Considering the behaviour of (34) 4 for θ greatermuch 1, we have ˙q ? N0 ～? radicalbig γ(T 2N?1 )T 2N?1 θ =? radicalBigg γ(T 2N?1 )p 2N?1 θ ρ 2N?1 for θ greatermuch 1. For γ(T)≡ γ 0 , we have the exact expression ˙q ? N0 = 1 4 parenleftBig 3 ?γ 0 ? radicalbig (3 ?γ 0 ) 2 + 16γ 0 θT 2N?1 + 8(γ 0 ? 1) parenrightBig . (33) 3.4. THE ALGORITHM 3.4.1.