1. 서 론

차세대 디젤 엔진, 로켓 엔진, 및 초임계 이산화탄소 발전시스템의 운전영역은 연료, 산화제, 또는 작동유체의 임계압력 이상의 압력과 임계온도 이상의 온도를 포함한다. 이러한 영역에서 유체는 초임계 상태로 액체의 밀도를 가지며 분자확산에 의해 혼합되는 기체의 특성을 동시에 가지게 된다. 따라서 발전시스템의 효율적인 설계를 위해서는 초임계 유체의 정확한 수치해석 및 분석이 요구된다. 그러나 초임계 유체의 수치해석 기술은 보편적인 기체의 수치해석 기술보다 상대적으로 미성숙한 상태이며 많은 부분에서 개선과 발전이 요구되고 있다.

초임계 유체는 밀도를 포함한 열역학적 물성치가 온도와 압력의 작은 변화에 큰 변동을 보이는 비선형적 특성을 보이고 있다. 이러한 비선형적 특성을 반영하기 위해서는 실제유체방정식을 적용해야만 한다. 열역학적 물성치의 비선형성은 수치해의 오차에도 영향을 주기 때문에, 기존의 이상기체방정식과 함께 기체 유동을 표현하는 데 어려움이 없었던 수치해석 기법이라도 초임계 유체의 유동을 표현하기 위해서는 수정이 필요하다.

Lacaze 등[1]은 초임계 유체 유동의 해석에서 밀도 차이가 큰 유체의 경계면에서 비정상적인 압력 섭동이 발생하는 것을 언급하며 세 가지 수치해석 기법(에너지 기반, 압력 기반, 엔탈피 기반)을 사용하여 속도장 및 압력장이 수치해석 방법에 따라 어떠한 영향을 받는지 보고하였다. Kitamura와 Shima[2]는 초임계 유체를 해석하기 위해 기존의 SLAU2(Simple Low-dissipation Advection Upstream spliting method 2)를 수정하였다. 에너지 보존방정식을 압력의 수송방정식으로 대체한 이 수치 기법은 압력의 비정상적 섭동 없이 초임계 유체 유동을 해석할 수 있었다.

Engine Combustion Network[3]의 Spray A는 도데케인을 전형적인 디젤엔진의 온도와 압력 조건의 정적연소실에 분사한 벤치마크로서 연소 모델 및 수치 모사를 검증하기 위해 제안되었다. 도데케인은 물리화학적 특성이 디젤과 유사하며 검증된 화학반응 모델이 존재하기 때문에 Spray A의 연료로 사용되었다. Spray A의 조건에서 연소실의 압력과 온도는 각각 6.0 MPa과 900 K이다. 이는 연료와 분위기의 임계압과 임계온도보다 큰 값이다. 따라서 연료는 연소실 내부에서 초임계 유체로 천이할 것이다. Lacaze 등[4]은 열역학적 상태를 분석하여 Spray A의 조건에서 연료 혼합이 액체의 미립화와 증발이 아닌 난류 혼합에 의해 지배되는 초임계 유체의 특성을 보일 것이라는 가설을 제시했다. 이는 분무 및 미립화 모델을 적용하여 연료 혼합과 연소 과정을 모사하는 것보다 실제유체방정식을 사용하여 난류 혼합으로 연료의 혼합 과정을 모사하는 것이 더 타당하다는 주장이다.

본 연구에서는 초임계 도데케인의 연소를 수치 모사하기 위해 오픈소스 유동 해석 라이브러리인 OpenFOAM[5]을 활용하여 압력기반 SLAU2 기법을 구현하였다. 더불어 열역학적 물성치와 화학반응의 효과를 반영하기 위해 실제유체방정식을 사용한 Flamelet Progress Variable (FPV) 모델[6]을 적용하였다. FPV 모델은 혼합율과 반응진행변수의 정보를 통해서 열역학적 물성치와 화학종의 농도에 관한 정보를 얻을 수 있기 때문에 효율적인 반응 유동 계산이 가능하다. 제안된 수치해석 기법을 사용하여 도데케인 제트의 혼합과 연소를 수치 모사하였으며 결과의 분석을 통해 제안된 수치 기법의 타당성을 검증하였다.

2. 수치해석 방법

2.1 압력기반 SLAU2 기법

압력기반 SLAU2 기법은 기존의 SLAU2에서 에너지 보존식 대신 압력의 수송방정식을 해석함으로써 비정상적 압력 섭동 없이 초임계 유체의 모사가 가능하도록 하였다. 압력기반 SLAU2 는 다상유동[7] 및 자기수력학적 유동[8]의 수치 모사에 사용되어 그 성능을 입증하였다. Kitamura와 Shima[2]에 의해 제안된 압력의 수송방정식은 다음 관계식으로부터 유도되었다.

질량 보존식과 에너지 보존식 그리고 열역학적 관계식을 이용하여 정리하면 아래와 같은 압력에 대한 수송방정식이 유도된다.

여기서 질량과 모멘텀의 선속 계산은 기존의 SLAU2 방법과 차이가 없으며 압력의 수송방정식에서 선속은 검사체적의 면방향의 속력 및 속력과 압력의 곱을 다음과 같이 근사하여 계산하였다.

2.2 Flamelet Progress Variable 모델

초임계 연소를 효율적으로 계산하기 위해 실제유체방정식에 기반한 FPV 모델을 사용하였다. FPV는 상압의 비예혼합연소 문제에 적용되어 점화 및 소화를 성공적으로 모사하였기 때문에 본 연구에서 연소 모델로 채택하였다. FPV 모델에 대한 보다 자세한 내용은 [6]에 상술되어 있다. 실제유체 기반의 화염편 모델의 구현은 Kim 등[9]을 참고하였다.

Fig. 1은 여러 삼차방정식을 사용하여 6.0 MPa에서 도데케인의 밀도를 계산한 것을 NIST 데이터[10]와 비교한 것이다. 혼합물의 임계 물성치는 Ribert 등[11]의 혼합 모델을 사용하여 계산하였다. Soave-Redlich-Kwong 및 Peng-Robinson 상태방정식을 사용하여 계산된 밀도는 계측값과분명한 차이를 보이며 특히 저온 영역에서 유체가 액체의 특성을 나타낼 때 더 큰 오차가 발생했다. RKPR(generalized Redlich Kwong-Peng Robinson) 상태방정식[12]을 사용하여 예측한 밀도가 NIST 데이터에 가장 근사한 결과를 보여 본 연구에 사용되었다.

Fig. 1.

Comparison between calculated dodecane density values using three different cubic equations of states. Symbols are NIST data [10]].

2.3 수치 기법 및 모델

SLAU2 기법과 열역학 물성치 라이브러리는 Open- FOAM6.0[3]의 유한체적 압축유동 해석코드인 rhoCentral- Foam을 수정하여 구현되었다. 대류항 계산을 위해서 2차 공간 정확도를 가지는 TVD 기법이 사용되었으며 확산항은 중앙차분법을 이용하여 이산화하였다. 시간에 대한 적분은 2차 정확도의 Runge-Kutta방법을 사용하여 계산하였다. 검증과정에 사용된 FPV 라이브러리는 일정한 압력(6.0 MPa)에서 계산되었다. 루이스 수는 1로 가정하였다. FPV 모델은 화염편 계산 결과를 바탕으로 PDF 적분을 통해 라이브러리를 생성한다. 본 계산에서는 혼합율의 평균과 분산 그리고 반응진행변수의 평균에 대한 3차원 라이브러리(100*25*100)를 사용하였다. 반응 조건에서의 FPV 라이브러리 계산은 Yao 등[13]의 화학반응 메커니즘을 사용하여 계산하였다. Kim 등[14]의 연구를 참고하여 비정상 상태의 화염편은 점화과정을, 정상상태의 화염편은 안정화된 화염을 모사할 수 있도록 라이브러리를 구성하였다. 난류 모델은 RANS $k-\epsilon$ 난류모델을 사용하였다. 기본 모델 상수를 사용할 경우 분사기 근방의 난류 효과를 과대 예측하였기 때문에 모델 상수를 수정하였다.

3. 결과 및 논의

3.1 1D 초임계 유체 대류 문제

Fig. 2는 초임계 유체의 1D 대류 문제를 도식화한 것이다. 계산영역은 5 mm이며 200개의 균일한 크기의 유한체적으로 구성되어있다. x=1 mm-1.5 mm영역에는 액체산소가 분포되어 있으며 나머지 영역에는 수소가 있다. 전 영역에서 초기 속도와 압력은 50 m/s와 100 bar이다. 입구에서 속도는 50 m/s로 고정되어있고 출구에서는 압력을 100 bar로 고정하였으므로 이론상 속도장과 압력장은 시간에 따라서 변하지 않고 50 m/s와 100 bar로 일정하게 유지되어야 한다.

Fig. 2.

Schematics of 1D supercritical fluid advection test.

Fig. 3은 본 연구에서 개발된 코드를 사용하여 계산한 시간에 따른 밀도 분포의 변화를 나타낸 것이다. 속도와 압력은 시간에 따른 변화 없이 50 m/s와 100 bar로 일정했다. 두 유체의 밀도 비가 80에 가까움에도 불구하고 수치적으로 안정적인 해를 구할 수 있었다. 다만 이산화에서 야기된 수치적 확산에 의해 시간이 지남에 따라 경계면에서 밀도의 기울기가 서서히 감소하는 경향을 보였다.

Fig. 3.

Computed density distribution using SLAU2 at t=0, 0.02, and 0.04 ms.

Fig. 4는 OpenFOAM에 구현한 압력기반 SLAU2와 상용 유동해석 코드인 Fluent[15]의 AUSM(Advection Upstream Splitting Method)[16]을 사용하여 계산한 1차원 대류문제의 결과를 함께 나타낸 것이다. 그래프는 시간 0.02 ms에서 속도와 압력으로 이론상 초기값과 같은 결과들을 보여야 한다. 압력기반 SLAU2를 사용하여 계산한 속도장의 경우 계산영역 전체에서 50 m/s의 일정한 값을 보였다. 반면에 AUSM을 사용하여 얻은 계산 결과는 유체의 밀도 차이가 큰 경계면에서 속도가 비정상적으로 크게 변하는 경향을 보였다. 압력 역시 압력기반 SLAU2는 이론값에 일치하는 결과를 보여주지만 AUSM은 비정상적인 압력의 섭동을 나타내었다.

Fig. 4.

1D advection simulation results: velocity (top) and pressure (bottom) distributions at t = 0.02 ms.

에너지 보존식을 통해 문제를 해석할 경우 경계면의 에너지가 상태방정식의 비선형성에 의해 정확하게 표현되지 않는다. 그로 인해 에너지의 오차가 발생하게 되고 이 오차에 의해 압력($p=p(\rho ,e)$)의 비정상적 섭동이 발생하게 된다. 상태방정식의 비선형성이 크지 않은 경우 경계면에서의 비정상적 압력 섭동이 큰 문제가 되지 않는다. 즉 비정상적 압력 섭동은 밀도차가 큰 유체의 경계면에서 비선형성이 큰 상태방정식을 사용해 밀도와 에너지를 계산할 경우 두드러진다[17]. 압력기반 SLAU2가 AUSM에 비해 더 정확한 유동해석 결과를 보이는 것은 궁극적으로 해석변수의 차이에 있다. AUSM이 질량, 모멘텀, 에너지 보존방정식을 해석하는 반면 압력기반 SLAU2의 경우 질량, 모멘텀의 보존방정식과 압력의 수송방정식을 해석한다. 하지만 SLAU2 기법이 항상 더 정확하다고 할 수는 없다. 압력기반 SLAU2는 에너지 보존식이 그에 상응하는 비보존형의 압력식으로 치환된 기법이다. 따라서 유동장이 매끄럽지 않은 경우(non-smooth) 에너지가 보존되지 않는다.

3.2 도데케인 제트의 수치 모사

제시된 수치 기법이 실제 문제에 적용되었을 때의 성능을 검증하기 위해서 초임계 도데케인 제트의 수치해석을 수행하였다. 도데케인은 디젤연료와 유사한 열역학적, 화학적 특성이 있어 디젤 엔진의 수치 모사에서 실제 디젤연료를 대신하는 경우가 많다. Sandia 연구소에서 수행한 정적연소실 내의 도데케인 분사[3]를 수치 모사한 후 해석 결과를 계측값과 비교 분석하였다.

실험에서 도데케인은 한 변의 길이가 105 mm인 정육면체의 정적연소실 내부로 분사된다. 도데케인이 분사되는 분사기의 직경은 D=0.09 mm이며 분사압력은 150 MPa이다. 유출계수 0.89를 사용하여 계산한 제트의 유속은 대략 570 m/s 이며 이 값이 수치 모사에서 입구의 경계조건으로 사용되었다. 비반응 및 반응 조건 모두 연소실 압력은 60 bar그리고 온도는 900 K이다. 비반응 조건에서 연소실 조성은 N₂:89.71, CO₂:6.52, H₂O:3.77이며, 반응 조건에서는 O₂:15, N₂:75.15, CO₂:6.23, H₂O:3.62이다. 도데케인의 임계점은 658 K, 18 bar로 연소실 내부 압력은 임계 압력보다 높고 분사 온도는 임계 온도보다 낮다. 도데케인은 액체로 분사된 후 혼합 및 열전달로 인해 초임계유체로 전이할 것이다.

Fig. 5는 RANS 계산에 사용된 계산영역을 도식적으로 나타낸 것이다. 효율적인 계산을 위해 축 대칭 격자를 사용하였다. 계산영역은 축 방향으로 60 mm 반경 방향으로 30 mm이며, 분사기와 대칭축을 제외한 모든 부분은 벽으로 단열과 노슬립 조건을 사용하였다. 격자의 최소크기는 0.05D이다. 격자를 분사기 근방에 조밀하게 분포하여 혼합층에서 속도 및 열역학적 물성치의 급격한 변화를 나타낼 수 있도록 하였다. 두 가지 격자가 수치 모사에 사용되었다. 격자는 각각 12만 개 및 18만 개의 유한체적들로 구성되어 있다. 수치해석 결과 격자에 의한 해석 결과의 차이가 미미하여 18만 개의 유한요소들로 구성된 격자를 사용한 결과가 논의에 사용되었다.

Fig. 5.

Schematic diagram of the RANS computation domain.

3.2.1 비반응 도데케인 제트의 수치 모사

Fig. 6은 예측된 기체 침투 길이의 시간에 따른 변화를 실험값과 함께 나타낸 것이다. 침투 길이는 분사기로부터 중심축을 따라서 도데케인의 질량분율이 0.01인 지점을 나타낸 것이다[4]. 초기에는 예측값과 계측값이 근사하였으나 0.3 ms 이후에는 RANS모델이 침투 길이를 과소예측하는 경향을 보였다. 기체 침투 길이의 증가량은 시간에 따라서 감소하는데 이는 제트와 연소실 내부의 분위기의 혼합에 의해 모멘텀과 에너지가 전달되는 효과에 의한 것이다. 해석 결과는 후류에서 난류혼합 효과를 과대예측한 것으로 사료된다.

Fig. 6.

Measured [3] and predicted vapor penetrations.

Fig. 7은 예측된 축 방향(a)과 반경 방향(b)의 혼합율(또는 도데케인의 질량분율)을 계측값과 함께 나타낸 것이다. 반경방향의 혼합율 분포는 x=25 mm의 결과값이다. 축방향과 반경방향의 결과 모두 주어진 범위 내에서 실험값에 근사하였다. 계측된 위치에서 x/D값이 매우 큼에도 불구하고(x/D>166) 제시된 수치 기법과 모델을 사용하여 초임계 제트에 대한 적절한 예측을 할 수 있었다. 이는 초임계 제트가 완전히 발달할 경우 상압의 실험 결과로부터 제안된 난류모델로 모사할 수 있음을 의미한다.

Fig. 7.

Central(a) and Radial(b) mixture fraction profiles for the non-reacting case. Lines are simulation results and symbols are measurement data [3].

3.2.2 반응 도데케인 제트의 수치 모사

Fig. 8은 반응 조건에서 계산된 시간에 따른 연소실의 최대 온도를 나타낸 것이다. 연료 분사 시각으로부터 단위 시간당 온도의 변화가 가장 큰 때를 착화지연시간이라 정의하였다. 예측된 착화지연시간은 0.4 ms로 실험에서 계측된 값(0.4 ms)과 일치했다. 이는 이상기체방정식을 사용하여 라이브러리를 구성한 선행연구 결과[14] 와도 유사하다. 이러한 이유로는 첫째, 착화반응이 일어나는 영역이 희박영역이기 때문에 상대적으로 온도가 높아(700~900 K) 혼합물의 열역학적 물성치가 상태방정식에 따라 큰 차이가 없기 때문이다. 또한 착화지연시간은 화학반응 메커니즘에 크게 의존하는데 두 연구에서 모두 같은 화학반응 메커니즘[13]을 사용하였다.

Fig. 8.

Temporal evolution of the peak temperature inside the chamber.

Fig. 8에서 연료 분사 직후에는 상대적으로 느린 발열반응이 진행되어 최대온도가 1100 K까지 상승하는데 대략 0.35 ms가 소요된다. 하지만 온도가 상승함에 따라 반응속도가 급격히 증가하여 1100-2000 K구간에 소요되는 시간은 0.15 ms이다. 일반적으로 도데케인 점화과정은 저온 연소와 고온 연소 두 단계로 구분한다. 그리고 각각의 단계가 활발한 영역은 특정 화학종의 분포를 통해 확인할 수 있다.

Fig. 9는 반응조건에서 계산된 온도, CH₂O, 그리고 OH의 분포를 나타낸 것이다. 화염이 본질적으로 비정상적이므로 상류의 화염이 정상상태에 이른 것으로 예상되는 1.0 ms에서의 결과를 나타내었다. 도데케인의 연소에서 저온 연소가 일어나는 영역은 CH₂O의 분포로 그리고 고온 연소가 일어나는 영역은 OH의 분포로 유추할 수 있다. 대칭축 근방에서 CH₂O는 저온 연소로 인해 생성되지만 고온 연소의 부재로 OH는 분포하지 않는다. 반면 고온 영역에서는 CH₂O의 농도가 급격히 감소하고 국소적으로 OH의 농도가 증가한다. 고온영역과 OH 농도의 관계는 주요 발열반응인 CO+OH$\iff$CO₂+H 에 기인한 것이다.

Fig. 9.

Temperature(a), CH₂O mass fraction(b), and OH mass fraction(c) distributions for the reacting case at 1.0 ms.

수치해석으로부터 얻어진 화염 부상 높이는 16.8 mm로 실험에서 계측된 화염 부상 높이인 16.5 mm에 근사하였다. 부상 높이는 유동과 화학반응의 상호작용에 의해 결정되므로 이 결과는 제안된 수치해석 기법이 초임계 유동과 화학반응의 상호작용을 효과적으로 고려할 수 있음을 반영한다. 그러나 초임계 조건의 반응 유동 해석은 난류 모델과 화학반응 모델의 높은 불확실성을 가지고 있어 다양한 조건에서 실험과 수치해석의 비교검증을 통한 지속적인 개선이 필요할 것으로 사료된다.

4. 결 론

강건하고 효율적인 초임계 반응 유동 해석을 위한 수치해석 프로그램이 OpenFOAM 라이브러리를 활용하여 구현되었다. 프로그램은 초임계 유체의 안정적이고 정확한 해석을 위해 압력기반 SLAU2 기법을 사용하였으며 연소반응을 효율적으로 나타내기 위해 실제유체방정식 기반의 FPV 모델을 사용하였다. 구현된 프로그램은 1D 초임계 유체 대류 문제와 도데케인 제트 분사의 수치 모사에 적용되어 그 성능을 검증하였다. 기존의 에너지 보존식 기반의 해석기법을 사용하여 1D 초임계 유체 대류를 수치 모사할 경우 비정상적 압력 섭동이 발생하였으나 본 연구에서 구현한 압력기반의 수치 기법에서는 그러한 비정상적 압력 섭동이 발생하지 않았다. 제안된 수치 기법은 높은 밀도 비를 가지는 유체 경계면의 수송과정을 안정적으로 모사할 수 있었다. 수치해석 방법은 디젤연료의 분사와 연소의 벤치마크로 잘 알려진 ECN의 Spray A 조건을 모사하여 그 정확도와 활용 가능성을 검증하였다. 수치 모사를 통해 얻어진 기체 침투 길이, 연료의 공간적 분포, 착화지연 시간, 그리고 화염 부상 높이는 실험에서 계측된 값에 근사하였다. 이를 통해 개발된 수치해석 프로그램이 초임계 반응 유동의 기본적인 특성을 구현할 수 있음을 확인하였다. 그러나 화학반응 모델과 난류 모델이 초임계 조건에서 많은 불확실성을 내재하고 있기 때문에 다양한 조건에서 실험과 계산의 비교 검증을 통한 개선이 필요할 것으로 사료된다.

기호설명

c: Speed of sound

D : Injector diameter

e : Internal Energy

p : Pressure

T : Temperature

u : Velocity vector

Z : Mixture fraction

κ : Thermal conductivity

ρ : Density

τ : Viscous stress tensor