54 Journal of the Korean Society of Propulsion Engineers Vol. 18, No. 5, pp. 54-61, 2014 Research Paper DOI: http://dx.doi.org/10.6108/kspe.2014.18.5.054 DSMC 방법을이용한로켓플룸의해석 전우진 a 백승욱 a, * 박재현 b 하동성 c Rocket Plume Analysis with DSMC Method Woojin Jeon a Seungwook Baek a, * Jaehyun Park b Dongsung Ha c a Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, Korea b Department of Aerospace and System Engineering, ReCAPT, Gyeongsang National University, Korea c Advanced Propulsion Technology Center, Agency for Defense Development, Korea * Corresponding author. E-mail: swbaek@kaist.ac.kr ABSTRACT In this study, a plume exhausted from rocket nozzle is investigated by using an unstructured 2-dimensional axisymmetirc DSMC code at various altitude. The small back-pressure to total-pressure ratio(p b /P o ) and large P b /P o represent low and high altitude condition, respectively. At low altitude, the plume shows a typical complicated structure (e.g. Mach disk) of underexpanded jet while the high altitude plume experiences plain expansion. The various features of exhaust plume is discussed including density, translational/rotational temperature, Mach number and Knudsen number. The results shows that even at 20 km altitude where the freestream Knudsen number is small as 1.5 10-5, the transitional and rarefied flow regimes can occur locally within the plume. It confirms the necessity of DSMC computation at low altitude. 초 록 본연구에서는비정렬격자계를사용하는 2차원축대칭 DSMC 법을사용하여로켓노즐에서사출되는플룸을해석하였다. 오리피스의출구전압에대한배압의비율이높은경우와낮은경우의플룸에대하여해석을실시하여저고도와고고도를대표하는두가지조건에서플룸유동의차이를관찰하였다. 저고도플룸은 Mach disc 등복잡한유동구조를보인반면고고도플룸은단순팽창만을보였으며, 유동이상류방향으로심하게꺾였다. 또한고도 20 km의대기조건에서소형로켓노즐에서사출되는플룸에대한해석을수행하여연속체해석결과와비교하였으며과소팽창되는로켓플룸의유동구조가잘나타났다. 또한, 플룸내부에국지적인천이유동이발생할수있음을확인하였다. Key Words: Rarefied Flow( 희박유동 ), Direct Simulation Monte Carlo( 직접모사법 ), High Altitude ( 고고도 ), Rocket Plume( 로켓플룸 ) Received 11 December 2013 / Revised 28 August 2014 / Accepted 5 July 2014 Copyright C The Korean Society of Propulsion Engineers pissn 1226-6027 / eissn 2288-4548 [ 이논문은한국추진공학회 2013 년도추계학술대회 (2013. 12. 4-5, 경주현대호텔 ) 발표논문을심사하여수정 보완한것임.] : relative speed Nomenclature This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License(http://creativecommons.org /licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
제 18 권제 5 호 2014. 10. DSMC 방법을이용한로켓플룸의해석 55 : molecular velocity vector : molecular diameter : normalized velocity distribution function in velocity space : number of real molecules represented by a single DSMC molecule : force vector : Knudsen number : characteristic length : number density : position vector : mean free path : total collision cross-section 매우심하게변화한다. 높은고도에서와같이기체의밀도가매우낮은영역에서는기체를연속체라고가정하여그거동을모사한 Navier-Stokes 방정식과같은전형적인해석방법이유효하지못하게되어주로통계역학적인방법인 Bird가개발한 Direct Simulation MonteCarlo(DSMC) 방법을사용하여유동을해석한다. 인공위성의추력기에서진공으로사출되는배기플룸의해석에도 DSMC 방법이주로이용되고있다 [3,4]. 본연구에서는 DSMC 방법을사용하여고고도와저고도를대표하는조건에서의오리피스를통해사출되는플룸의유동을계산하여비교하였고, 실제로켓의노즐출구의조건을적용한플룸영역의계산을실시하였다. 1. 서론로켓은고온의배기플룸을사출함으로써운동량보존의원리에의하여추진력을얻는다. 그렇지만, 이러한플룸의형상은일정하지않고로켓의고도가상승함에따라급격하게팽창한다. 또한로켓의추력을정확하게해석하기위해서는고도에따른플룸의형상변화를정확하게예측해야한다 [1,2]. 또한플룸이로켓저부면에충돌하는현상과플룸으로의열전달을파악하기위하여서도플룸형상에대한정확한예측은필요하다. 이와더불어, 최근에는로켓의피탐지성능이강조됨에따라플룸에서방사되는복사강도의예측역시중요한설계요소의하나이기때문에고도에따른플룸형상의정확한해석의중요성이더욱강조되고있다. 이러한로켓플룸의형상변화는고도상승에따른대기의급격한압력과밀도하강이원인이다. 성층권에해당하는고도 20 km의경우에는대기밀도가해수면에서의약 0.07배, 압력은약 0.05배로떨어지고, 고도가 80 km에이르게되면, 밀도는해수면의약 1.6 10-5 배, 압력은 1.0 10-5 배로매우급격하게감소한다. 지표면근처에서의전형적인플룸구조와고고도에서의플룸구조의개략도를 Fig. 1에나타냈는데, 이와같이고도에따른플룸의구조는 2. DSMC 해석 높은고도에서는대기의밀도가낮기때문에특성길이 (characteristic length) 에대한분자의평균자유행로 (mean free path) 로정의되는 Knudsen 수, Kn (Eq. 1) 가해수면고도에서보다상당히커지게되어이때대기의유동은연속체라고가정할수없고희박유동 (rarefied flow) 이된다. (1) Fig. 2에나타낸바와같이 Eq. 2와같은 Boltzmann 수송방정식은큰 Knudsen 수영역에서도유동을모사할수있다 [5]. Fig. 1 Detailed plume structures at (a) low- and (b) high-altitude[2].
56 전우진 백승욱 박재현 하동성한국추진공학회지 (2) 그렇지만방정식의우변에있는충돌항으로인하여충돌을고려하지않아도되는극히제한적인경우를제외하고는 Boltzmann 방정식을직접적으로푸는것은대단히어렵다. 따라서희박유동장을계산하는문제에있어유동장내에있는기체의분자를직접모사하는 DSMC 법이일반적으로적용되고있으며 DSMC 법을사용하여얻은해가 Boltzmann 방정식을풀어서구한해와동일하다는것이이미수학적으로증명된바있다 [6,7]. DSMC 방법에서는다수의실제입자를하나의모사입자 (simulated particle) 로대표하며하나의모사입자가대표하는실제입자수는경우에따라적절하게조정된다. 먼저격자계를구성한후, 주어진초기조건에따라유동장내에모사입자를생성한다. 모사입자를이동시킨후입자의새로운위치를계산한다. 이때모사입자와경계영역 (boundary) 사이의상호작용도고려한다. 다음으로분자모델을적용하여분자사이의충돌과, 충돌시입자사이의에너지교환및내부에너지재분포를계산한다. 이러한분자의이동과충돌의과정을충분히반복하게한후에각셀내에있는개별모사입자가지닌미시적인 (microscopic) 물성치를샘플링 (sampling) 하여유동장의거시적인 (macroscopic) 물성치를구한다. 입자간충돌모델로는분자의직경이분자의상대속도의함수로주어지는 VHS(variable hard sphere) 모델을적용하였고 [8] 충돌쌍 (collision pair) 을샘플링하는방법으로는 Bird의 NTC(No Time Counter) 방법 [5,9] 을적용하였다. 충돌의일부만을비탄성충돌로간주하고나머지는탄성충돌로가정하는 Larsen-Borgnakke 모델 [10] 을사용하여충돌하는분자의에너지교환을모델링했다. 분자간의충돌이비탄성적이면병진 (translational) 에너지와내부 (internal) 에너지는평형분포로부터샘플링된다. 3. DSMC 해석결과 3.1 코드검증본연구를수행하기위하여비정렬 2차원축대칭 DSMC 코드를사용하였으며, 이는 Bird의정렬격자계를사용한축대칭 DSMC 코드인 DSMC2A[5] 프로그램을기반으로작성된프로그램이다. 플룸을해석하기에앞서본프로그램을검증하기위하여저밀도노즐에대하여결과가공개되어있는 Rothe[11] 의실험결과를사용하였다. Rothe 노즐의형상은 Fig. 3과같다. 사용된기체는질소이고, 노즐입구의경계조건으로전온도 300 K, 전압력 474 Pa이다. 또한배압은 1.5 Pa로이때노즐목직경을기준으로한 는 0.023이다. Fig. 4는노즐대칭축을따른밀도분포이다. x-축은노즐목부터노즐출구까지노즐목의반경으로정규화된축방향거리이고, y-축 Fig. 2 The Knudsen number limits on the mathematical models[4]. Fig. 3 Rothe's nozzle[10].
제 18 권제 5 호 2014. 10. DSMC 방법을이용한로켓플룸의해석 57 Fig. 4 Density variation along axis. Fig. 5 Temperature variation along the axis. 은 로정규화된밀도이다. 계산결과가실험결과와잘일치되고있음을확인할수있다. Fig. 5는노즐중심축을따른온도분포이고 로온도값들을정규화하였다. 계산결과는측정결과와잘일치하고있으며, 급격한팽창으로인하여축을따라서회전온도는병진온도보다높게나타난다. 3.2 오리피스에서사출되는플룸 3.2.1 해석조건오리피스를통한플룸사출의경우에대하여 2차원축대칭 DSMC 방법을사용한계산을수행하였다. 저고도와고고도를대표하는두가지경우를가정하고저고도를대표하는경우에는출구전압에대한배압의비 (P b /P o ) 를 0.1로, 고고도를대표하는경우에는압력비를 0.001로설정하였다. 연소생성물인 CO 2 와 H 2 O가질량비각각 0.55, 0.45로균일하게혼합되어있고, 오리피스출구전압이 1,000 Pa이며, 온도가 300 K인기체가직경 2 cm인오리피스를통하여온도가 300 K, 압력이 100 Pa( 저고도 ), 1 Pa( 고고도 ) 인외부로사출된다. 경계조건으로는대칭축에는모멘텀이보존되는정반사조건을부여하였으며, 노즐출구를포함한외부의열린경계에는자유류조건을부여하였다. 계산에사용된격자계는축방향으로 0.2 m 반경방향으로 0.12 m이고 28,405개의노드와 56,188개의셀로구성되어있다 (Fig. 6). Fig. 6 Grid system for orifice plume flow. 3.2.2 결과정상상태에서의모사입자수는압력비가 0.1인경우에는 89,000여개, 0.001인경우에는 6,300여개다. 압력비에따른밀도분포를 Fig. 7에비교하여나타내었다. 배압이출구전압의 0.1인경우 ( 저고도 ) 에는 Mach disc 등의불연속성같은복잡한유동구조를포함한과소팽창 (underexpanded) 한플룸이관찰되었고 Fig. 8 (a) 에서나타난바와같이주로축방향의방향성을가지고플룸이이동하였다. 배압이출구전압의
58 전우진 백승욱 박재현 하동성한국추진공학회지 Fig. 7 Density profile : (a) P b /P o,exit =0.1; (b) P b /P o,exit =0.001. Fig. 8 Stream line : (a) P b /P o,exit =0.1; (b) P b /P o,exit =0.01. 0.001인경우에는특별한유동구조를보이지않고단순하게팽창하였으며, Fig. 8 (b) 와같은유선을보였다. 이경우오리피스에서사출된플룸은 Fig. 8 (a) 에비하여플룸이거의모든방향으로팽창하며, 일부는 90 도이상으로꺾여역류하기도한다. 이와같은유동경향은 Fig. 1 (b) 에서간략하게도시된고고도에서의플룸팽창에대한경향과일치하고있다. 3.3 로켓노즐플룸 3.3.1 해석조건연소실출구전압이 172 MPa, 전온도가 2,700 K 이고, CO 2 와 H 2 O가질량비각각 0.55, 0.45로균일하게혼합되어있는기체가출구직경 6.14 cm 의노즐을통해서사출되는로켓의배기플룸에대하여 2차원비정렬축대칭 DSMC 해석을수행하였다. 이때노즐외부는고도 20 km인대기조건으로압력은 5.5 kpa, 온도는 217 K이고전압에대한배압의압력비는약 3.2 10-4 이다. 이때대기의밀도는 0.089 kg/m 3 으로해수면에서의밀도보다 0.07배만큼작다. 노즐내부는기체의압력과밀도가높아명백한연속체로볼수있으므로노즐내부유동은연속체전산유체 Fig. 9 Grid system for rocket plume flow. 역학을사용하여노즐의출구조건을계산하였고, 이렇게계산한노즐출구단면의밀도, 속도, 온도분포값을 DSMC 해석을위한배기플룸의입구조건으로사용하였다. 노즐출구에서의조건은압력이평균약 220 kpa, 밀도가평균약 0.5 kg/m 3, 온도가평균약 1,400 K이다. 계산영역은반경방향으로 0.5 m 축방향으로 6 m로잡았으며, 사용된노드는 7,911개, 셀은 15,090개다 (Fig. 9). 경계조건은축과로켓동체
제 18 권제 5 호 2014. 10. DSMC 방법을이용한로켓플룸의해석 59 Fig. 10 DSMC solution of rocket plume at altitude 20 km : (a) density; (b) overall temperature; (c) translational temperature; (d) rotational temperature; (e) Mach number; (f) cell Kundsen number. Fig. 11 Temperature contour comparison: (top) DSMC; (bottom) CFD. (under-expanded) 된플룸이생성되었다. Fig. 10 에밀도, 온도, 회전온도병진온도, 마하수, 셀 Knudsen 수분포를축방향으로 3 m까지만확대하여해석결과를나타냈다. 국소 Knudsen 수는분자의평균자유행로 () 를해당셀면적의제곱근으로나누어구하였는데, 평균자유행로는 Eq. 3에의하여계산되었다. 표면을제외한모든경계를자유류조건으로설정하였으며, 매시간단계에노즐출구에서는노즐유동해석에의해얻은결과에해당하는모사입자의숫자분포가출구단면을따라투입되고나머지자유류경계에서는고도 20 km에서의밀도에해당하는수의모사입자가입사된다. 대칭축에는정반사조건을적용하였고, 로켓동체표면에서의경계조건은난반사 (diffusive reflection) 모델을적용하였다. 또한로켓표면의온도는대기온도와동일하게부여하였다. 3.3.2 결과정상상태에서의모사입자수는 10만 2천여개다. 노즐출구에서의압력이약 220 kpa로대기의 5.5 kpa보다훨씬크기때문에과소팽창 [5] (3) 플룸의직경은약 0.5 m까지노즐출구직경의 8배이상팽창하였으며축을따라 6 m이상진행하였다. 또한, Mach disc나 barrel shock 등의불연속적인유동패턴이반복적으로나타나는것이관찰된다. 고도 20 km의대기는연속체영역에포함되기때문에동일고도에서동일한조건으로연속체해석을수행하여 DSMC 결과와비교하였다. Fig. 11에플룸의온도분포나타내었으며축상단이 DSMC 결과하단이연속체 CFD 결과이다. 경계조건을부여하는방식이나기체의물성치를적용하는방식이상이하여 Mach disc의위치에다소차이가있으나플룸이팽창한크기나내부의구조, 온도분포경향이매
60 전우진 백승욱 박재현 하동성한국추진공학회지 우유사하였으며, 이프로그램을사용하여계산한결과가타당함을보여준다. 특이한점은본고도에서대기의국소 Kn는 1.5 10-5 정도로연속체유동의영역에해당하지만급격한팽창으로인하여 0.5 m 지점이후에서는 0.001 이상의연속체와천이영역경계정도에해당하는 Kn가나타나기도하였는데, 이는연속체해석에영향을줄가능성이있으므로추후이부분에대한확인이필요할것으로사료된다. 4. 결론본연구에서는비정렬격자계를사용한 2차원축대칭 DSMC 방법으로저고도와고고도로대표되는높은배압과낮은배압으로의오리피스플룸의유동과고도 20 km조건에서의로켓플룸유동을해석하였다. 1. P b /P o =0.1인오리피스플룸의경우 Mach disc 등이포함된복잡한유동구조를형성하였지만 P b /P o =0.001인경우에는단순팽창하는유동만을보였다. 2. P b /P o =0.1일때, 주로축방향의유동을보였으나, P b /P o =0.001인일때에는유동의반경방향성분이크게증가하였으며일부는 90도이상으로심하게꺾기기도하였다. 3. 고도 20 km에서의로켓플룸은플룸내에형성되는 barrel shock과 Mach disc에의한반복적인패턴이나타나는복잡한유동구조가잘나타났다. 4. 고도 20 km에서대기는 Kn이작은연속체영역에해당하지만조건에따라서팽창으로인한국부적인천이영역이플룸내부에형성될수있음을확인하였다. 후기본연구는국방과학연구소국방개별기초연구사업의지원으로수행되었으며, 이에깊이감사드립니다. References 1. Sutton, G.P. and Biblarz, O., Rocket Propulsion Elements, 7th ed., John Wiley & Sons, Inc., New York, N.Y., U.S.A., 2001. 2. Simmons, F.S., Rocket Exhaust Plume Phenomenology, The Aerospace Press, El Segundo, CA., U.S.A., 2000. 3. Park, J.H., Kang, S.J., Kim, J.S., Baek, S.W. and Yu, M.J., "DSMC Analysis of Satellite Thruster Plume," Journal of the Korean Society for Aeronautical and Space Sciences, Vol. 29, No. 8, pp. 111-118, 2001. 4. Lee, K.H., Lee, S.N., Yu, M.J., Kim, S.K. and Baek, S.W., "Combined Analysis of Thruster Plume Behavior in Rarefied Region by Preconditioned Navier-Stokes and DSMC Methods," Transactions of the Japan Society for Aeronautical and Space Sciences, Vol. 52, No. 177, pp. 135-143, 2009. 5. Bird, G.A., Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Carlendon Press, Oxford, Oxfordshire, U.K., 1994. 6. Bird, G.A., "Direct Simulation and the Boltzmann Equation," Phyisics of Fluids, Vol. 13, No. 11, pp. 2676-2681, 1970. 7. Wagner, W., "A Convergence Proof for Bird's Direct Simulation Monte Carlo Method for the Boltzmann Equation," J. Statistical Physics, Vol. 66, Nos. 3/4, pp. 1011-1044, 1992. 8. Bird, G.A., "Monte Carlo Simulation in an Engineering Context," Progress in Astonautics and Aeronautics, Vol. 74, pp. 239-255, 1981. 9. Bird, G.A., "Perception of Numerical Methods in Rarefied Gas Dynamics," Progress in Astronautics and Aeronautics, Vol. 118, pp. 211-226, 1989. 10. Borgnakke, C. and Larsen, P.S., "Statistical
제 18 권제 5 호 2014. 10. DSMC 방법을이용한로켓플룸의해석 61 Collision Model for Monte Carlo Simulation of Polyatomic Gas Mixture," J. Comput. Phys., Vol. 18, Issue 4, pp. 405-420, 1975. 11. Rothe, D.E., "Electron-Beam Studies of Viscous Flow in a Supersonic Nozzle," AIAA J., Vol. 9, No. 5, pp. 804-810, 1971.