[ Original Article ]

Journal of the Korean Society for Marine Environment and Energy - Vol. 20, No. 1, pp.37-44

ISSN: 2288-0089 (Print) 2288-081X (Online)

Print publication date 25 Feb 2017

Received 31 Jan 2017 Revised 06 Feb 2017 Accepted 13 Feb 2017

DOI: https://doi.org/10.7846/JKOSMEE.2017.20.1.37

다상유동형 입자법을 이용한 Rayleigh-Taylor 불안정성의 수치해석

김경성¹

; 구본국²

; 김무현³

; 박종천⁴^{, †}

; 최한석⁵

; 조용진⁶

1동명대학교 조선해양공학과
2창원대학교 조선해양공학과
3Texas A&M University Civil Engineering
4부산대학교 조선해양공학과
5포항공과대학교 엔지니어링대학원
6동의대학교 조선해양공학과

Numerical Study on Rayleigh-Taylor Instability Using a Multiphase Moving Particle Simulation Method

Kyung Sung Kim¹

; Bonguk Koo²

; Moo-Hyun Kim³

; Jong-Chun Park⁴^{, †}

; Han-Suk Choi⁵

; Yong-Jin Cho⁶

1Department Naval Architect & Ocean Engineering, Tongmyong University, Busan 48520, Korea
2Department Naval Architect & Marine Engineering, Changwon National University, Changwon, 51140, Korea
3Department Civil Engineering, Texas A&M University, College Station 77843, TX, USA
4Department Naval Architect & Ocean Engineering, Pusan National University, Busan 46241, Korea
5Graduate School of Engineering Mastership, POSTECH, Pohang 37673, Korea
6Department Naval Architect & Ocean Engineering, Dong-Eui University, Busan 47340, Korea

Correspondence to: ^† jcpark@pnu.edu

초록

하나의 시스템 내에 2개 이상의 상이 다른 유체가 존재할 시에는 다상유동에 의한 복장성이 존재하며, 이는 해석의 어려움이 따른다. 두 개 이상의 상이 다른 다상유동은 유동 및 경계면에 영향을 끼치지 때문에, 불안정성과 같은 비선형 유동이 나타나게 된다. 여러 종류의 불안정성 중 레일리히-테일러 불안정성은 대표적인 예로 알려져 있다. 본 연구에서는 밀도차가 레일리히-테일러 불안정성에 미치는 영향을 조사하기 위해 다양한 Atwood 수를 선정하였으며, 초기 경계면 형상 역시 다양한 형태를 설정하고 시뮬레이션 하였다. 본 연구에서 사용된 입자법인 MPS(Moving particle simulation)은 이러한 다상유동에서 널리 쓰이지는 않았으나, 다상유동을 위한 입자간 상호 연성 모델인 자가-부력 항, 표면 장력 항과 경계면 경계 조건 항을 추가로 사용하여 수치해석이 가능하게 하였다. 본 연구에서 새로이 개발된 다상유동형 입자법을 이용하여 고려된 경우들에 대해 수치해석을 수행하였으며, 각각의 결과들을 비교 분석하였다. 또한 레일리히-테일러 불안정성에 기인한 유동의 속도를 측정하여 포텐셜 기반의 이론값과의 비교를 통해 경향성이 일치함을 알 수 있었다. 이론값과의 크기의 차는 포텐셜 기반의 이론값에서는 고려가 힘든 비선형성에 기인한다고 사료된다.

Abstract

Complexity of multiphase flows due to existence of more than two interface including free-surface in one system, cannot be simulated easily. Since more than two fluids affect to flows and disturb interface, nonlinearities such as instabilities can be appeared. Among several instabilities on multiphase flows, one of representative is Rayleigh-taylor instability. In order to examine in importance of density disparity, several cases with numerous Atwood number are set. Moreover, investigation of influence on initial disturbance were also considered. Moving particle simulation (MPS) method, which was employed in this paper, was not widely used for multiphase problem. In this study, by adding new particle interaction models such as self-buoyance correction, surface tension, and boundary condition at interface models, MPS were developed having more strength of physics and robust. By applying newly developed multiphase MPS, considered cases are performed and compared each other. Additionally, though existence of disagreement of magnitude of rising velocity between theoretical values from linear potential theory and that of numerical simulation, agreement of tendency can be proved of similarity of result. the discordance of magnitude can be explained due to non-linear effects on numerical simulation which was not considered in theoretical result.

Keywords:

Moving particle simulation, Multiphase flow, Rayleigh-Taylor Instability

키워드:

입자법, 다상유동, 레일리히-테일러 불안정성

1. 서 론

일반적으로 해양공학에서 저장조(reservoirs)로부터 프로세싱 플랜트로 수송되는 오일과 가스 흐름(Oil and gas flow)은 다상유동이며, 특히 관내에서의 다상류에 관한 예측과 컨트롤은 안전성과 경제성 측면에서 매우 중요하다.

이러한 다상유동을 해석하기 위한 유체역학은 크게 이론적 방식, 실험적 방식, 그리고 수치 해석적 방식으로 크게 분류될 수 있다. 이 중 수치 해석적 방식은 다시 3가지 방식으로 나눌 수 있으며, 이는 격자 방식의 오일러 접근법(Eulerian approach), 입자 방식의 라그란지 접근법(Lagrangian approach), 그리고 그 둘을 혼합한 Hybrid 방식으로 분류될 수 있다. 첫 번째 접근법인 오일러 접근법은 정지 또는 이동 격자계를 사용하며, 각 상(Phase)의 경계면(interface)을 표현하는 스칼라 함수의 수송방정식을 시간 적분에 의해 풀어가는 방법이다. 이러한 접근법은 비교적 간단한 알고리즘으로 계산시간이 빠르다는 장점이 있지만, 경계면의 수치적 번짐(Numericallydiffused)으로 인해 높은 정확도를 요구하는 경계면 탐색에 어려움이 있으며, 특히 경계면의 붕괴 혹은 대변형 문제에 있어 제약이 따른다. 기존에 주로 사용되는 방법으로서는 VOF(Volume-of-Fluid)법 (Chen, 2011)이나 Level-Set법(Sussman et al., 1994) 등이 있다. 두 번째로 Hybrid 방식은 오일러 방식의 격자를 사용하면서, 입자나 무질량의 마커를 이동시켜 경계면을 후술하는 라그란지 접근법으로 추적하는 방식이다. 이 방식은 오일러 접근법과 라그란지 접근법을 효과적으로 혼합하여 비교적 빠른 계산 시간과 높은 정확도를 동시에 기대할 수 있지만, 그 접근법의 연성이 복잡하고 잘못된 연성으로 인해 수치 불안정성이 발생할 가능성이 높다. 대표적인 방법으로 PIC(Particle-In-Cell)법(Harlow, 1955)과 MAC(Marker-and-Cell)법 (Harlow & Welch, 1965)을 들 수 있다. 마지막으로 라그란지 접근법은 오일러 방식의 격자 대신 계산 영역의 모든 곳을 자유롭게 이동할 수 있는 입자를 사용하는 방법이다. 제안된 방법으로는 SPH(Smoothed Particle Hydrodynamics)법 (Monaghan, 1994), MPS(Moving Particle Semi-implicit)법 (Koshizuka and Oka, 1996), DPD(Dissipative Particle Dynamics)법 (Moeendarbary et al., 2009) 등이 있다.

입자법의 경우, 경계의 대 변형과 붕괴에 높은 정확도를 나타내며, 특히 각각의 입자가 고유의 인식 번호와 함께 물리량을 함께 갖고 이동하는 특성은 다상 유동에서 경계면 탐색 및 이동에 유리함을 보인다. 입자법 중 하나인 SPH법은 초기 천문학 분야에서 압축성 유체의 해석에 적용하기 위해 개발되었으며 이 후 고체 및 비압축성 유체의 유동해석에 그 적용범위를 넓히고 있다. SPH법은 유체를 이산화된 입자들 사이의 상호작용을 커널(Kernel)이라는 가중치 함수를 도입하여 물리량의 적분과 편미분 연산자를 수학적인 모델을 유도하여 풀게 된다. 또 다른 입자법 중 하나인 MPS법은 연속체의 지배방정식에 포함된 편미분 연산자들의 이산적인 계산을 위해 차분 근사와 분자동력학의 통계적 개념을 바탕으로 하여 입자간 상호작용으로 치환하며, 개발 초기부터 자유표면을 갖는 비압축성 유체해석을 위해 개발되었다. MPS법은 입자를 사용한다는 것에서 SPH법과 유사성이 있으나, MPS법에서는 가중치 함수인 구배와의 결합이 배제된 평균 가중치 과정을 독자적으로 사용하는 편미분 방정식을 사용한다는 차이점이 있다. 이러한 차이는 압력 계산 방식에 큰 차이를 가져온다. 일반적으로 압력 계산에서 물리우주론에서 사용하는 상태방정식(Equation of state)을 양적으로 풀어내는 SPH법은 계산 속도에서 우위를 점하고 있다. 반면 비압축성 자유표면 유동해석에 사용할 목적으로 개발된 MPS는 포아송 압력 방정식(Poisson Pressure Equation)으로 압력을 매 시각 반복 계산하므로 연산 속도는 상대적으로 느리지만 SPH법에 비해 비압축성 유동 해석에 비교적 높은 정확도를 가진다.

특히, 상기 언급한 입자법 중 MPS법은 Koshizuka et al.(1996)에 의해 처음 제안되었으며, 이 후 Lee et al.(2010)에 의해 압력 안정화와 자유 표면 탐색의 정확도를 높였다. 또한 기체-액체 해석을 위해 Shirakawa et al.(2001)에서는 부력 수정 모델, Nomura et al.(2001)은 표면 장력 모델을 추가 도입하여 다상 유동에서 발생하는 급격한 물리량의 변화에 따른 수치 불안정성을 경감시키는 노력을 하였다. 이 후 Kim et al.(2014)은 제안된 부력 수정 모델과 표면장력 모델을 물리적 관점에서의 강화를 통해 이론적 견고함을 높혔다. 입자법을 이용한 다상 유동에 대한 연구는 Shakibaeinia et al.(2012)에서 Kelvin Helmholtz 불안정성 해석을 시뮬레이션 하였고, Khayyer et al.(2013)은 밀도의 완충지역 설정으로 밀도차이가 큰 기체-유체간의 동적 해석을 수행하였다. 또한 Kim et al.(2014)에서는 입자법의 활용범위를 다상유동에까지 확장시켜 3상을 가진 슬로싱 문제에 적용하였다.

일반적으로 단일 유체의 유동에 비해 다수의 상이 다른 유체의 유동은 복잡함을 띤다. 특히 유체-유체간의 경계면에서 급격한 물리량의 변화로 인해 수치적 불안정성에 의한 경계면의 붕괴 및 대변형 등에 대한 해석의 어려움이 존재한다. 본 연구에서는 해양공학 분야에서 필수적으로 요구되어지는 다상 유동의 해석을 위해 입자법을 채용하였으며, 이에 따라 압력 계산에 정확도를 보이는 MPS 계열의 방법을 채용한다. Kim et al.(2014)에서 제안한 다상 유동형 MPS법을 사용하여 다상유동시 발생하는 불안정성 중 하나인 Rayleigh-Taylor 불안정성에 대해 다양한 밀도 조합을 이용하여 수치해석을 수행하였으며, 다양한 초기 경계면의 형상에 따른 불안정성으로부터 기인한 유동의 발달 양상에 대해 연구를 수행하였다.

2. 수치해석법

2.1 단상 유동형 MPS 법

일반적인 MPS법은 단상(Single phase)의 유체 해석에 사용되어왔다. MPS법의 지배방정식은 연속방정식과 Navier-Stokes 방정식이다.

D ρ D t = 0

(1)

D u ⃗ D t = - 1 ρ ∇ p + v ∇ 2 u ⃗ + σ k n ⃗ + F ⃗

(2)

여기서 ρ는 밀도, t는 시간, $u ⃗$ 는 속도, p는 압력, v는 동점성계수, $F ⃗$ 는 외력, ∇는 구배, σ는 표면 장력 계수, k는 포면장력시의 작용곡률, 그리고 $n ⃗$ 는 법선 벡터를 나타낸다. 격자기반의 전산유체역학에서는 연속방정식 (1)이 속도의 발산(∇· $u ⃗$ )로 표현되는 것에 반해 완전 라그란지 접근법인 입자기반의 전산유체역학은 밀도의 시간에 대한 전미분의 형태(Dρ/DT)로 표현되므로, 이는 질량보전의 법칙이 완벽하게 만족되는 특징이 있다. Navier-Stokes 방정식 (2)도 격자기반의 방식과는 차이가 있다. 격자 기반의 Navier-Stokes 방정식은 좌변이 이산항을 포함하는 2개의 항[ $∂ u ⃗ / ∂ t + u ⃗ ∂ u ⃗ / ∂ x$ ]으로 표현되는 것에 반해, 입자 기반의 방식에서는 좌항이 속도에 대한 시간의 전미분 형식 [ $D u ⃗ / D t$ ]이 되며, 이로 인해 이산항에서 발생할 수 있는 수치 오류가 없다는 장점을 지닌다. Navier-Stokes 방정식의 우변은 각각 압력 구배항, 점성항, 외력항으로 크게 구분이 되며 본 연구에서는 표면장력을 고려한 표면장력 항을 우변의 3번째 항에 추가하여 시뮬레이션 하였다.

입자 기반의 비압축성 점성 유동의 구현을 위해 지배방정식의 모든 항들은 MPS법에 적합한 형태로 변환되어야 하며, 이를 위해 입자간 상호 모델을 이용한 형태로 변환하였다. 표면장력항을 제외한 모든 항들은 Lee et al.(2010)에서 제시한 개량된 PNU-MPS법을 따랐으며, 표면장력 모델은 Kim et al.(2014)에서 제안된 형태를 사용하였다.

2.2 다상 유동형 MPS법

앞서 서술된 입자법은 단일 유동 해석에는 높은 정확도를 보인다. 그러나 입자가 각각의 물리량을 가지고 이동을 한다는 특징에 의해 입자와 입자사이의 인터페이스에서에 급격한 물리량 변화가 있는 다상 유동에서는 매우 제한된 범위 내에서만 구현이 가능하였다. 특히 전방향 확산 형태인 MPS법의 가중치 함수로 인해 밀도 차이에 의한 부력의 방향 및 크기가 잘 못 계산되어지는 문제가 발생할 수 있다. 이를 해결하기 위해 Shirakawa et al.(2001)에서는 부력 항을 추가하였고, Kim et al.(2014)에서는 물리적 의미를 강화한 형태의 자가 부력 수정 항을 도입하여 발생가능한 오류를 제거 혹은 완화하였다. 단일 유체(특히 자유표면을 갖는 유동) 문제 해석의 경우, 기체는 빈 공간으로 간주하여 기체 입자를 배치하지 않고 시뮬레이션을 하는 것이 MPS법에서는 일반적이다. 이 경우 자유 표면에서의 표면장력은 그 비중이 낮아 계산 결과에 영향을 미치지 않는 것에 반해, 다상 유동의 경우는 유체-유체간의 경계면에서는 밀도가 기체에 비해 상대적으로 높은 유체들이 배치되어 있기 때문에 표면장력 효과를 무시할 수 없다. 그 결과 표면장력 모델이 Nomura et al.(2001)에 의해 도입되었고, Kim et al.(2014)에서는 법선 벡터 탐색 방식을 수정한 표면장력 모델을 이용하였다. 본 연구에서는 앞서 언급된 개선한 자가 부력 수정 모델과 표면장력 모델을 사용하였으며, 경계입자 탐색법과 유체-유체간 경계 조건 모델을 추가하여 다상 유동 해석에 대한 MPS법의 수치 정확도를 높였다.

본 연구에서 사용된 표면장력모델은 Nomura et al.(2001)의 방식에서 법선벡터 탐색법을 저자들에 의해 수정된 방식을 사용하였다. 표면장력 모델을 위해 2개의 새로운 입자수 밀도, n_i^st1과 n_i^st2를 도입하였으며 그 식은 아래와 같다.

n i s t 1 = ∑ j ≠ i w s t 1 | r ⃗ j - r ⃗ i w h e r e w s t 1 = 1 0 < r ≤ r e s t 0 r e s t < r

n i s t 2 = ∑ j ≠ i w s t 2 | r ⃗ j - r ⃗ i w h e r e w s t 2 = 1 0 < r ≤ r e s t 0 r e s t < r

여기서 w^st1과 w^st2는 가중치 함수, r은 입자간 거리, r_e^st는 유효거리이며, 수치적 방식의 검증을 통해 유효거리는 1.3으로 설정하였다. 새로우 구한 입자수밀도를 이용하여 표면장력을 위한 곡률은 다음의 식으로 구할 수 있다.

k = 1 R = 2 r e s t cos ⁡ 0.5 π n i s t 2 n 0 s t 1

여기서 n₀^st1은 초기 배치에서의 입자수밀도이다. 곡률을 구한 뒤 법선벡터는 다음의 식으로 구할 수 있다.

a ⃗ i = n ⃗ x ∑ j ≠ i x n w s t 1 | r ⃗ j - r ⃗ i + n ⃗ y ∑ j ≠ i y n w s t 1 | r ⃗ j - r ⃗ i n ⃗ i = a ⃗ i a ⃗ i

여기서 x_n과 y_n은 중심입자와 주변입자간의 수평 및 수직거리이다.

3. 수치해석 조건 및 결과

다상 유동과 같이 두 개 이상의 상이 다른 유체가 혼합되어 하나의 시스템에 존재할 때, 복잡한 유체의 운동이 발생하며, 단일의 현상이 아닌 다양한 불안정성에 기인한 경계면 붕괴 및 유체간의 혼합이 일어나기 때문에 이를 해석하기가 쉽지 않다. 이를 위해 반복적으로 발생하는 현상을 관찰하여, 반복성의 확인 후 각각의 불안정성을 명명할 수 있으며, 대표적으로는 시스템의 붕괴 시 경계면의 형상에 따라 나뉠 수 있다. 경계면이 버섯 형태의 유동을 만드는 경우는 Rayleigh-Taylor Instability(RTI)라 불리며, 이는 상이 다른 유체의 내부로 다른 유체가 침투하며 발생을 한다(Sharp, 1984). 또 다른 대표적인 불안정성은 Kelvin-Helmholtz Instability(KHI)가 있는데, 이 경우는 평행하게 배치된 유체가 서로 반대방향으로 운동을 하거나 경계면에서 유체간의 상대 속도가 일정 이상 높을 시 경계면에서 쇄파(Breaking Wave)형태의 유동을 생성할 경우를 일컫는다.

본 연구에서는 다상 유동에서 발생할 수 있는 불안정성 형태 중 RTI에 대한 연구를 수행하였다.

RTI의 경우는 일반적으로 수직관에서 2개 이상의 상이 다른 유체를 배치하여 구현하는 것이 일반적으로 많이 쓰이는 모델이다. 이 경우 수직관의 상부에는 무거운 유체를 하부에는 가벼운 유체를 배치하여 중력에 의한 자발적인 유동이 생성되도록 하여, 추가적인 외력이 존재하지 않을 때를 가정한다. 본 연구에서도 같은 방식의 일반적인 방식의 RTI 생성 모델을 설정하였으며, 기하학 집중 부를 형성하여 RTI 시뮬레이션의 반복성을 확보하고, 원활한 유동의 생성을 가능하도록 유체-유체간의 경계면을 설정하였다. 시뮬레이션에서 사용된 초기 경계면의 형상은 Fig. 1에서 보이는 것과 같이 3가지 형태의 경계면을 설정하였으며, 그 경계면은 다시 3개의 각기 다른 진폭으로 형성하여, 총 9개의 경계면의 형상에 대해 시뮬레이션 하였다. 또한 본 시뮬레이션은 중력에 의한 자발적인 유동이 생성되도록 설정하였기 때문에, 유체의 밀도에 대해서도 3가지 밀도비의 조합을 선택하였으며, 이는 밀도비를 나타내는 Atwood수인 AT = (ρ₂−ρ₁)/(ρ₂+ρ₁)의 설정에 상응하도록 선택하였다. 본 연구에서 사용된 시뮬레이션 조건은 Table 1에 나타나 있으며, 3가지의 밀도비 변수에 9가지의 초기조건을 갖는 경우의 조합을 고려하여 총 27개의 시뮬레이션을 수행하였다. 각 조합을 읽는 방식은 첫 번째 항은 밀도비인 Atwood수, 두 번째 항은 경계면의 형태, 세 번째 항은 경계면의 진폭을 나타낸다. 예를 들어, Case 1-II-(c) 의 경우는 Atwood수가 0.2이며 경계면의 형상은 η = Acos(4x−0.5)π을 따르며 그 경계면의 진폭(A)은 0.075 m임을 나타낸다.

Fig. 1.

Schematic model of Rayleigh-Taylor Instability with various interface shape.

Fig. 2는 Case (1, 2, 3)-I-((a), (b), (c)), 즉 경계면의 형태가 I타입인 경우에 3가지 Atwood수를 적용한 유체를 이용하였으며, 사용된 경계면의 진폭은 0.025 m, 0.050 m, 그리고 0.075 m를 각각 사용한 시뮬레이션의 결과를 임의의 시간인 1.0초, 2.0초, 3.0초 그리고 4.0초에 대해 그림으로 나타낸 것이다. I 타입의 경계면은 하부유동이 상부로 향하는 하나의 기하학 집중을 가지는 형태이다. 예측되는 결과는 기하학 집중부로부터 상승 유동이 발달하며, 상승유동이 일어난 후 빈 공간을 상부의 유동이 채워지며 급속한 하강유동을 만든다. 달리 말하면, 하부에 생긴 빈 공간을 채우기 위해 내려간 상부의 유체역시 빈 공간을 형성하게 되고, 이로 인해 상승유동이 상부의 유동의 방해를 적게 받으며 원활히 상승할 수 있다고 판단된다. Fig. 2에서 보이는 시뮬레이션 결과는 앞서 예측한 결과를 잘 따르고 있는 것으로 보인다. 그러나 Case 1의 경우는 2와 3의 경우와는 다른 양상을 보인다.

Fig. 2.

Snapshots of Simulation results from Case [1,2,3]-I-[(a),(b),(c)].

Case 1-I-(a)의 경우는 기하학 집중부에서 상승 유동이 발생하지 못하고 그 주변부인 기하학 특이점에서 상승 유동이 발생하는 것으로 보인다. 이는 에너지의 집중을 위해 기하학 집중부를 인위적으로 형성하였으나, 두 유체간의 밀도차가 크지 않아, 낮은 부력에 의한 낮은 상승 속도가 생성되었고, 이는 상부의 무거운 유체를 관통하며 상승할 수 있는 상승 유동을 만들 수 없었다고 보인다. 그로인해 경계면에서 기하학 특이점에서의 경계면 붕괴가 먼저 일어나는 결과를 만들었고, 붕괴가 일어난 특이점에서 상승 유동이 발생하였다. Case 1-I-(b), (c)의 경우는 상대적으로 높은 진폭을 가진 경계면에서 에너지 집중이 Case 1-I-(a)에 비해 높게 일어났기 때문에 기하학 집중부에서의 상승 유동이 발행하였다. 그러나 Case 2와 3의 다른 경우들과 비교를 하였을 때, Case 2와 3은 RTI의 전형적인 형태인 버섯형태의 유동을 만들며 상승하였으며 경계면의 진폭이 RTI 유동 발달의 시간에만 변화를 준 것에 반해, Case 1-I-(b)와 (c)는 상승 유동만 발생하고 버섯형태의 유동을 유발하는 회전력을 가진 유동을 발생시키지 못하였으며 진폭이 낮은 Case 1-I-(a)는 특이점에만 상승유동이 발생하는 것을 보였다. 이는 역시 Atwood수에 따른 밀도 차에 의해 발생하였다고 볼 수 있는데, Case 2와 3의 경우 충분히 확보된 밀도 차에 의해 생성된 부력으로 만들어진 상승 유동이 종단속도에 도달한 뒤에도 잉여 에너지가 남아 횡 방향으로 에너지를 발산하여 이 후 하강 유동을 가질 수 있다고 판단된다. 그러나 작은 밀도 차에 의해 Case 1에서는 대류로 인식될 수 있는 횡방향 유동과 하강 유동이 발견되지 않았다. 이후 계산 시간이 증가함에 따라 충분히 생성된 에너지에 의해 다른 유체를 관통하는 유동을 만들 수 있었고 이는 RTI의 전형적인 형태인 버섯형태의 유동을 생성하였다. 그러나 Case 2와 3에서 보이는 것과는 다르게 경계면의 형상이 매끄럽지 않다. Case 1의 경우는 상대적인 에너지의 차가 크지 않아 경계면에서 두 유체가 지속적으로 충돌을 일으키며 경계면을 흐트러트렸기 때문으로 판단된다.

Table 1.

Simulation conditions and properties of fluids

수행 된 RTI의 수치해석 결과 중 진폭이 0.025m 인 Case1-[I,II,III]-(a)의 유동 발달을 측정하여 이론값과 비교하였다. 선형이론에 의하면 불안정성에 기인한 경계의 발달은 지수적으로 증가하며 경계면의 형상을 식으로 나타내면 아래의 식으로 표현이 가능하다.

η t = A e ω t

(3)

여기서 η는 경계의 위치 함수이며, t는 시간, ω는 발달율을 나타내며 다음과 같이 해석적으로 구할 수 있다(Mikaelian, 1996).

ω 2 + 2 v k 2 ω - A T g k 1 - k 2 k c 2 = 0

(4)

여기서 k는 파수, k_c는 참고 파수, 그리고 v는 동점성 계수를 나타낸다. 위 식에서 나타나 있듯이 이론값은 점성효과를 고려하였다. 또한 참고 파수 k_c는 아래의 수식으로 구할 수 있다.

k c = ρ 2 - ρ 1 g T

(5)

단, T는 표면장력을 나타낸다.

식 (4)의 해를 식 (3)에 대입하여 구하여진 경계의 위치와 시간은, 아래에 기술된 방식으로 무차원화 할 수 있다.

τ R T * = t g λ, η R T * = η λ

(6)

여기서, λ는 초기 경계면 형상의 파장을 나타낸다.

Fig. 3은 RTI 유동 발달에 대해 수치해적 결과와 이론값의 무차원화 결과를 비교하여 나타낸다. Atwood수의 변화에 따라 유동 발달 속도는 변화하며, 그 변화량에 있어 이론과 수치해석 결과가 동일한 경향성을 띈다. Atwood수가 0.2와 0.4사이의 변화정도에 비해 0.4와 0.6 사이의 변화정도가 이론값과 수치해석 결과 모두에서 상대적으로 크게 나타나며, 전부 지수적 증가 형태를 띠고 있다. 경향성의 일치와는 달리 RTI의 크기는 일치하지 않는 것이 나타나는데, 이는 타 수치해석 결과에서도 동일한 양상을 나타낸다(Jeong et al., 2013). 크기의 불일치는 포텐셜 이론을 기반으로 한 이론값에서는 고려가 힘든 비선형성에 기인한다고 사료된다.

Fig. 3.

Comparison of RTI evolution between numerical and theoretical values.

Fig. 4에서는 경계면의 형태가 II인 경우에서 각기 다른 Atwood 수과 경계면의 진폭을 적용한 시뮬레이션, 즉 Case [1,2,3]-[II]-[(a),(b),(c)]의 결과들을 나타낸 것이다. 선택된 계산 시간은 이전과 동일한 1.0초, 2.0초, 3.0초 그리고 4.0초이다. 이 경우는 앞서 비교되었던 초기 변형이 I인 경우와는 다른 형태를 가진다. 초기 변형 II 형태는 하부 유체가 2개의 상부로 향하는 기하학 집중부를 가진다. 다시 말하면, 이는 상부 유체가 하부로 향하는 1개의 기하학 집중을 가진다고 볼 수 있다. 그로인해 상승 유동으로 인한 RTI 불안정성에 기인한 유동보다는 하강 유동에 의한 유동발달이 두드러진다. 이에 따라 하부로 향하는 유동에서 RTI 유동이 뚜렷하게 나타나며, 이는 Atwood수가 낮은 Case 1에서도 보여진다. 그러나 낮은 Atwood수에 기인한 낮은 부력으로 인해 초기에는 경계면의 출렁거림이 나타났고, 출렁거림에 의해 발생한 취약점으로부터 하강 유동이 발달함을 보였다. Case 2와 3의 비교에서는 이전의 경우인 Case [2,3]-[I]-[(a),(b),(c)]의 경우들과 같이, 밀도 차에 의한 부력차이에 기인한 유동 발달 시간이 다른 것 외에 크게 다른 점을 나타나지 않았다. 앞서 보인 I타입의 결과와 II 타입의 경계 형태와의 가장 큰 차이는 Case 1에서 보였다. 이전의 경우는 낮은 부력으로 인해 상부 유체와 하부 유체간의 내부 에너지의 차가 크지 않아 경계면이 혼잡한 모습을 보였으나, II 타입에서는 2개의 상승 유동으로 인하여 하부에 생긴 공간으로 상부 유체가 밀려들어가는 형태를 띄었으며, 이는 관통력이 낮은 경우지만 하강 유동에서 회전 유동이 일어날 수 있는 충분한 하강 유동을 발달시켰다. 그로인한 에너지 집중 및 분산이 원활히 일어나서 경계면에서 유체간의 부딪힘이 적어 깨끗한 경계면이 형성될 수 있다고 판단된다.

Fig. 4.

Snapshots of Simulation results from Case [1,2,3]-[II]-[(a),(b),(c)].

Fig. 5는 Case [1,2,3]-[III]-[(a),(b),(c)]의 시뮬레이션 결과를 그림으로 나타낸 것이다. III 타입의 경계를 가진 경우에 대한 시뮬레이션 결과는 II 타입의 경계면을 가진 경우와 다르지 않은 결과를 보여주고 있다. 이 경우들의 비교에서 주목할 점은 Atwood 수가 낮은 Case 1의 경우에는 앞선 경우와 같이 경계면에서의 교란이 보였으나, II 타입에서는 기하학 집중부에서 유동이 발달한 것과는 다르게, I 타입과 비슷한 형태인 기하학 특이점에서 경계면의 붕괴가 발생하였다. 그리고 좌측 벽면에 생성된 기하학 집중부에서 무거운 유체의 하강이 아닌 가벼운 유체의 상승 유동이 일어났으며, 이는 Case 2와 3에서도 함께 관찰되었다. 좌측 벽면에서 발생된 상승유동과 바로 그 옆에서 발생한 하강 유동은 높은 상대속도를 생성하였고, 그로인해 Kelvin-Helmholtz 불안정성에 기인한 유동에서 흔히 나타나는 쇄파(Breaking Wave) 형태의 유동이 발생함을 확인하였다. 본 결과에는 보이는 또 다른 현상은 Case [1]-[III]-[(a)]에서 나타난 것과 같이 시뮬레이션 시간 4.0초에서 하강 유동에서 발생된 버섯형태의 유동이 다른 경우들과는 달리 머리 부분이 가늘고 넓게 퍼진 형태를 보이고 있다. 이는 에너지 집중에 의해 RTI에 기인한 버섯형태의 유동을 만들었으나 하부 유체를 관통할 수 있는 충분한 힘이 확보되지 못하여 유동이 하부로 나아가지 못하고 양 측면으로 퍼저나가는 것으로 보인다. 이는 Case [1]-[III]-[(b),(c)]에서도 볼 수 있는데, 하부 유동의 머리 부분에서 회전 운동을 하는 유동이 Case 2와 3에 비해 크게 나타나지 못하고 하강하는 것을 볼 수 있다. 이로써 RTI에 기인한 비선형적인 유동의 발생은 경계면의 형상과 그 진폭에 따른 변화보다는 밀도 차에 의한 변화가 주된 요소임을 확인할 수 있다.

Fig. 5.

Snapshots of Simulation results from Case [1,2,3]-[III]-[(a),(b),(c)].

4. 결론 및 고찰

일반적으로 단일 유체의 유동에 비해 다수의 상이 다른 유체의 유동은 복잡함을 띤다. 특히 유체-유체간의 경계면에서 불안정성에 의한 경계면의 붕괴 및 대변형 등에 대한 해석의 어려움이 존재한다. 본 연구에서는 다상 유동이 가능한 개선된 MPS법(Kim et al., 2014)을 이용하여 밀도가 서로 다른, 즉, Atwood수가 서로 다른, 두 유체간 경계면에서 발생하는 복잡한 문제 중 하나인 RTI에 대해 시뮬레이션 하였다. 밀도차가 서로 다른 경우를 가정하여, 밀도차가 RTI유동의 발달 속도 및 과정에 큰 영향을 미치는 것을 확인하였다. 또한 초기 경계면의 형태에 따른 RTI 유동 발달의 양상을 비교한 결과, RTI 유동 발달의 차이는 다소 있었으나, 밀도차에 의한 결과에 비해 원활하지 않은 유동의 발달을 보였다. 이러한 불안정성의 수치 해석적 정확도의 확인을 통해 향 후 다양한 형태의 다상 유동에 대한 MPS법의 확장성을 확인할 수 있었으며, 이를 통하여 해양공학 분야에 있어 계속적으로 요구되는 복잡한 다상 유동을 포함한 응용문제 해석에 크게 기여할 수 있을 것으로 보인다.

Acknowledgments

본 연구는 산업통산자원부 산업핵심기술개발사업(10051 136, 해양플랜트 구조물의 부유거동 예측을 위한 입자기반 유체-다물체동역학 통합 해석 소프트웨어 개발)의 일환으로 수행하였음.

References

Chen, H. C., (2011), “CFD Simulation of Compressible Two-phase Sloshing Flow in a LNG Tank”, Ocean Sys. Eng., 1(1), p29-55.
Harlow, F. H., (1955), A Machine Calculation Method for Hydrodynamic Problems, Los Alamos Sci. Lab., report LAMS-1956.
Harlow, F. H., and Welch, J. E., (1965), Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Phys. Fluids, 8, p2182-2189. [https://doi.org/10.1063/1.1761178]
Joeng, S. M., Nam, J. W., Hwang, S. C., Park, J. C., and Kim, M. H., (2013), “Numerical prediction of oil amount leaked from a damaged tank using two-dimensional moving particle simulation method”, Ocean Eng., 69, p70-78. [https://doi.org/10.1016/j.oceaneng.2013.05.009]
Khayyer, A., and Gotoh, H., (2013), “Enhancement of Performance and Stability of MPS Mesh-free Particle Method for Multiphase Flows Characterized by High Density Ratios”, J. Comput. Phys., 242(1), p211-233. [https://doi.org/10.1016/j.jcp.2013.02.002]
Koshizuka, S., and Oka, Y., (1996), “Moving-Particle Semi-implicit Method for Fragmentation of Incompressible Fluid”, Numer. Sci. Eng., 123, p421-434.
Kim, K. S., Kim, M. H., and Park, J. C., (2014), “Simulation of Multiliquid-Layer Sloshing with Vessel Motion by using Moving Particle semi-implicit Method”, J. Offshore Mech. Arctic Eng., 137, p051602-1, 2015. [https://doi.org/10.1115/1.4031103]
Lee, B. H., Park, J. C., and Kim, M. H., (2010), “Numerical Simulation of Impact Loads Using a Particle Method”, Ocean Eng., 37, p164-173. [https://doi.org/10.1016/j.oceaneng.2009.12.003]
Moeendarbary, E., Ng, T. Y., and Zangeneh, M., (2009), Dissipative Particle Dynamics: Introduction, Methodlogy and Complex Fluid Applications – A Review, Int. J. Appl. Mech. Eng., 1(4), p737-763. [https://doi.org/10.1142/S1758825109000381]
Monaghan, J. J., (1994), Simulating free surface flows with SPH, Journal of computational physics, 110(2), p399-406. [https://doi.org/10.1006/jcph.1994.1034]
Nomura, K., Koshizuka, S., Oka, A., and Obata, H., (2001), “Numerical Analysis of Droplet Breakup Behavior Using Particle Method”, J. Nucl. Sci. Technol., 38(12), p1057-1064. [https://doi.org/10.1080/18811248.2001.9715136]
Shakibaeinia, A., and Jin, Y.C., (2012), “MPS Mesh-free Particle Method for Multiphase Flows”, Comp. Methods Appl. Mech. Eng., 229, p13-26. [https://doi.org/10.1016/j.cma.2012.03.013]
Sharp, D.H., (1984), “An overview of Rayleigh-Taylor instability”, Physica D: Nonlinear Phenomena, 12(1), p3-18. [https://doi.org/10.1016/0167-2789(84)90510-4]
Shirakawa, N., Horie, H., Yamamoto, Y., and Tsunoyama, S., (2001), “Analysis of the Void Distribution in a Circular Tube with the Two-Fluid Particle Interaction Method”, J. Nucl. Sci. Technol., 38(6), p392-402. [https://doi.org/10.1080/18811248.2001.9715045]
Sussman, M., Smereka, P., and Osher, S., (1994), A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow, J. Comput. Phys., 114, p272-280. [https://doi.org/10.1006/jcph.1994.1155]

Case	ρ₁ [kg/m³]	ρ₂ [kg/m³]	v [s/m²]	σ [N/m]	AT	Case	Initial Perturbation	Case	Amp. (A)
1	600	400	0.01	0.0723	0.2	I	η = Acos(2πx)	(a)	0.025m
2	700	300	0.01	0.0723	0.4	II	η = Acos(4x−0.5)π	(b)	0.050m
3	800	200	0.01	0.0723	0.6	III	η = Acos(πx)	(c)	0.075m