Journal Archive

일반적으로 해양공학에서 저장조(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 불안정성에 대해 다양한 밀도 조합을 이용하여 수치해석을 수행하였으며, 다양한 초기 경계면의 형상에 따른 불안정성으로부터 기인한 유동의 발달 양상에 대해 연구를 수행하였다.

다상 유동과 같이 두 개 이상의 상이 다른 유체가 혼합되어 하나의 시스템에 존재할 때, 복잡한 유체의 운동이 발생하며, 단일의 현상이 아닌 다양한 불안정성에 기인한 경계면 붕괴 및 유체간의 혼합이 일어나기 때문에 이를 해석하기가 쉽지 않다. 이를 위해 반복적으로 발생하는 현상을 관찰하여, 반복성의 확인 후 각각의 불안정성을 명명할 수 있으며, 대표적으로는 시스템의 붕괴 시 경계면의 형상에 따라 나뉠 수 있다. 경계면이 버섯 형태의 유동을 만드는 경우는 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의 경우는 상대적인 에너지의 차가 크지 않아 경계면에서 두 유체가 지속적으로 충돌을 일으키며 경계면을 흐트러트렸기 때문으로 판단된다.

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

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

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

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

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

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

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)].

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