1D 네트워크 모델을 이용한 항공용 가스터빈 연소기에서의 음향장 해석

1. 서 론

최근 배출가스 규제가 심화됨에 따라 배출 및 연료 소비율의 감소를 모두 만족시키는 희박 예혼합 연소시스템(Lean Premixed Combustion System)이 주목을 받고 있다. 그러나 이와 같은 연소 시스템은 연소가 희박 가연 한계에서 일어나기 때문에, 외부의 작은 인자에도 당량비 또는 공기 및 연료의 유량 변화가 일어나게 되고, 이로 인해 연소기 내부의 열방출파의 변화를 유도한다. 이와 같은 열방출파의 변동은 다시금 연소기 내부의 압력파 변동을 일으키게 되고, 최종적으로 이 변동들 간의 상호작용으로 인해 연소불안정(Combustion Instability) 현상이 발생하며, 이것은 반드시 해결되어야 하는 큰 문제점으로 제기되고 있다. 그렇기 때문에 연소불안정 현상에 대한 정확한 이해와 그에 대한 정확한 특성을 예측할 수 있는 기술 개발이 요구된다[1-3].

연소불안정 현상을 예측하는 방법에는 여러 기법들이 존재하는데, LES(Large Eddy Simulation), RANS(Reynolds Averaged Navier Stokes) 등의 경우에는 질량, 운동량, 에너지 보존 방정식을 이용하여 연소장과 음향장을 직접 연립하여 풀고 유동장까지 해석하게 되므로, 시간과 비용면에 있어서 비교적 큰 소비를 하게 된다. 반면 본 연구에 적용된 열음향 해석 모델(Thermoacoustic Analysis Model, TA Model)은 실제 연소 시스템을 간소화 시켜 유동장 및 연소장을 직접 풀지 않고, 보다 빠르고 단순하게 음향장을 해석하는 기법으로서 앞선 방법들에 비해 시간과 비용을 절약할 수 있고, 또한 실제 해석 결과에 있어서 신뢰도가 확보되어 널리 사용되고 있다[4,5].

본 연구에서는 열음향 해석 모델 중, 1D 네트워크 모델(1D Network Model)을 개발 및 연구에 적용하였고, 본 해석 모델은 간단한 다단 덕트의 형상뿐만이 아니라, 실제 환형 연소기의 형상과 같이 다중경로를 포함한 복잡한 형상에도 적용 가능하다.

본 논문은 개발된 1D 네트워크 모델링 기법의 소개와 더불어 실제 개발단계에 있는 환형 연소기의 형상에서 사전 검증된 3D 유한 요소 해석(Finite Element Method) 기법을 기반으로 하는 Helmholtz solver[6,7]의 음향 해석 결과와 비교하여, 본 모델의 해석 가능성을 확인하였다.

2. 해석 모델 및 방법

2.1 1D 네트워크 모델링

Fig. 1은 가장 일반적인 환형연소기의 반경방향으로 매우 작게 분할된 형상의 개략도를 나타낸 것으로서, 여기서의 압력 섭동, 길이방향 및 접선방향의 속도 섭동 그리고 밀도 섭동은 아래와 같이 표현될 수 있다.

Fig. 1

Diagram of a thin annular sector with an area change for the thermoacoustic network model.

$\[ p^{\mathrm{\text{'}}}\left(x,\theta ,t\right)=A_{\pm }\exp \left(i{k}_{\pm }+in\theta +i\omega t\right)=\widehat{p}\exp \left(in\theta +i\omega t\right) \]$

(1)

$\[ {\rho }^{\mathrm{\text{'}}}\left(x,\theta ,t\right)=\frac{1}{\overline{c^2}}{A}_{\pm }\exp \left(i{k}_{\pm }+in\theta +i\omega t\right)=\widehat{\rho }\exp \left(in\theta +i\omega t\right) \]$

(2)

$\[ u^{\mathrm{\text{'}}}\left(x,\theta ,t\right)=-\frac{k_{\pm }}{\overline{\rho }{\alpha }_{\pm }}{A}_{\pm }\exp \left(i{k}_{\pm }+in\theta +i\omega t\right)=\widehat{u}\exp \left(in\theta +i\omega t\right) \]$

(3)

$\[ w^{\mathrm{\text{'}}}\left(x,\theta ,t\right)=-\frac{n}{R\overline{\rho }{\alpha }_{\pm }}{A}_{\pm }\exp \left(i{k}_{\pm }+in\theta +i\omega t\right)=\widehat{w}\exp \left(in\theta +i\omega t\right) \]$

(4)

여기서 A는 전달되는 음파의 크기를 나타내며, 상첨자 ‘ 및 ^는 각각 시간 및 Fourier 영역에서의 섭동량을 의미한다. 그리고 하첨자 ±는 음향장에서의 상류방향과 하류방향을 의미하고, n과 k는 각각 원주방향 및 길이방향의 파수를 의미한다. 또한 $$ \overline{c}$$는 평균 음속이며, α_±및 k₀, k_±는 아래와 같다.

$\[ {\alpha }_{\pm }=\omega +\overline{u}{k}_{\pm } \]$

(5)

$\[ k_0=-\frac{\omega }{\overline{u}} \]$

(6)

$\[ k_{\pm }=\frac{\overline{M}\omega \mp {\left({\omega }^2-{\omega }_c^2\right)}^{\frac{1}{2}}}{\overline{c}\left(1-\overline{M^2}\right)} \]$

(7)

$\[ {\omega }_c=\frac{n\overline{c}}{R}{\left(1-\overline{M}\right)}^{\frac{1}{2}} \]$

(8)

이때 Fig. 1에서 보는 것과 같이 주된 해석대상의 연소기의 형태가 반경방향(r-direction)의 특성 길이가 길이방향(x-direction) 및 원주방향(θ-direction) 대비 짧은 형상을 갖고 있으므로, 반경 방향에 대한 공진 모드는 계측되지 않는 것으로 간주하여, 반경방향의 변화를 무시하고 모델링을 진행하였다. 이와 같은 가정은 모델링 과정 및 해석 과정의 시간과 비용을 크게 단축시킬 수 있다.

본 연구에서의 1D 네트워크 모델은 점성과 열확산을 무시한 이상기체 상태의 질량, 운동량, 에너지에 대한 보존방정식으로 표현할 수 있다. 또한 섭동량이 평균값에 비해 크기가 매우 작다는 가정을 통해 비제차 파동방정식(inhomogeneous wave equation)으로 표현할 수 있다[1,4]. 여기서 Eq. 1-4의 압력, 속도, 밀도에 대한 식을 변환한 형태로서 보존방정식을 표현할 수 있고, 아래와 같은 식으로 구성될 수 있다. 이때 Eq. 9-10은 면적변화, Eq. 14-15는 길이변화 상태에서 활용될 수 있는 행렬식이다.

$\[ {\left(\widehat{p},\widehat{u},\widehat{\rho },\widehat{w}\right)}^T=F·W\left(x\right) \]$

(9)

$\[ {\left(\widehat{m},\widehat{f_x},\widehat{f_{\theta }},\widehat{e}\right)}^T=G·{\left(\widehat{p},\widehat{u},\widehat{\rho },\widehat{w}\right)}^T \]$

(10)

이때 행렬 W(x), F, G는 아래와 같다.

$\[ W\left(x\right)={\left(A_{+}{e}^{i{k}_{+}x},A\_e^{ik\_x},\mathrm{0,0}\right)}^T \]$

(11)

$\[ F=\left(\begin{array}{cccc}1& 1& 0& 0\\ -\frac{k_{+}}{\overline{\rho }{\alpha }_{+}}& -\frac{k_{-}}{\overline{\rho }{\alpha }_{-}}& 0& \frac{\eta \stackrel{-}{u}}{\overline{\rho }\overline{c^2}}\\ \frac{1}{\overline{c^2}}& \frac{1}{\overline{c^2}}& -\frac{1}{\overline{c^2}}& 0\\ -\frac{\eta }{R\overline{\rho }{\alpha }_{+}}& -\frac{\eta }{R\overline{\rho }{\alpha }_{-}}& 0& \frac{k_{0}R}{\overline{\rho }\overline{c}}\end{array}\right) \]$

(12)

$\[ G=\left(\begin{array}{cccc}0& \overline{\rho }& \overline{u}& 0\\ 1& 2\overline{\rho }\overline{u}& \overline{u^2}& 0\\ 0& 0& 0& R\overline{\rho }\overline{u}\\ \eta \overline{u}& \eta \overline{p}+\frac{3\overline{\rho }\overline{u^2}}{2}& \frac{\overline{u^3}}{2}& 0\end{array}\right) \]$

(13)

$\[ W\left(x_0+x_1\right)=P\left(x_1\right)W\left(x_0\right) \]$

(14)

$\[ P\left(x_1\right)=\left(\begin{array}{cccc}{e}^{i{k}_{+}{x}_1}& 0& 0& 0\\ 0& e^{i{k}_{-}{x}_1}& 0& 0\\ 0& 0& e^{i{k}_0{x}_1}& 0\\ 0& 0& 0& e^{i{k}_0{x}_1}\end{array}\right) \]$

(15)

여기서 Fig. 1과 같은 면적변화가 있는 형상에서는 질량, 각운동량 그리고 에너지 보존방정식은 일반적으로 만족하지만, Stow 등[8]이 언급한 것과 같이 축방향의 운동량이 증가 혹은 감소하기 때문에, 다음과 같은 변경된 보존 방정식을 적용하여 모델링을 진행하였다.

$\[ A_2{\rho }_2{u}_2^2-A_1{\rho }_1{u}_1^2=A_2\left(p_1-p_2\right) \]$

(16)

이때 하첨자 1, 2는 각각 면적변화 전과 후를 의미한다.

2.2 해석 방법

Fig. 2는 본 연구에서 개발된 모델을 검증하기 위해 해석 대상으로 선정한 연소기의 형상으로서, 한화에어로스페이스에서 실제 개발단계에 있는 25 kN급 엔진에 적용되는 연소기의 형상이다. 기체의 흐름은 노즐 이전에 케이싱과 이후 라이너 및 12개의 노즐로 분리되어 연소실로 흐르고, 이때 노즐에서는 연료와 공기가 혼합되어 연소실에 분사된다. 실제 개발 연소기에서 정확한 음향 경계조건을 측정하기에는 한계가 있는 관계로, 다양한 음향 경계조건에 대한 영향성 평가가 필요하다. 그러나 현재의 초기 연구에서는 케이싱의 입구, 연소 후에 기체가 나가는 연소실의 출구 그리고 연소실과 라이너 사이의 벽면에서의 음향경계조건은 닫힌 경계조건으로 분류하였다.

Fig. 2

Sectional schematics of the aero gas turbine combustor under development.

Table 1은 엔진 사이클 해석을 통하여 도출된 연소기 주요 부위의 물성치로서, 동일한 값들이 현재의 모델링에 적용되었다.

Table 1.

Operating conditions and gas properties.

앞서 설명된 면적변화가 있는 구간에서 사용되는 보존방정식 Eq. 9-10과 수정된 축방향 보존방정식인 Eq. 16을 조합하여, 각 네트워크 모듈의 면적이 감소하는 영역을 행렬식으로 구성하면 아래와 같다.

$\[ \left(\begin{array}{c}\widehat{m_2^n}\\ \widehat{f_{x_2}^n}\\ \widehat{f_{{\theta }_2}^n}\\ \widehat{e_2^n}\end{array}\right)=\left(\begin{array}{ccc}1/12& 0& \begin{array}{cc}0& 0\end{array}\\ E^{*}& F^{*}& \begin{array}{cc}0& G^{*}\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}\begin{array}{cc}1/12& 0\end{array}\\ \begin{array}{cc}0& 1/12\end{array}\end{array}\end{array}\right)·\left(\begin{array}{c}\widehat{m_1^n}\\ \widehat{f_{x_1}^n}\\ \widehat{f_{{\theta }_1}^n}\\ \widehat{e_1^n}\end{array}\right) \]$

(17)

$\[ E^{*}=\left[\frac{{\gamma }_1+\frac{1}{2}\overline{M_1^2}\left({\gamma }_1-1\right)}{\overline{M_1^2}-{\gamma }_1}\right]\overline{u_1}\left(12\beta -1\right) \]$

(18)

$\[ F^{*}=\left[1-\frac{{\gamma }_1\overline{M_1^2}}{\overline{M_1^2}-{\gamma }_1}\right]\left(12\beta -1\right) \]$

(19)

$\[ G^{*}=\left[\frac{\overline{M_1^2}\left({\gamma }_1-1\right)}{\overline{M_1^2}-{\gamma }_1}\right]\frac{1}{\overline{u_1}}\left(12\beta -1\right) \]$

(20)

면적이 증가하는 영역에 해당하는 모듈의 행렬식은 아래와 같다.

$\[ \left(\begin{array}{c}\widehat{m_4^n}\\ \widehat{f_{x_4}^n}\\ \widehat{f_{{\theta }_4}^n}\\ \widehat{e_4^n}\end{array}\right)=\left(\begin{array}{ccc}1/12& 0& \begin{array}{cc}0& 0\end{array}\\ E^{**}& F^{**}& \begin{array}{cc}0& G^{**}\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}\begin{array}{cc}1/12& 0\end{array}\\ \begin{array}{cc}0& 1/12\end{array}\end{array}\end{array}\right)·\left(\begin{array}{c}\widehat{m_3^n}\\ \widehat{f_{x_3}^n}\\ \widehat{f_{{\theta }_3}^n}\\ \widehat{e_3^n}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ 0\\ \widehat{q}\end{array}\right) \]$

(21)

$\[ E^{**}=\left[\frac{{\gamma }_3+\frac{1}{2}\overline{M_3^2}\left({\gamma }_3-1\right)}{\overline{M_3^2}-{\gamma }_3}\right]\overline{u_3}\left(12{\beta }^{*}-1\right) \]$

(22)

$\[ F^{**}=\left[1-\frac{{\gamma }_3\overline{M_3^2}}{\overline{M_3^2}-{\gamma }_3}\right]\left(12{\beta }^{*}-1\right) \]$

(23)

$\[ G^{**}=\left[\frac{\overline{M_3^2}\left({\gamma }_3-1\right)}{\overline{M_3^2}-{\gamma }_3}\right]\frac{1}{\overline{u_3}}\left(12{\beta }^{*}-1\right) \]$

(24)

여기서 β 및 β^*는 면적의 감소 및 증가할 때의 면적변화율을 의미한다.

위의 과정에 있어서 면적이 증가하는 노즐과 연소실 사이의 관계에 있어서 열의 유입 및 열발생 섭동을 고려하는 경우, Eq. 21의 열발생 섭동($$ \widehat{q}$$)항에 실험 및 시뮬레이션을 통한 화염전달함수 데이터를 적용하여 해석을 진행할 수 있다. 그러나 본 연구에서는 열발생 섭동을 고려하지 않았으므로, 열발생 섭동항을 0으로 처리하였다.

앞서 설명된 여러 행렬식을 통해 각각의 네트워크 모듈을 구성한 개략도는 Fig. 3(a)와 같고, 이를 해석 형상으로 구성하면 Fig. 3(b)과 같이 나타낼 수 있으며, 이는 입구, 케이싱, 12개의 노즐, 라이너 그리고 연소실 및 출구로 구성할 수 있다. 이때 복잡한 형상을 갖는 노즐은 실제 노즐과 동일한 면적의 12개의 단일 노즐로 구성하였다.

Fig. 3

(a) Representative network elements and (b) simplified sectional view of the target combustor geometry for the 1D network model.

본 연구에서 개발한 1D 네트워크 모델을 비교검증하기 위해 사용된 3D Helmholtz solver의 주요 서브루틴은 Fig. 4와 같다. 격자계의 구성은 사면체, 프리즘, 피라미드, 육면체의 4가지 요소를 혼합한 비정렬 격자계를 사용하였고, ARPACK(ARnoldi PACKage)에서 제공하는Arnoldi[9] 방법을 사용하여 고유치 문제를 해결하였다. 또한 MUMPS(Multi-frontal Massively Parallel Solver)를 활용하여 대규모 복소 행렬식의 계산 수렴성을 향상시켰다. 이러한 해석 모델에 대해서는 이전 연구[6,7]에서 자세하게 설명하였다.

Fig. 4

Structure and Subroutines of the 3D Helmholtz Solver[6].

3. 해석 결과 및 고찰

Fig. 5는 노즐 및 라이너 등을 제외하고, 오직 연소실만을 고려한 결과로서, 비교적 간단한 형상부터 우선적으로 Helmholtz solver 모델과 비교하여, 모델을 1차적으로 검증하였다. 그림에서 설명한 모드는 첫 번째 원주 방향(1st circumferential)에 해당되는 주파수로서, 1D 네트워크 모델이 643.8 Hz, 3D Helmholtz solver가 644.4 Hz로서 약 0.1%의 오차범위를 갖고 있는 것으로 확인되었다. 또한 음향 모드 형상에서 확인할 수 있듯이, 동일한 원주방향(θ)으로의 모드 형상을 확인할 수 있었다.

Fig. 5

Modeshape calculation results of combustor only using (a) the 1D network model and (b) the 3D helmholtz solver.

Fig. 6과 Fig. 7은 해석 대상 연소기의 전체 형상을 고려하여 본 네트워크 모델과 3D Helmholtz Solver의 음향 해석을 진행한 결과로서, 다양한 모드 중에서 길이 방향 모드(longitudinal mode)와 복합 모드 대비, 상대적으로 낮은 주파수에서 도출된 1C(Circumferential mode)의 계산 결과를 나타낸 것이다. Fig. 5에서의 연소실만 고려하였을 때보다, 외부 라이너 공간이 반영된 관계로 (즉, 원주 방향의 길이 증가), 1C 모드의 주파수가 모두 감소하였음을 알 수 있다. 특히, 이 주파수에서 1D 및 3D 결과 모두 외부 라이너와 케이싱 사이에서 상대적으로 높은 압력 진폭이 존재하는 것으로 나타났다. 그러나 해당 1C에서 1D 모델의 주파수는 553.5 Hz였고, 3D 모델링의 해석 결과는 527.2 Hz로서, Fig. 5에서의 연소실만 해석하였을 때보다 차이가 크게 계산되었다.

Fig. 6

1C modeshape(553.5 Hz) calculation results of the overall geometry using the 1D Network model. (a) 3D plot and (b) 2D plot at the dump plane.

Fig. 7

1C modeshape(527.2 Hz) calculation results of the overall geometry using the 3D helmholtz solver.

Table 2는 1,000 Hz 미만에서 존재하는 3가지 모드에 대하여 1D와 3D 모델링 계산 주파수를 비교한 것이다. 현재의 1D 모델은 순수 길이 방향(1L) 및 원주 방향(1C)와 복합 모드(1C1L)를 모두 예측하는 것은 성공하였으나, 해당 모드의 주파수에 있어서는 3D 해석 결과와 약 5% 수준의 오차가 존재하는 것으로 나타났다.

Table 2.

Comparison of resonance frequencies between the 1D Network model and the 3D Helmholtz solver calculation results.

이러한 오차의 주된 원인으로는 두 모델에서 모두 라이너를 경계로 하는 냉각 공기 통로 및 연소실 사이에서의 두 공간의 음향장간의 상호 영향이 정확하게 반영되지 않았기 때문이라고 판단된다. 실제 연소실과 냉각공기 통로 사이의 라이너에서는 수많은 냉각 공기공들이 존재하며, 이를 통한 음향장의 상호 간섭 및 음향 에너지의 댐핑과 흡수 현상이 발생하게 된다[1]. 그러나 현재의 3D 모델은 라이너가 완전한 벽면으로 연소실 측면에서 닫힌 음향 경계 조건으로 설정되고, 1D 모델에서는 연소실 입구에서 냉각 공기 통로로 분리된 유량이 연소실 후단에서 합쳐지는 형태로 모델이 구성되어 있다 (Fig. 3(b) 참고). 냉각 공기공에 의한 음향파의 전달 및 댐핑의 영향을 반영한 모델 개발이 현재 진행 중에 있으며, 이에 대한 내용은 향후 연구 결과를 통하여 소개할 예정이다.

4. 결론 및 향후계획

본 연구에서는 실제 개발단계에 있는 항공용 가스터빈 연소실의 음향장 해석을 위하여 1D 네트워크 모델을 개발하였고, 3D 기반의 Helmholtz solver의 해석 결과와 비교하였다. 본 해석 모델은 다단의 면적 변화 및 동시에 많은 경로로 구성된 형상까지 해석 가능하도록 개발되었고, 이를 통하여 높은 계산 비용 및 시간이 요구되는 3D 해석을 이용하지 않고도(예, 3D 해석 : 약 1시간 vs. 1D 해석 : 약 1분), 비교적 복잡한 형상을 갖는 시스템에서의 음향장 해석이 가능함을 확인하였다. 현재의 연구는 복잡한 시스템에서의 네트워크 모델 개발에 대한 1차적인 연구 결과물로써, 향후 연구에서는 냉각 라이너와 같은 다공판에서의 유동의 분리 및 합류와 음향장의 댐핑 등을 반영할 수 있도록, 추가적인 신뢰성 확보를 위한 모델 개발을 진행할 예정이다.

Acknowledgments

본 연구는 2018년도 정부(교육부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구 사업(NRF-2018R1D1A3B07044440) 및 산업통상자원부 항공우주부품기술개발사업(10067074)의 지원을 받아 수행된 결과입니다.

References

• Lieuwen, T., and Yang, V., “Combustion Instability in Gas Turbine Engines“, The American Institute of Aeronautics and Astronautics, 1801 Alexander Bell Drive, Reston, VA, U.S.A, (2005).
• Huang, Y., Sung, H.G., Hsieh, S.Y., and Yang, V., “Large-Eddy Simulation of Combustion Dynamics of Lean-Premixed Swirl-Stabilized Combustor“, Journal of Propulsion and Power, 19(5), (2003), p782-794, (2003). [https://doi.org/10.2514/2.6194]
• Ducruix, S., Schuller, T., Durox, D., and Candel, S., “Combustion Dynamics and Instabilities: Elementary Coupling and Driving Mechanisms“, Journal of Propulsion and Power, 19(5), p722-734, (2003). [https://doi.org/10.2514/2.6182]
• Lieuwen, T., “Modeling Premixed Combustion-Acoustic Wave Interactions: A Review“, Journal of Propulsion and Power, 21(4), p591-599, (2003).
• Wolf, P., Balakrishnam, R., Sraffelbach, G., Gicquel, L., and Poinsot, T., “Using LES to Study Reacting Flows and Instabilities in Annular Combustion Chambers“, Journal of Flow, Turbulence and Combustion, 88(1-2), p191-206, (2012). [https://doi.org/10.1007/s10494-011-9367-7]
• Kim, S.K., Kim, D., and Cha, D.J., “Finite Element Analysis of Self-Excited Instabilities in a Lean Premixed Gas Turbine Combustor”, International Journal of Heat and Mass Transfer, 120, p350-360, (2018). [https://doi.org/10.1016/j.ijheatmasstransfer.2017.12.021]
• Lim, J., Kim, D., Kim, S.K., and Cha, D.J., “Effects of Acoustic Boundary Conditions on Combustion Instabilities in a Gas Turbine Combustor“, Journal of the Korean Society of Propulsion Engineers, 19(4), p15-23, (2015). [https://doi.org/10.6108/kspe.2015.19.4.015]
• Stow, S.R., and Dowling, A.P., “Thermoacoustic Oscillations in an Annular Combustor“, ASME Turbo Expo 2001, New Orleans, Louisiana, U.S.A., ASME GT2001-0037, June), (2001. [https://doi.org/10.1115/2001-gt-0037]
• Lehoucq, R., and Sorensen, D., “ARPACK. User’s Guide: Solution of Large Scale Eigenvalue Problems with Implicitly restarted Arnoldi Methods“, retrieved 5 March 2017, from http://www.caam.rice.edu/software/ARPACK/ 8), October, (1997.