Non-equilibrium 2상 유동 모델을 적용한 OTEC 터빈의 수치 해석
Abstract
Carbon neutrality is emerging as an irresistible mission of the times. To realize this, the proportion of renewable energy must be expanded, and ocean energy can be an alternative. OTEC (closed type), one of ocean energy, has the largest theoretical potential among ocean energy. In the case of the closed type of OTEC, two-phase flow may occur in the turbine, and it is necessary to model in CFD (Computational Fluid Dynamics) for more accurate design and analysis. In this study, a non-equilibrium (Non-eq) two-phase model was applied instead of the equilibrium (Eq) model assuming thermodynamic equilibrium. The result was compared and analyzed with single-phase model, two-phase model, and experimental results. In addition, the GCI method was used to estimate the discretization error. As a result of the CFD simulations, the Non-eq model was the most similar to the experimental values compared to the single-phase and Eq model. Qualitatively, the Non-eq model showed some differences in velocity field and mass fraction at spanwise 50 percent. In conclusion, the model estimated to be closest to reality is the Non-eq model, and the rationale for this is estimated quantitatively and qualitatively.
Keywords:
Ocean Thermal Energy Conversion, Turbine, Computational Fluid Dynamics, Two Phase, Non-equilibrium Model1. 서 론
세계는 지금 탄소중립 실현을 위해 다양하고 구체적인 방안을 제시하고 이행하고 있다. 우리나라의 경우 탄소중립 시나리오 최종(안)을 지난 2021년 10월 18일에 정부에서 발표하였다. 본 시나리오에는 배출과 흡수 및 제거의 대분류 하에 전환, 산업, 건물, 수송, 수소와 흡수원 및 CCUS 등으로 부문을 나누어, A안과 B안 총 2개의 방안을 구체적으로 제시하였다. 수소 생산 부문에 있어서 수소 공급량의 대부분은 해외 수입에 의존하게 되는데, 현재 호주, 미국, 사우디 등의 풍력 및 태양광 발전을 통한 그린 수소 수입이 가장 유력하다. 그러나 다양한 발전원을 이용하여 생산 및 수입 다변화를 모색하고 수요량에 대응하기 위해서는 추가적인 발전원을 통한 그린수소가 필요하다. 이와 같은 요구에 대응하기 위해 해양에너지(조류, 파력, 해수온도차발전 등)가 대안이 될 수 있다.
세계에너지기구(IRENA)1)의 보고서에 따르면 해양에너지 잠재량은 연간 해수온도차발전(OTEC)이 44,000 TWh, 파력발전 29,500 TWh, 조류발전 1,200 TWh 등으로 전망하고 있다. 해수온도차발전은 1881년 d’Arsonval에 의해 제안된 이후 미국, 한국, 일본, 프랑스 등의 나라에서 지속적으로 개발되고 있다. 현재 최고 용량의 실증은 국내의 선박해양플랜트연구소에서 수행한 2019년 1MW 바아지 탑재 폐쇄형 OTEC 플랜트이다. OTEC 플랜트는 해수를 감압·증발시킨 증기를 작동유체로 이용하는 개방형과 냉매를 이용하는 폐쇄형으로 구분할 수 있다. OTEC은 해수를 이용하기 때문에 수온에 영향을 받고, 파일럿 플랜트 규모에서 수온의 1℃ 변화 당 약 1.3 kW의 출력 변화가 보고된 바 있다.2) 폐쇄형 OTEC 플랜트의 가장 큰 특징은 개방형에 비해 상대적으로 높은 열효율이다.
고효율의 터빈 설계 및 해석은 높은 열효율 달성을 위해 중요한 요소이고, 냉매에 따른 OTEC 터빈의 특성도 고려되어야 한다. OTEC 터빈은 냉매를 작동유체로 사용하는 유기랭킨사이클(ORC) 터빈과 그 특성이 동일하지만 OTEC만의 특성(낮은 열원과 열침의 온도차이)으로 인해 Wet type의 냉매를 사용하게 될 경우 출구에서 액적이 발생하는 2상 유동을 갖는다. 따라서 OTEC용 터빈을 설계하고 전산유동해석(CFD) 프로그램으로 해석할 때에는 2상 유동 모델을 적용하는 것이 적절하다.
본 연구에서는 이러한 OTEC 터빈의 특성을 고려하여 CFD 해석 시 2상 유동 모델을 적용하여 그 결과를 분석하고자 한다. 특히 상변화(응축)가 매우 빠르게 일어나고 국부적으로 기체와 액체의 온도가 동일하다고 가정하는 평형 모델(Equilibrium model)이 아닌 Wilson point 이하에서는 기체와 액체의 온도가 다른 것을 모델링한 비평형 모델(Non-equilibrium model)과 단상(sing phase) 모델을 이용하여 해석을 수행하였다.
2. 2상 유동 모델
본 연구에서는 CFD 해석프로그램으로 Ansys의 CFX3)를 사용하였다. 본 프로그램에서는 2상 유동 해석을 위해 homogeneous 모델과 in-homogeneous 모델 모두 제공하고 있다. 전자는 두 개의 다른 상이 서로 상대적인 속도가 발생하지 않는 경우 사용하고, 후자는 두 개의 상의 속도가 다른 경우 사용한다. 예를 들면, 터빈에서 발생하는 액적의 경우, 그 질량이 충분히 작다고 볼 수 있으므로 발생하는 즉시 주변 기체와 같은 속도로 운동한다고 가정할 수 있다. 따라서 homogeneous 모델을 적용하는 것이 적절하다. In-homogeneous 모델을 적용해야 하는 대표적인 예는 정지 상태의 물이 담긴 용기를 끓여서 기체가 발생하는 경우이다.
평형 및 비평형 모델은 2상 유동이 발생할 때 두 개의 서로 다른 상의 온도에 따라서 구분된다. 상기 1절에서 서술한 바와 같이 비평형 모델은 두 상의 온도를 다르도록 모델링하여 해석하기 때문에 이 모델이 실제에 더 근접한 해석 모델로 볼 수 있다. 그 이유는 Chaplin R. A.4)에 따르면 다음과 같다. 기화 조건에서 자유 표면에 배치된 기포는 다른 기포에 비해 더 많은 분자에 둘러싸이게 되고(Fig. 1), 이로부터 벗어나기 위해서는 더 많은 에너지를 요구하게 된다. 따라서 일반적인 포화온도에 비해 순간적으로 기포의 온도가 높아지게 된다. 응축 조건에서는 이와 반대로 온도가 포화 온도보다 낮아지게 된다. 이러한 현상은 기포(혹은 물방울)의 성장이 끝나는 때까지 지속되고 이는 준평형 상태(Metastable equilibrium)로 지칭된다. 이 상태는 일시적으로 나타나고 팽창에 따른 압력 감소가 계속됨에 따라 열역학적 평형상태로 바뀌게 된다.
Effect of surface tension on small drops and bubbles4)
이와 같은 응축 과정을 과포화 팽창(supersaturated expansion) 과정이라 하고, 이 과정에서 발생하는 손실은 다음과 같이 추정할 수 있다. 등엔트로피 팽창이 천천히 그리고 충분히 일어날 경우에는 H-s 선도(Fig. 2)의 1지점에서 3s 지점으로 일어난다. 그러나 응축이 없는 준평형 상태로 팽창이 일어날 경우에는 2s에서 2로 프로세스가 진행되고 이 과정에서 엔트로피의 상승이 발생한다. 즉, 손실은 2s 지점과 2지점의 엔탈피 차이가 된다(참고로, 이 지점에서의 엔탈피와 같은 물성은 일반 물성 표가 아닌 metastable 표에서 구할 수 있다). 2지점에서 팽창이 더 일어나면 열역학적 평형상태에 도달하는 4지점으로 프로세스가 진행된다.
Comparison between normal and metastable expansion
3. 지배방정식
3.1 Equilibrium model
본 모델은 응축이 빠르게 일어나고, 두 개의 상의 온도를 동일하다고 가정한다. Equilibrium 모델의 지배방정식은 단상의 방정식과 매우 유사하다. 물성(ϕ)만 혼합물 형태로 대체되는데 이는 건도(χ)를 이용하여 계산한다. 물성치는 엔탈피, 엔트로피, 비열 등이다.
| (1) |
| (2) |
3.2 Non-equilibrium model
본 연구에서는 Ansys CFX의 The droplet condensation 모델을 이용하여 2상 유동을 모델링하였다. 본 모델은 정확도 향상을 위해 droplet 수와 지름을 계산한다. 또한 두 개의 서로 다른 상이 열역학적으로 평형을 이루는 것을 가정하지 않는 non-equilibrium 모델을 이용한다.
지배방정식은 continuous phase와 condensed phase를 포함하고, 둘의 운동 속도는 동일하다(homogeneous 모델 가정). 각 상의 연속방정식은 아래와 같다.
| (3) |
| (4) |
우측의 항들은 source 항으로서 첫 번째 항은 droplet growth 모델을 따르고, 두 번째 항은 nucleation 모델을 따른다. 첫 번째 항은 식 (5)와 같고, 접촉면(interfacial area)의 밀도와 droplet growth rate의 곱이다. dRd/dt는 Knudsen(Kn) 수의 영향을 고려하고 있는데, 이는 droplet의 크기가 연속성(continuum)이 인정되지 않는 영역에서부터 생성되기 때문이다. Jd는 핵형성(nucleation) 모델이고 단위 시간과 부피에 대해 생성되는 droplet의 수이다. 이때 A는 모델에 따른 상수, G*는 깁스 자유에너지 변화량, k는 Boltzman 상수, Tg는 과냉각 증기(supercooled vapor)의 온도이다.
| (5) |
| (6) |
| (7) |
운동 방정식은 아래 식 (8)과 같다. conservation 형태의 Navier-Stokes 방정식에 근거하여 부피 분율을 고려하고, momentum source(Sα), interfacial force(Mα), 상간의 질량 전달에 의한 운동량 전달(momentum transfer induced by interphase mass transfer) 항(Γαβ)이 포함된다.
| (8) |
continuous phase의 에너지방정식을 전엔탈피 형태로 나타내면 다음과 같다.
| (9) |
| (10) |
| (11) |
| (12) |
condensed phase에도 유사한 형태의 에너지 방정식이 적용될 수 있다. 다만 droplet의 크기에 따라 열전달 및 상변화 모델이 달리 적용되어야 하는데, 본 연구의 모델은 일반적으로 지름 1 um 이하의 액적의 발생이 예상되어 small droplet heat transfer 모델을 적용하였다. 열전달 source 항(SH)은 연속방정식에 이용되었던 Sd와 두 개의 상 사이의 열전달(Qd)에 의해 계산될 수 있다. Td는 포화온도에서 과냉각 온도와 droplet의 반지름 비율에 따라 정의된다.
4. CFD 해석 모델
4.1 Mesh and boundary conditions
CFD 해석에 사용된 터빈의 형상과 격자는 Fig. 3에서 보는 것과 같다. 무차원 벽 거리인 y+는 노즐과 로터의 블레이드 면과 교차하는 span 50%면인 선 위에서 각각 1.5와 2.1의 값을 가지며, 그 분포는 Fig. 4에서 볼 수 있다. 주요 경계 조건도 해당 그림에서 볼 수 있으며, 그 외의 사용된 경계 조건은 다음과 같다.
Generated mesh and boundary conditions
Boundary conditions for the turbine analysis
y+ distribution on blades at span 50%
Fig. 5에서 보는 것과 같이 실제 제작용 도면에서 터빈은 Volute가 별도로 없고, Diffuser의 형상은 periodic하므로 반경 방향을 따라 유동 조건이 균등한 가정으로 periodic 조건을 적용하여 해석 시간의 최소화를 위해 1개의 로터와 노즐에 대해 해석하였다. 난류 모델은 SST 난류 모델을 사용하였고, 작동 유체의 물성은 Peng & Robinson 상태방정식을 사용하여 계산하였다. 이 상태방정식은 Soave에 의해 개선된 Redlich-Kwong 상태방정식5)보다 포화액의 부피 예측에 우수한 성능을 나타내는 것으로 검증되었다.6) 압력과 속도 커플링은 Rhie and Chow에 의해 제안된 알고리즘이 수치해석 계산 프로그램에 내재적으로 사용되었다. 2상 유동 모델은 3절에서 소개된 두 모델을 각각 적용하였다.
Drawing of the 20 kW OTEC turbine and generator (Jinsol Turbo Inc.)
4.2 격자 의존성 및 이산화 오차
본 연구는 격자 의존성 및 이산화 오차를 추정하기 위해서 비평형 모델을 적용한 수치해석에 대해 Grid convergence index(GCI) method7)를 사용하였다. 이 방법은 이산화 오차를 정량화하여 퍼센트 단위의 불확도와 외삽법에 의한 값을 제시한다. 오차 추정 절차는 1) 격자 크기 정의, 2) 3개 이상의 수치해석 수행, 3) 해석결과를 이용한 GCI 법의 값(s, p, ε 등) 계산, 4) 외삽값(extrapolated value) 및 불확도 계산 순으로 진행된다. GCI와 관련된 구체적인 수식은 Celik 등8)의 연구에서 볼 수 있다. 본 연구의 GCI 법의 적용한 결과는 Table 3에서 보는 것과 같다. 격자 정제 비율(grid refinement factor, r)은 실험을 근거로 1.3 이상이 권장된다. 이 값은 약 1.4가 적용되었고, 각 해석 모델의 격자 수(Table 2)를 기반으로 계산되었다. ϕ는 각 모델 별 변수 값이고 ε은 각 값의 차이이다. ε21에 대한 ε32 비의 부호 함수(sign function)인 sf가 음수인 경우에는 진동 수렴(oscillatory convergence)을 의미한다. 국부 정확도(local order of accuracy)인 p는 0.7에서 5.2로 과냉각 온도에서 가장 낮았다. Mesh1과 Mesh2를 이용하여 외삽된 출력 값은 26,466이고 이를 이용한 상대 오차(eext21)는 1.79%이다. 이는 직접 오차 값(ea21)인 0.84% 대비 높게 나타난다. 이를 통해 조밀한 격자의 점근선과의 오차(numerical uncertainty)는 GCI 값을 통해 추정할 수 있고, 이는 각각 질량 유량 0.53%, 출력 3.35%, 과냉각 온도가 5.83%로 나타났다.
Number of nodes for different models
Grid convergence index and extrapolated values to estimate discretization error
5. 수치해석 결과
설계점(Table 1)에서 Mesh1 기준으로 단상(R32 기체) 모델, 2상 유동 모델(평형 및 비평형 모델)을 이용하여 수치해석을 수행하였다. 각 CFD 해석은 Table 4에서 보는 것과 같이 수렴하였고, 단상과 비평형 모델에서 더 낮은 잔차(Residual)로 수렴하였다. 주요 해석 결과는 Table 5에서 보는 것과 같다.
Convergence result by different phase models
Performance summary for different conditions at design point
실험 결과는 선박해양플랜트연구소에서 보유하고 있는 20 kW OTEC Pilot 플랜트(Fig. 6)를 이용하였다.9) 출력 밀도는 출력을 질량 유량으로 나눈 값으로 단위 질량당 출력을 비교하기 위해 사용하였다. 실험에 의한 결과는 설계 및 해석 결과보다 낮은 출력을 가지는데, 이는 일반적인 결과이며 설계 및 해석 조건에서 반영되지 않은 요소(디퓨저 손실, 표면 조도, 2차 유로에 의한 손실 등)에 의한 영향이다. CFD 해석 결과 간에 비교를 통해 실험의 출력 및 출력 밀도와 가장 유사한 결과를 나타내는 것은 비평형 모델인 것을 확인하였다. 그 이유는 작동 유체가 터빈을 통과하면서 팽창하는 과정에서 2상 유동이 예상되고, 이때 두 개의 상이 서로 열평형이 이루는 가정을 배제한 결과(상변화에 의한 엔트로피 상승의 반영)가 실제에 더 근접하기 때문이다. 2상 모델 중 평형 모델은 단상 모델보다 미소하게 더 높은 출력을 나타내는데, 모델의 안정성이 상대적으로 낮은 것에 기인하는 것으로 추정된다. 이는 Table 4에서와 같이 상대적으로 잔차가 큰 점과 해석 과정에서 격자 수가 5~10%만 변화해도 해석 중 발산하는 등 해석 안정성이 떨어졌기 때문으로 보인다.
KRISO`s 20 kW OTEC pilot plant installed in Goseong, Republic of Korea
정성적인 결과의 비교를 위해 압력 및 속도 분포를 비교하였다. 압력 분포는 해석 모델별 큰 차이가 없었고, 속도 분포는 Fig. 7에서 보는 것과 같이 로터 전방(leading edge)에서 국소적으로 저속부가 나타나는데, 이 영역이 비평형 모델에서 미소하게 크게 나타나는 것을 확인하였다.
Velocity contour for (a) single phase, (b) two phase model at span 50%
로터 입구 영역에서 발생하는 저속은 터빈 출력 저하의 원인이 될 수 있고, 이 영역이 확장된 것은 비평형 모델이 상대적으로 출력 밀도가 낮게 나타나는 것에 영향을 준 것으로 보인다. Fig. 8은 비평형 모델에 의한 과냉각 온도의 분포로 약 -1℃에서 6℃로 나타나는데, 이는 노즐에서 약 1도의 과열도를 가지고 있던 증기가 터빈을 통과하면서 팽창함에 따라 2상 유동이 발생하게 되고, 액체와 기체가 약 6℃ 정도의 온도 차이에 의한 것으로 보인다. 이와 같은 열역학적 불평형을 반영하여 모델링하는 것이 평형 모델과 비평형 모델의 차이이며, 과냉각 정도와 발생한 액적의 양이 두 모델 간의 출력 차이에 반영되는 것으로 보인다. 액적이 발생한 정도를 질량 분율을 보면 평형 모델은 약 2%를 차지(Fig. 9)하고, 비평형 모델은 0.1% 이하로 매우 낮게 나오는 것으로 확인하였다.
Supercooling distribution by non-equilibrium model at span 50%
Liquid mass fraction contour by equilibrium model at span 50%
Baumann10)에 따르면 습도(1-χ)가 1% 증가할 때마다 건조 조건에서 등엔트로피 터빈 효율은 1%가 감소하는 것으로 알려져 있다. 그러나 이후 연구에서는 Baumann factor가 수증기에 대해 0.4~2.5로 변화하는 것으로 보고되고 있다.11)
6. 결 론
본 연구에서는 폐쇄형 OTEC 플랜트의 터빈에서 발생할 수 있는 2상 유동을 CFD 해석을 수행하고 단상 및 2개의 2상 모델, 실험값을 이용하여 정량적, 정성적으로 비교 분석하였다. 해석 결과를 정리하면 다음과 같다.
1) 20 kW OTEC의 터빈(노즐, 로터)를 해석하기 위해 수치해석 모델을 구축하고 단상, 2상(평형 및 비평형 모델)을 적용하였다.
2) 비평형 모델의 격자의존성 및 오차를 분석하고자 GCI 법을 이용하였고, 4가지 변수(출력, 질량유량, 증기 질량 분율, 과냉각 온도)에 대해 진동이 없는 것을 확인하고 상대오차 및 GCI 값을 제시하였다.
3) 단상, 2상 유동 모델은 모두 수렴하였고, 가장 해석 시간을 크게 요구한 평형 모델에서 잔차 또한 가장 높게 나타났다.
4) 설계점에서 CFD 해석 결과와 실험 결과를 비교하였고, 그 결과 비평형 모델의 출력 밀도가 7.98 kJ/kg으로 가장 실험(6.46 kJ/kg)과 근접한 결과를 나타내었다.
5) 단상 모델과 2상(비평형) 모델의 속도 분포를 비교하였고, 미소하게 로터 전방에서 속도가 다르게 나타나는 것을 확인하였다.
6) 2상 모델(평형 및 비평형 모델) 간의 질량 분율을 비교하였고, 약 2%로 평형 모델에서 질량 분율이 높게 나타나는 것을 확인하였다.
― 기 호 설 명 ―
| A : | 핵형성(nucleation) 모델 상수 |
| c : | 음속 |
| e : | 오차 |
| G : | 깁스 자유에너지 변화량 |
| GCI : | Grid convergence index |
| H : | 엔탈피 |
| J : | 핵형성 모델 |
| k : | Boltzman 상수 |
| Kn : | Knudsen 수 |
| nd : | 액적 수(number of droplet) |
| m : | 질량 |
| M : | 표면력 |
| p : | 국부 정확도 |
| P : | 압력 |
| Q : | 단위면적당 열전달량 |
| r : | 부피 분율 |
| R : | 반지름 |
| s : | 엔트로피 |
| S : | Source 항 |
| sf : | 부호 함수(sign function) |
| t : | 시간 |
| T : | 온도 |
| U : | 속도 벡터(u, v, w) |
| y+ : | 무차원 벽 거리 |
그리스 문자
| Γ : | 운동량 전달량 |
| ϵ : | 차이 값 |
| ϕ : | 물성, 매개 변수 |
| χ : | 건도 |
| ρ : | 밀도 |
상첨자
| 1,2,3 : | Mesh number |
| * : | critical property |
하첨자
| 1,2,3 : | Mesh number |
| α : | α phase |
| β : | β phase |
| αβ : | from α phase to β phase |
| a : | absolute value |
| c : | continuous phase |
| d : | condensed phase, droplet property |
| ext : | extrapolated value |
| fine : | fine grid value |
| H : | heat transfer property |
| l : | liquid phase |
| m : | mass property |
| M : | momentum property |
| sat : | saturated condition |
| sc : | supercooled level in gas phase |
| t : | total property |
| v : | vapor phase |
Acknowledgments
본 논문은 해양수산부의 국가 R&D 사업인 “선박폐열을 활용한 발전 시스템 및 RCI 기술 실증(1/6) (PMS5380)”에 의해 수행되었습니다.
Author contributions
J. B. Seo; Conceptualization, Writing-original draft, Investigation, Methodology. H. S. Lee; Writing-review & editing. H. J. Kim; Supervision.
References
- International Renewable Energy Agency, 2020, "Innovation Outlook ocean energy technologies", International Renewable Energy Agency.
-
J. H. Moon, H. S. Lee, S. T. Lim, J. B. Seo and H. J. Kim, 2020, "Plant Performance Evaluation Experiment to Expand the Applicability of Ocean Thermal Energy Conversion", Journal of Power system Engineering, Vol. 24, No. 5, pp. 78-85.
[https://doi.org/10.9726/kspse.2020.24.5.078]
- C. F. X. Ansys, 2009, ANSYS CFX-solver theory guide. Ansys CFX Release, 15317, pp. 724-746.
- R. A. Chaplin, 2009, Thermal power plants. Encyclopedia of Life Support Systems.
-
G. Soave, 1972, "Equilibrium constants from a modified Redlich-Kwong equation of state", Chemical engineering science, Vol. 27, Issue 6, pp. 1197-1203.
[https://doi.org/10.1016/0009-2509(72)80096-4]
-
C. T. Lin and T. E. Daubert, 1980, "Estimation of partial molar volume and fugacity coefficient of components in mixtures from the soave and Peng-Robinson equations of state", Industrial & Engineering Chemistry Process Design and Development, Vol. 19, No. 1, pp. 51-59.
[https://doi.org/10.1021/i260073a009]
-
P. J. Roache, 1994, "Perspective: a method for uniform reporting of grid refinement studies", Journal of Fluids Engineering, Vol. 116, Issue 3, pp. 405-413.
[https://doi.org/10.1115/1.2910291]
-
I. B. Celik, U. Ghia, P. J. Roache and C. J. Freitas, 2008, "Procedure for estimation and reporting of uncertainty due to discretization in CFD applications", Journal of fluids Engineering, Vol. 130, Issue 7, 078001.
[https://doi.org/10.1115/1.2960953]
- H. S. Lee, Y. G. Jung, H. J. Kim and B. S. Choi, 2014, "Design and Experiment of the 20kW OTEC", The korean Association of ocean science and technology societies, pp. 2448-2452.
-
K. Baumann, 1921, "Some recent developments in large steam turbine practice", Journal of the institution of electrical engineers, Vol. 59, Issue 302, pp. 565-623.
[https://doi.org/10.1049/jiee-1.1921.0040]
- M. J. Moore and C. H. Sieverding, 1976, "Two-phase steam flow in turbines and separators: theory, instrumentation, engineering", United States.