[ Original Article ]

Journal of the Korean Society for Marine Environment and Energy - Vol. 19, No. 2, pp.111-119

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

Print publication date May 2016

Received 15 Jan 2016 Revised 14 Mar 2016 Accepted 14 Mar 2016

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

Lattice-Boltzmann Method를 이용한 2차원 기체-액체간 거동 기초 연구

정노택^†

울산대학교 조선해양공학부

Feasibility Study on the Gas-Liquid Multiphase by Lattice-Boltzmann Method in Two-Dimensions

Rho-Taek Jung^†

School of Naval Architecture & Ocean Engineering, University of Ulsan, Ulsan 44610, Korea

Correspondence to: ^† rtjung@ulsan.ac.kr

초록

기체-액체 이상유동의 거동 시뮬레이션을 위해 Lattice Boltzmann방법(LBM)을 이용하였다. 기체-액체사이의 경계면에서 상호포텐셜 모델인 Shen-Chan방식과 Carnahan-Starling 상태방정식을 도입하였다. 또한 외력항의 처리는 Exact Difference Method를 사용하였다. 개발된 코드를 통하여 상태방정식 특성파악, 기체-액체의 상분리, 표면장력 및 기체액체 경계면 거동 특성, Homogeneous와 Heterogeneous 캐비테이션, 기포 붕괴등의 시뮬레이션을 수행하였다.

Abstract

Gas-Liquid multiphase flow simulation has been carried out using the Lattice boltzmann method. For the interface treatment, pseudo-potential model (Shan-Chen) was used with the Carnahan-Starling equation of state. Exact Difference Method also applied for the treatment of the force term. Through the developed code, we simulated coexsitence structure of high and low density, phase separation, surface tension effect, characteristics of moving interface, homogeneous and heterogeneous cavitation and bubble collaps.

Keywords:

Lattice-Boltzmann Method, Phase separation, Bubble Cavitation, D2Q9

키워드:

격자볼츠만방법, 상분리, 버블 캐비테이션, 2차원9방향

1. 서 론

격자 볼츠만법(Lattice Boltzmann Method: LBM)은 마이크로스케일에서 메조스케일을 거쳐 마크로 스케일까지, 단일상 뿐만 아니라 다상류의 유체유동등 응용사례가 확장되고 있다(Succi, S., [2008, 2015], Chen[2003]). Kinetic 방정식(Boltzmann)으로부터 유체입자의 움직임을 입자밀도함수로 표현한 LBM은 Cellular Automata와의 정합성 뿐만 아니라 연속체로 가정하여 유도되는 Hydrodynamic 방정식(Navier-Stokes)으로 전개가 가능하기 때문에 공학적으로 유용한 전산해석기법이 되고 있다(He and Luo[1997]).

LBM은 입자밀도함수의 충돌항과 전진항으로 이루어져 있어 비교적 간략한 식으로 구성되어 이산화가 용이하다는 점, 대류항이 없다는 점, 기체-액체간 경계 처리가 용이하다는점, 그리고 병렬화 수행이 뛰어나다는 점등의 특징이 있다.

기체-액체 이상간의 시뮬레이션은 LBM 영역에서도 많은 연구자들의 관심사안으로, 기체-액체 유체입자가 경계면에서 상호 작용하는 힘을 물리적으로 모델링하여 진행하게 된다. 다상류 시뮬레이션에서 가장 중요한 요소중의 하나가 그 경계면의 처리일 것이다. 이 또한 Navier-Stokes기반 시뮬레이션처럼 LBM에서도 다양한 방법들이 제안되고 있다. 다상류 해석시 세가지 정도의 모델이 독자적으로 연구개발되고 있다(Chen and Doolen [1998]).

첫째, Color gradient mode로 Gunstensen et al.[1991]에 의해서 개발되었으며, i방향성을 가지는 (red)(blue) 두 개의 국소밀도분포 함수 f_i^(red)와 f_i^(blue)로 나누어서 기체-액체의 시뮬레이션을 수행하였다. 경계면 가까이에서는 두함수의 충돌항의 일부가 혼합 존재하며, 경계에서 국소 구배를 만족하도록 모델을 개발하였다. Lishchuk et al.[2003]는 외력항에 밀도의 구배를 추가하여 모델 정도를 향상시켰다. 이 Color gradient model 모델은 질량과 모멘트가 보존된다. 둘째, Pseudo-potential model로서 Shan-Chen모델(Shan, X et al.[1993, 1994])로 알려져 있으며, 기체-액체 입자간의 충돌을 기체-액체간의 상호충돌력으로 표현하고, 이 상호 충돌력은 그린함수와 유효질량함수의 관계로 설정한다. 여기서 상태방정식으로 표현되는 유효질량함수는 밀도의 함수이며, 임계값 G에 의해서 상분리의 유무를 판단하게 된다. 이 방식은 국소적으로는 모멘텀이 보존되지 않으나 영역전체에 대해서는 보존된다. Falcucci et al.[2006]은 주변 셀 정보를 추가한 two-belt로 모델의 안정성을 더 높혔다. 셋째, Free energy model(Swift et al.[1995,1996])은 Cahn-Hilliard[1958]의 자유에너지를 압력텐서에 직접적으로 대입한 것으로 자유에너지는 밀도, 밀도구배의 제곱승과 표면장력으로 구성되게 된다. 질량과 모멘텀이 보존되며, 모델 상에 나타난 Galilean invariance의 문제는 Inamura [2000]에 의해 논의되었다. Lee et al.[2005]는 Free energy 모델을 바탕으로 이산화를 엄밀히 하여 기체-액체 경계면에서 발생하는 이상 속도(parasitic velocity)가 없는 방식을 개발했다. 1990년대 초반부터 2000년대 중반에 걸쳐 개발된 LBM 다상모델은 개별 특징을 가지며 지속적으로 발전하고 있다(Yu et al. [2003]).

본 논문은 위 세가지 LBM 다상류 시뮬레이션 중에서 Pseudopotential model을 이용하여, 기체-액체 캐비테이션 연구의 활용가능성을 확인하기 위하여 코드를 개발하였으며, 기본 검증을 연구결과로 나타내었다. 기체-액체 Shan-Chen LBM(SC LBM)방법을 이용하여 다음과 같은 시뮬레이션을 수행하였다. 먼저 본 논문에 사용된 지배방정식의 공간정확성을 파악하고, Carnahan-Starling(CS) EOS의 특성, 기체-액체의 상분리(3.1절), 표면장력 특성(3.2절), 기체액체 경계면 특성(3.3절), Homogeneous(3.4절)와 Heterogeneous 캐비테이션(3.5절), 기포 붕괴(3.6절)등의 순으로 논한다.

2. Lattice-Boltzmann Method

이상유체 유동을 해석하기 위해 Shan-Chen모델을 사용하였다. 국소분포함수 f_i와 격자의 이동속도 e_i를 사용하여 외력항 F_i가 존재시에 볼츠만 방정식을 격자 볼츠만 법으로 표현하면 식 (1)과 같다(Huang et al.[2011]).

f i x + e i, t + ∆ t = f i x, t + Ω i f i x, t + F i

(1)

오른쪽 두 번째항은 분포함수 f_i의 i방향간에 국소 충돌함수이며, 세 번째항은 외력항이다. 본 시뮬레이션에서 Δt는 시간스텝으로 1로 설정하였으며, D2Q9(2차원 9방향 단위속도성분)방식을 채택했다. 국소충돌함수로서 Bhatnagar et al.[1954]이 제안한 Single-timerelaxation process를 적용하면 Ω_i= -τ^-1 (f_i - f_i^eq)로 대처할 수 있으며, f^eq는 국소등가분포함수로서 유체변수값 즉, 밀도와 속도의 함수이다. 유체동점성계수를 결정하는 τ = v/c_s² + Δt/2, sound of speed c_s= c/ $3$ , c = Δx/Δt = 1.0이며, Δx는 계산시의 단위격자로서 1로 설정하였다. τ 값을 적절히 선택함에 따라 2차정도의 비압축성유체를 해석할 수 있다. Fig. 1은 2차원 사각형 격자에 9방향(D2Q9)을 택하였으며, LB법은 분포함수의 충돌스텝 -Ω = F_i과 i방향으로의 전진스텝 f_i(x + e_i,t + Δt) = f_i(x,t)으로 나누어지며, 국소등가분포함수 f_i^eq는 다음과 같이 주어진다.

f i e q = ρ w i 1 + 3 e i ∙ u / c 2 + 9 e i ∙ u 2 / 2 c 4 - 3 u ∙ u / 2 c 2

(2)

Fig. 1.

Momentum distribution in two dimensional nine directions (D2Q9).

w_i값은 i가 0일 경우는 4/9, 1~4일경우는 1/9, 5~8일 경우는 1/36의 값을 가진다. 마크로 스케일에서의 물리량인 밀도와 속도는 아래와 같이 주어진다.

ρ = ∑ i f i, ρ u = ∑ i f i e i

(3)

물-공기 이상류의 인터페이스 상에 존재하는 국소분포함수의 힘을 외력으로 정의하면 다음과 같은 식을 사용한다.

F i x, t = - G ψ x, t ∑ i w i ψ x + e i ∆ t, t e i

(4)

여기서 G는 인터페이스상의 세기를 나타내는 인자이며, w_i의 값은 국소등가분포함수식 내의 값과 동일하다. ψ는 상호 포텐셜로 밀도의 함수이며, 밀도가 증가할수록 단조증가하는 것으로 가정한다. 또한 이는 온도가 변하지 않는 항온상태내에서만 가능하다(Shan and Chen[1994]). 상호 포텐셜 ψ는 다음과 같이 정의한다.

ψ ρ = ψ 0 e - ρ 0 / ρ

(5)

식 (4)와 (5)에 의해서 시뮬레이션을 위한 비이상기체(nonideal gas) 방정식이 유도된다(He and Doolen [2002]).

p = ρ R T + G R T 2 ψ 2 ρ

(6)

여기서, RT = 1/3이며, 식 (6)의 오른쪽 첫 번째 항은 이상기체, 두번째항은 비이상기체로 입자사이에 주어지는 힘과 관련된 항으로 G가 0보다 작을 경우에 압력이 감소되는 효과를 가져온다. 비선형 상태방정식인 식 (6)은 G의 임계값 이하가 되면 상분리가 발생하는 것으로 가정한다. G는 상호작용시에 강세를 나타내며 임계치를 기준으로 기체와 액체의 영역이 달라진다. 밀도의 함수로 주어지는 ψ(ρ)는 effective mass의 역할을 하는 변수로, 밀도가 증가하면 ψ(ρ)도 증가된다. 다시 상호포텐셜항은

ψ ρ = 2 p - ρ c s 2 G c s 2 c s = R T, R T = 1 / 3

(7)

로 나타낼 수 있으며, 본 논문에서는 압력항 p를 Carnahan-Starling 상태방정식을 사용하였다.

p = ρ R T 1 + k + k 2 + k 3 1 - k 3 - a ρ 2, k = b ρ r,

(8)

여기서 a = 0.4963(RT_c )² /p_c, b = 0.1873RT_c/p_c 이며, T_c와 p_c는 임계온도, 임계압력을 나타낸다. 본 연구에서는 (a,b,R) = (1,4,1)을 적용하였다. 본 격자 시뮬레이션에서 적용된 임계값은 (T_c, ρ_c, p_c) = (0.0943, 0.1136, 0.0044)이다.

식 (4)와 달리 기존 상호포텐셜에 2e_iΔt 위치에서의 상호포텐셜을 추가적으로 포함하여 다음과 같은 식을 제안 하였다(Sbragaglia et al.[2007]).

F i = - ψ x, t ∑ i w i g 1 ψ x + e i ∆ t, t + g 2 ψ x + 2 e i ∆ t, t e i

(9)

이때의 $ψ ρ = 2 p - ρ c s 2 / c s 2 g 1 + 2 g 2$ , 기체-액체 경계면에서 작용하는 상호포텐셜의 g₁, g₂ 값에 대해서는 3.3절에 나타내었다.

LBM에 중력장, 가속도항, 전자기장등과 같은 외력항 처리에 대해서 많은 연구가 진행되었다(Huang, et al.[2011]). 운동방정식 F = mdu/dt에 의해서 Δu = FΔt/ρ으로 변경가능하며, 이것은 전진 스텝시에 외력으로 주어진 F는 격자내에서 Δu 만큼 이동한다는 것을 의미한다. 또한 LBM상에서의 외력항의 구성은 속도장내에서 f_i^eq의 변화량 Δf_i = f_i^eq(ρ,u +Δu) - f_i^eq (ρ, u) 으로 표현한 Exact Difference Method (EDM, Kupershtokh et al. [2009]) 방식을 채택했다. 식 (4)를 이용하여 외력벡터를 구한 후 식 (11)에 따라 외력항을 구성하게 된다.

F i = f i e q ρ, u + F ρ ∆ t - f i e q ρ, u

(11)

본 LBM논문에서 사용되는 단위는 Lattice 단위로 사용된다. 숫자뒤에 단위가 표기 되어 있지 않을 경우에는 Lattice 단위로 사용된 것이며, 이것은 길이[lu], 시간[ts], 질량[mu], 밀도[mu/lu³], 압력[mu/(lu×ts²)], 표면장력 [mu/ts²] 온도[tu] 로 나타낼 수 있다. 실제 물리량과 비교시에는 무차원값으로 비교가능하다. 밀도의 경우ρ = ρ^la/ρ_c^la = ρ^ph/ρ_c^ph, 압력의 경우는 p = p^la/p_c^la = p^ph/p_c^ph, 온도의 경우는 T = T^la/T_c^la = T^ph/T_c^ph으로 윗첨자 la는 LBM 격자 단위, 윗첨자 ph는 실제 물리량, 아래 첨자 c는 임계값을 의미한다.

3. 계산결과

일반적으로 잘 알려진 LBM의 지배방정식(1)은 1차정도의 이산화형태로 나타나 있음에도 불구하고 공간적인 정확도는 2차정도로 알려져 있다(Huang et al. [2011]). 먼저 지배방정식의 신뢰성을 확보하기 위하여 공간상의 정확성을 위해 기체-액체의 다상류 현상을 공간상에 구현할 때 2차정도의 정확도를 가지는지 확인하기 위해 Fig. 2에 격자수에 따른 밀도공간의 정확도를 확인하였다.

격자는 Nx×Ny로 50×50, 100×100, 200×200, 300×300. 400×400에 대해 높은 밀도(ρ_h)와 낮은 밀도(ρ_l)의 정확도를 400×400격자에서의 밀도값으로 비교해 보았다. 해당 에러율은 Error(Nx)=|밀도(Nx)-밀도(400)|로 정의하였다. Fig. 2에서 보이듯이 ρ_h와 ρ_l 각각의 에러율이 포물선형태인 것을 알 수 있고, 격자수가 증가할수록 ρ_h의 경우가 ρ_l경우보다 급격히 에러율이 낮아진다. 본 논문에서는 100×100의 격자수로도 유의한 결과를 얻을 수 있다고 판단하여 이후 동일한 격자수에서 계산을 수행하였다.

본 논문에서는 식 (7)의 자유에너지 CS EOS를 사용하고 있다. 상태방정식(Equation of state)은 상분리를 재현하기 위해 사용되어진다. 3차곡선으로 표현되는 상태방정식은 압력에 대한 밀도의 구배가 음수일 때는 물리적이지 않으며, 양수일때만 물리적 의미를 가진다. 그 구배 값이 0이 되는 두 변곡점 외부의 밀도값이 그 상태에 따라 낮은 밀도(ρ_l)과 높은 밀도(ρ_h)로 재배치 되면서 상분리가 이루진다. 이때의 해당 압력값(ρ_l, ρ_h)도 구해지게 된다. Fig. 3은 온도가 주어진 상태에서 2차원 밀도 랜덤 시뮬레이션을 통하여 10,000 ts에서 얻은 ρ_l와 ρ_h 값을 밀도 및 온도 임계값에 대해 무차원화(ρ = ρ/ρ_c, T = T/T_c) 한 값이다(Yuan et al.[2006]). 실선은 ES EOS의 Maxwell 등면적 구성으로 구한 값이다. T의 값에 따라 등가의 밀도값이 ρ_l와 ρ_h이 각각 정해지며, Maxwell 실선과 비교하면 기체 쪽은 잘 따르고 있으나, 액체 쪽이 약간 밀도가 크게 나타났다. 계산상 τ는 1, g₁, g₂는 -0.33과 0.0을 사용하였다.

3.1 2차원 기체-액체 상분리

저밀도(ρ_l)와 고밀도(ρ_h)가 2차원적으로 상분리 되는 시뮬레이션을 수행하였다. 주어진 온도조건하에서 랜덤 (밀도)값을 초기값으로 분포시키고, 계산 진행시 그 밀도값은 해당 EOS 방정식을 통하여 비물리적이고 불안정한 상태에 있는 밀도값이 각각 안정한 값을 찾으면서 고밀도와 저밀도로 상분리가 일어나게 된다. 상분리가 진행되는 동안 응축과 증발은 없다고 가정한다. 100×100격자, (g₁, g₂) = (-0.33, 0)으로 설정하였으며, 이 상분리는 인터페이스의 면적이 최소화하는 방향으로 진행이 된다(Fig. 4). 계산이 진행되면서 입자 간의 상호작용에 의해서 기포가 합쳐지는 현상이 일어나며, 최종적으로 한 개의 원으로 안정화 된다. 이것은 자유에너지(free energy)의 최소화 과정에 있으며, 시뮬레이션 대상이 기포인지 액적인지는 초기밀도값 분포시에 결정된다(본 논문에서의 기포는 기체-액체 경계면 안쪽이 ρ_l, 바깥쪽이 ρ_h를 나타내며, 액적은 반대인 경우이다).

특히 Fig. 4에서 경계면의 두께가 계산이 진행되면서 일정함을 알 수 있고 세 개의 등고선은 (ρ_l, ρ_h) = (0.02, 0.33) 경계면 (ρ_l+ ρ_v)/2의 값을 나타내었다. 주어진 온도 T 0.75에서 두상이 고밀도와 저밀도로 안정화되어 평형상태에 도달한다.

Fig. 2.

spacial accuracy of ρh and ρl.

Fig. 3.

Coexistance density of gas and liquid of CS EOS along temperatures with Maxwell equal-area construction shown by solid line.

3.2 표면장력

기체-액체 상분리 시뮬레이션을 이용하여 상분리로 인해 액적 또는 기포가 유지 될 때에 표면장력을 추정하였다. CS LBM에서 표면장력과 액적 또는 기포의 반지름과 압력차의 관계(라플라스 법칙)를 파악하였다. 단일 액적과 단일 기포 구성을 위해 액적의 경우는 초기분포를 (0.0~0.15), 기포의 경우에는 (1.0~4.0)의 범위로 랜덤 수를 형성하였다. T 0.75으로 고정하고 격자수를 증가시키며 최종 R을 100,000 ts에서 측정하였다. 이때 R는 경계면 중앙(ρ_l +ρ_h)/2)에서 원중심까지의 거리이며, 50×50, 70×70, 100×100, 200×200을 테스트 하였다. 다양한 크기의 액적 또는 기포로 표면장력을 추정할 수 있으며, 시뮬레이션으로 얻어지는 경계면 내외의 압력과 그때의 액적 또는 기포의 크기를 시뮬레이션을 통해 얻을 수 있다. Fig. 5(a)에서 보이는 바와 같이 계산은 표면장력의 라플라스식을 잘 따르고 있음을 확인할 수 있다. 액적 또는 기포의 반경에 따라 EOS 압력차를 그래프화 하였다. 실선은 y = 0.01x + 0.00001을 나타내며, T 0.75에서 표면장력값은 0.01이다. Fig. 5(a) 흰색원은 기포를 검은색원은 액적을 나타낸다. 동일 반경에 대해서 기포가 액적보다 약간 큰 압력차값을 나타낸다. Fig. 5(b) 온도에 따른 표면장력의 변화를 확인하기 위하여 기포 반경 R₀를 40으로 고정시킨 상태로 시뮬레이션 하였다. 온도 값이 내려 갈수록 표면장력이 더욱 커진다. 이때의 단일기포의 밀도분포를 식 (12)와 같이 나타내었다(Huang et al. [2015]).

ρ x, y = ρ l + ρ h 2 + ρ l - ρ h 2 1 - tanh ⁡ 2 x - x c 2 + y - y c 2 - R o W

(12)

ρ_l와 ρ_h는 기포와 액체의 밀도, x_c, y_c는 기포의 중심, W는 경계면의 두께를 나타낸다. R₀는 중심에서 (ρ_l+ρ_h)/2 까지의 거리를 나타낸다.

3.3 기체-액체 경계면 특성

SC 모델에서 기체-액체 시뮬레이션에서 경계면의 거동을 2차원 100×100격자내에 단일 기포의 밀도분포 식 (12)와 같이 구성하고 초기 경계면 두께 W를 4로 두었으며, W의 값이 작을수록 계산이 불안정하게 된다. 사각형의 계산영역에서 사방 모두 주기경계조건으로 설정되었다. 정적상태에서 기체-액체 경계에서 밀도값의 변화를 식 (9)내의 g값의 역할에 대해 조사하였다. g₁= −0.33으로 동일하며, g₂= 0, -0.33, -50 세가지의 경우에 대해 인터페이스 거동 정도를 파악하였다. 모든 경우에 대하여 초기형태를 제외하고는 10,000 ts에서의 값을 플롯하였다. Fig. 6(a)는 g₂ 값이 작아짐에 기체-액체 경계면을 잘 따르는 것 같아 보이나 초기 ρ_h 값과는 차이가 나며 (Fig. 6(b)) 아래면에서는 경계면 뿐만 아니라 초기 ρ_l 값과도 차이가 나타났다(Fig. 6(c)). 반면, g₂가 0일 경우는 오히려 경계면의 구배가 더욱 커지며, 초기 ρ_h과 ρ_l을 잘 추적하고 있다. 따라서 본 논문에서는 (g₁, g₂) = (-0.033, 0)을 사용한다.

Fig. 4.

Phase separation to ρl (inside) and ρh (outside) from a random value of the density (average density value is 0.175).

Fig. 5.

Surface tension ((a) laplace relation, (b) tension with non dimensional temp).

Fig. 6.

Interface consistency and construction in one dimension ((a) density profile along Nx 100 (b) and (c) density zoom part).

Fig. 7.

(a) Homogeneous cavitation settling down to the equilibrium density value (b) density profile along the Nx.

3.4 Homogeneous 캐비테이션

시뮬레이션 대상 유체가 순수 유체라고 가정하고 단일 밀도를 분포시킨 후 일정속도로 양측에서 잡아당길 경우 유체내에 캐비테이션이 어떻게 발생하는 지를 확인하였다. 유체 장력의 한계로부터 발생하는 캐비테이션(homogeneous)은 핵 캐비테이션(nucleus cavitation) 연구에 중요하다. 유체장내에서 형성되는 기포의 임계 반경 R_c과 그 압력차 Δp가 핵 캐비테이션을 특성짓기 때문이다(Sukop et al.[2005]).

격자는 100×100, T/T_c는 0.75에서 CS EOS의 spinodal 밀도값보다 조금 높은 값인 0.3을 초기밀도로 설정하고 시뮬레이션을 수행하였다. 위 아래면에서 속도값(u_w= 0.005)과 Zou and He[1997] 경계조건을 적용하고, 좌우는 반복조건을 주었다. 즉 식 (3)에 의해 계산영역내로 들어오는 방향밀도 함수인 f₄, f₇, f₈을 위 경계면에서 구하고 f₂, f₅, f₆는 아래 경계면에서 계산하게 된다. 계산 중 밀도변화의 모니터링을 위해서 높은 밀도가 예상되는 속도 경계면 근방[50,1]과 낮은 밀도가 예상되는 지점[50,50]을 택하여 모니터링 하였다. Fig. 7(a)는 밀도 0.3인 유체는 계산이 진행됨에 따라 spinonal 밀도값(0.25)으로 접근하고 이 때 spinodal 최대 표면장력을 이기면서 고밀도와 저밀도로 급격하게 분리된다. 바로 초기진동이 발생하며, 이후 진동이 감쇄되면서 등가의 밀도로 수렴한다. 고밀도 ρ_h의 경우가 비압축성의 영향으로 저밀도 ρ_l의 경우보다 진동이 더 커지는 결과를 나타내고 있다.

3.5 Heterogeneous 캐비테이션

에너지 크기에 따라 어떠한 형상으로 변화하는지를 시뮬레이션하였다. 마이크로 기포는 기포의 임계반경과 압력에 의해서 그 특징이 결정된다. 기포 에너지는 2차원의 경우 표면장력에 의한 에너지와 기포 내외 압력차로 인한 에너지로 나타낼 수 있다(Or et al.[2002]).

∆ E = 2 π r σ + π r 2 ∆ p

(13)

Fig. 8(a)는 식 (13)을 플롯한 것이며, 버블의 임계반경이 각각 20, 30, 40이며, 이때의 임계에너지 ΔE_c는 0.454, 1.353, 2.892이다. 유체가 음의 압력을 가질 때 경계면의 에너지는 임계 반경까지 버블반경과 동시에 증가한다. 이는 버블 표면 에너지가 증가하기 때문이다. 반면 임계반경을 넘어서면 반경은 증가하나 에너지는 감소한다. 음의 압력이 상대적으로 더 크게 작용하기 때문이다. 예를 들어 Δp가 0.00388119로 설정하고 T 0.75일때 R₀는 25.76인 된다. R₀보다 작은 경우 R₀=20와 큰 경우 R₀=30을 대상으로 ΔE가 증감하는 시뮬레이션 결과를 Fig. 8(b)(c)에 각각 나타내었다. 초기 반경이 20일 경우 Δp는 0.005이며, 이에 해당하는 기포의 밀도와 액체장의 밀도를 (ρ_l, ρ_h) = (0.00806502, 0.331026)으로 설정하였고 초기반경이 30일 경우 Δp=0.0003333이며, ρ_l만 변경하여 (ρ_l, ρ_h) = (0.01097479, 0.331026)로 설정하였다. 결과치를 보면 음의 압력으로 초기반경 20의 경우는 임계반경이내임에 따라 점점 줄어들고, 30의 경우는 증가한다. 특히 두 경우 증감이 가속화 되는 것은 동일한 현상이나 버블이 기체 쪽으로 감속하는 속도가 버블이 액체쪽으로 증가하는 속도 보다 훨씬 빠르게 나타났다.

3.6 기포 붕괴

2차원 기포붕괴 시뮬레이션을 다음과 같이 수행하였다. 격자는 400×400, T/T_c=0.75, Δp=0.0572이며, 이는 Δp/p_c=1.3에 해당한다. 1400ts까지 200ts간격으로 플롯하였다. 아래면의 경계조건은 고체벽 조건을 주고, 나머지 면의 경계조건은 동일하게 Zou and He[1997]경계조건을 주었다. 기존 논문과 기포 붕괴 형상비교를 위해서 b/R_max를 1.5로 두었다(Plesset et al. [1997]). 여기서 b는 초기 반경 R_max의 원 중심부터 아래 벽면까지의 거리이다. 결과값은 1400ts 까지 200ts 간격으로 Fig. 9(b)에 나타내었다(a, b, c, d, e, f, g, h, i)=(initial, 200, 400, 600, 800, 1000, 1200, 1400). 좌우와 위경계면으로부터 주어진 일정 압력으로 인하여 기포 상부에서 점점 아래로 줄어들어가는 형상을 보인다. g라인 부터는 기포 경계면에 외곡이 진행된다. 기존 논문의 형상과 비교해보면 LBM 시뮬레이션에서도 고체경계면이 존재시에 고체면과 반대편의 기포가 붕괴되는 현상을 제현할 수 있는 것으로 판단된다. 기포 표면은 외각의 높은 압력과 유체에서 기포 경계면을 통하여 전파된 속도가 기포 내부의 더 큰 속도로 바뀌어 표면의 외곡을 유발하는 것을 Fig. 9(b)를 통하여 확인할 수 있다. 기포가 붕괴시점이 되면 고체면으로 향하여 가속되어 지는 힘이 발생한다. Fig. 10(a)(b)에 나타낸 바와 같이 g와 i 시간에서 속도장을 표현한 것으로 속도가 커진 것을 알 수 있다. 이는 3.7절의 Heterogeneous 캐비테이션의 기체 움직임이 액체움직임보다 빠른 것과 동일하다. 그리고 기포 경계면에서 속도가 존재함을 확인 할 수 있으며, 이는 기포 붕괴와 같은 hydrodynamic 속도가 빠른 경우에는 무시한다고 가정해도 유용한 결과 값을 얻을 수 있을 것으로 판단된다(Yang et al.[2015]).

Fig. 8.

(a) ΔE vs. radius (b) bubble shrink onset R0 = 20 and (c) bubble enlarge on set R0 = 30 with time evolution.

Fig. 9.

A collapsing bubble in theory Plesset et al., [1997] presented in solid line (a) and dots in experience by Lauterborn et al. [1975] (a) and LBM present simulation in (b) (a, b, c, d, e, f, g, h, i)=(initial, 200, 400, 600, 800, 1000, 1200, 1400).

Fig. 10.

Vector fields around a collapsing bubble.

4. 결 론

격자볼츠만법의 비이상기체모델 중 Psuedo potential 방법(Shan-Chen)을 이용하여 개발된 기체-액체간 시뮬레이션 코드로 기본 예제를 통하여 코드의 신뢰성을 확보하였다. 특히 경계면에서는 Carnahan-Starling상태방정식의 특성에 따라 기체-액체 분리가 잘 이루어지는 것을 확인하였고, homogeneous 캐비테이션과 heterogeneous 캐비테이션 시뮬레이션에도 밀도의 진동이나 경계면의 확산 없이 적절히 수행되었다. 특히 버블 붕괴 시뮬레이션에서는 기존 실험값의 비교해 본 결과 붕괴 형태가 정성적 잘 따라 간다는 것을 파악했다.

본 Psuedo potential 모델이 다양한 확장성과 모델의 이산화 이식성이 좋음에도 불구하고, 국소적으로 모멘텀이 보존하지 않는 것과 이상속도(spurious velocity, 10^-2 정도)값이 비교적 크게 나타나는 것이다. 버블 붕괴와 같은 비교적 큰 유속을 다루는 시뮬레이션에서는 영향이 적을 것으로 판단하나, 본 이상속도가 수치해에 영향을 끼치지 않는다고 볼수는 없기 때문에 시뮬레이션시 주의가 요구된다. 향후 파라메터 연구를 통해서 시뮬레이션 결과가 보다 정량적인 비교가 될 수 있도록 계속 검증해 나갈 필요가 있다.

Acknowledgments

LBM 시뮬레이션의 확장가능성을 보여주신 이태훈(The City College of New York)교수님께 감사드립니다. 또한 울산대학교의 지원에 감사드립니다.

References

Bhatnagar, P.L., Gross, E.P., and Krook, M., (1954), “A model for collision processes in gases. I: small amplitude processes in charged and neutral one-component system”, Phys. Rev, 94, p511-525. [https://doi.org/10.1103/physrev.94.511]
Cahn, J.W., and Hilliard, J.E., (1958), “Free energy of nonuniform system. I. Interfacial Free Energy”, J. Chem. Phys, 28, p258. [https://doi.org/10.1063/1.1744102]
Chen, H., Kandasamy, S., Orszag, S., Shock, R., Succi, S., and Yakhot, V., (2003), “Extended Boltzmann Kinetic Equation for Turbulent Flows”, Science, 301, p633-636. [https://doi.org/10.1126/science.1085048]
Chen, S., Chen, H., Martinez, D., and Matthaeus, W., (1991), “Lattice Boltzmann model for simualtion of magnetohydrodynamics”, Phys. Rev. Lett, 67, p2776. [https://doi.org/10.1103/PhysRevLett.67.3776]
Chen, S., and Doolen, G.D., (1998), “Lattice Boltzmann Method for fluid flows”, Annu. Rev. Fluid Mech, 30, p329-364. [https://doi.org/10.1146/annurev.fluid.30.1.329]
Falcucci, G., Bella, G., Chiatti, G., Chibbaro, S., Sbragaglia, M., and Succi, S., (2006), “Lattice Boltzmann models with mid range interactions”, Commun. Comput. Phys, 1(1), p1-13.
Gunstensen, A.K., and Rothman, D.H., (1991), “Lattice Boltzmann model of immiscible fluids”, Phys. Rev. A, 43(8), p4320-4327.
He, X., and Luo, L.-S., (1997), “Lattice Bolzmann Model for the Incompressible Navier-Stokes Equation”, J. Stat. Phys, 88, p927-944. [https://doi.org/10.1023/B:JOSS.0000015179.12689.e4]
Huang, H., Krafczyk, M., and Lu, X., (2011), “Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models”, Phys. Rev. E, 84, p046710. [https://doi.org/10.1103/PhysRevE.84.046710]
Huang, H., Sukop, M., Lu, X., (2015), “Multiphase Lattice Boltzmann Methods: Theory and Application”, WILEY Blackwell. [https://doi.org/10.1002/9781118971451]
Inamuro, T., Konishi, N., and Ogino, F., (2000), “A Galiean invariant model of the lattice Boltzmann method for multiphase fluids using free-energy approach”, Comput. Phys. Commun, 29, p32-45. [https://doi.org/10.1016/S0010-4655(00)00090-4]
Kupershtokh, A.L., Medvedev, D.A., and Karpov, D.I., (2009), “On equations of state in a lattice Boltzmann method”, Comput. Math. Appl, 58, p965-974. [https://doi.org/10.1016/j.camwa.2009.02.024]
Lauterborn, W., and Bolle, H., (1975), “Experimental investigations of cavitation bubble collapse in the neighborhood of a solid boundary”, J. Fluid Mech, 72, p391-399. [https://doi.org/10.1017/S0022112075003448]
Lee, T., and Lin, C.-L., (2005), “A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio”, J. Comput. Phys, 206, p16-47. [https://doi.org/10.1016/j.jcp.2004.12.001]
Lishchuk, S.V., Care, C.M., and Halliday, I., (2003), “Lattice Boltzmann algorithm for surface tension with greatly reduced microcurrents”, Phys. Rev. E, 67, p036701. [https://doi.org/10.1103/PhysRevE.67.036701]
Nourgaliev, R.R., Dinh, T.N., and Sehgal, B.R., (2002), “On lattice Boltzmann modeling of phase transition in an isothermal non-ideal fluid”, Nucl. Eng. Des, 211, p153-171. [https://doi.org/10.1016/S0029-5493(01)00435-6]
Or, D., and Tuller, M., (2002), “Cavitation suring desaturation of porous media under tension”, Water Resour. Res, 38(5), p1061. [https://doi.org/10.1029/2001WR000282]
Plesset, M.S., and Prosperretti, A., (1997), “Bubble dynamics and cavitation”, Annu. Rev. Fluid Mech, 9, p145-185. [https://doi.org/10.1146/annurev.fl.09.010177.001045]
Sbragaglia, M., Benzi, R., Biferale, L., Succi, S., Sugiyama, K., and Toschi, F., (2007), “Generalized lattice Bolzmann method with multirage psudopotential”, Phys. Rev. E, 75, p026702. [https://doi.org/10.1103/PhysRevE.75.026702]
Shan, X., and Chen, H., (1993), “Lattice boltzmann model for simulating flows with multiple phased and components”, Phys. Rev. E, 47, p1825-1819. [https://doi.org/10.1103/PhysRevE.47.1815]
Shan, X., and Chen, H., (1994), “Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation”, Phys. Rev. E, 49, p42941-2948. [https://doi.org/10.1103/PhysRevE.49.2941]
Swift, M.R., Osborn, W.R., and Yeomans, J.M., (1995), “Lattice Boltzmann Simulation of Nonideal Fluids”, Phys. Rev. Lett, 75, p830-833. [https://doi.org/10.1103/PhysRevLett.75.830]
Swift, M.R., Orlandini, E., Osborn, W.R., and Yeomans, J.M., (1996), “Lattice Boltzmann simulations of liquid-gas and binary fluid systems”, Phys. Rev. E, 54, p5041-5052. [https://doi.org/10.1103/PhysRevE.54.5041]
Succi, S., (2008), 2008, “Lattice Boltzmann across scales: from turbulence to DNA translocation”, Eur. Phys. J. B, 64, p471-479. [https://doi.org/10.1140/epjb/e2008-00067-3]
Succi, S., (2015), “Lattice Bolzmann 2038”, EPL, 109(5), p50001. [https://doi.org/10.1209/0295-5075/109/50001]
Sukop, M.C., and Or, D., (2005), “Lattice boltzmann method for homogeneous and heterogeneous cavitation”, Phys. Rev. E, 71, p046703. [https://doi.org/10.1103/PhysRevE.71.046703]
Yang, J., Shen, Z., Zheng, X., and Li, L., (2015), “Simulation on Cavitation Bubble collapsing with Lattice Boltzmann Method”, J. Appl. Math. Phys, 3, p947-955. [https://doi.org/10.4236/jamp.2015.38116]
Yu, D., Mei, R., Luo, L.-S., and Shyy, W., (2003), “Viscous flow computaions with the method of lattice Boltzmann equation”, Progr. Aero. Sci, 39, p329-367. [https://doi.org/10.1016/S0376-0421(03)00003-4]
Yuan, P., and Laura, S., (2006), “Equations of State in a Lattice Boltzmann Model”, Phys. Fluids, 18, p042101. [https://doi.org/10.1063/1.2187070]
Zou, Q., and He, X., (1997), “On pressure and velocity boundary conditions for the lattice boltzmann BGK model”, Phys. Fluids, 9, p1597. [https://doi.org/10.1063/1.869307]