About Journal of the Korean Society for Marine Environment & Energy

공간과 시간에 따라 변화하는 유체현상을 시뮬레이션 할 경우, 관심영역(공간스케일, 시간스케일)을 정의하게 된다. 주로 유체 시뮬레이션시 활용되는 Navier-Stokes방정식은 비교적 작은 Knudsen수¹⁾의 범위에 해당하며, 이것은 실제 공학적으로 관심을 가지는 스케일이다. 컴퓨터 파워가 증가됨에 따라 보다 작은 계산영역에서 유체해석을 시도하고 있지만 입자 운동학적 효과가 주요한 현상을 시뮬레이션하기 위해서는 한계영역이 존재 한다.

Lattice Boltzmann Method(격자볼츠만법)는 Naiver-Stokes방정식의 공학적 스케일을 포함할 뿐만 아니라 입자 스케일까지도 시뮬레이션이 가능하다(분자단위 제외). Lattice Boltzmann Equation(볼츠만방정식)으로부터 이산화 된 격자볼츠만법은 보다 넓은 Knudsen수의 범위를 포함하고 있고 단순한 이산화과정으로 인해 전산코드개발이 쉬운 장점이 있어 원자력 분야의 비등현상(Kato et al.[1997]), 폼 기하학([Korner et al.[2005]), 결정체동력학(Lu et al.[2009])등의 Meso 스케일 분야 뿐만 아니고, 윈드터빈 블레이드 공력해석(Kim et al.[2012])등에도 유용한 전산해석법으로 활용되고 있다.

격자볼츠만법이 비압축성 유동을 풀수 있는 근거가 마련된 것도 1990년대 후반이며(He et al.[1997]), 다상류 유동 중에서 격자볼츠만법을 이용한 자유수면 해석법은 비교적 최근에 연구된 분야이다(Thürey[2003]).

본 논문에서는 2차원 9방향(D2Q9) 방식으로 중력장하에서 물-공기 영역 중 액체 영역만을 고려하여 단상류로 자유수면을 시뮬레이션 하기위하여 개발된 코드를 소개한다. 실험 결과와의 비교를 통하여 격자볼츠만법의 활용가능성을 검증 하였다. 자유표면의 재구성방법, 분포함수의 조합으로 이루어진 경계조건, 중력장의 안정화, 격자의존성, 밀도함수의 정량평가, 댐 붕괴 시뮬레이션 및 경사면 시뮬레이션으로 구성되어 있다. 특히, 격자볼츠만법의 밀도함수에 대해 구체적인 평가를 시도하였다. 댐 붕괴 시뮬레이션에서는 그 결과치가 실험치의 데이터와 일치함을 보였다.

격자볼츠만법을 이용한 유체 시뮬레이션은 직교격자를 기반으로 사용되어져 왔다. 그러나 임의 형상에서 유체 시뮬레이션이 요구됨에 따라 유체-물체 경계면 처리에 대해서 연구되고 있다[Mei et al.(2002)]. 본 연구에서는 고전적 댐 붕괴시뮬레이션으로 검증단계를 거쳐, 자유수면 격자볼츠만법의 활용도를 넓히기 위해 경사면 형태의 물체가 존재할 경우 자유수면 시뮬레이션을 수행하였다.

마크로(macro)와 마이크로(micro)의 중간영역인 메조스케일(Meso scale)을 대상으로 하는 볼츠만방정식은 입자군의 위치, 속도, 시간의 함수인 분포함수 f(x,u,t)로 정의하고 있으며, 그 분포함수의 “충돌”과 “전진”등의 조합 및 분포함수의 합으로 속도장 및 압력장을 해석한다. Bhatnagar et al.[1954]가 볼츠만 기본방정식에 충돌모델을 추가함으로 격자볼츠만시뮬레이션이 유용한 툴로서의 역할을 가지게 되었으며, 동일면(phase)내 밀도의 변동이 무시할 수 있을 정도로 작다는 가정하에서 Chapman-Enskog전개를 거치게 되면 결국 비압축성 유동의 Navier-Stokes방정식으로도 변환되게 된다(He et al.[1997]).

입자의 분포를 함수화한 단일 밀도분포함수 f_i, 격자의 이동속도 e_i를 사용하여 2차원 9개의 방향(D2Q9)의 Single-Time-Relaxation Lattice-Boltzmann 방정식을 표현하면 식 (1)과 같다.

식 (1)은 i의 방향으로 전진하는 국소분포함수값 f_i의 전진단계(등호 왼쪽)와 충돌단계(등호 오른쪽)로 이루어져 있고, 다시 충돌단계에는 국소분포함수와 국소등가함수(f ^eq)로 구성되며, 충돌항에 완화계수(τ^-1)가 곱해져 있다. 시간스텝동안 입자는 속도벡터방향으로 가장 가까운 격자점까지 이동하게 된다. D2Q9의 이산화 방향은 그 격자점을 포함하여 총 9개(i=0~8)로 구성되어 있으며, 계산의 단순화를 위해서 시간스텝은 ∆t=1로, 공간스텝 ∆x=e_i∆t로 두었다. 또한 입자의 질량을 1로 설정하면 격자점 x, 시간 t에 있어 i방향의 운동량은 e_i f_i(x,t)이 된다. e_i는 8개방향의 격자단위속도를 나타낸다. 즉, i=0~8일 때 주방향의 속도성분과 (0,0), (1,0), (0,1), (-1,0), (0,-1), 대각방향의 속도성분 (1,1), (-1,1), (-1,-1), (1,-1)으로 구성된다(Fig. 1 참조). 본 논문에서는 비압축성유동을 다루고 있기 때문에 유체변수값에 의존하는 f_ieq식은 아래와 같다(He et al.[1997]).

Fig. 1 
Two dimensional D2Q9 model.

w_i는 계수로서 i가 가로세로 방향인 1~4일때는 1/9, 대각선방향인 5~8일때는 1/36이 쓰였다. 또한 0 일 때 국소등가함수 f₀^eq=4/9(ρ - 3/2u²)를 도입하였다. 비압축성유동의 범위(incompressible limit)로서 Mach수가 낮은 경우에 식 (2)가 사용가능하다(He et al.[1997]). 또한 식 (2)에 나타나 있는 마크로 스케일의 물리량이며, 관심 대상인 속도와 밀도는 다음과 같이 표현된다.

즉, 밀도는 밀도분포함수의 합산으로 얻을 수 있으며, 속도는 분포함수에 격자단위속도를 곱한 것을 합산하고 밀도를 나누어 계산할 수 있다.

단일완화시간모델(Single-time-relaxation approximation)의 운동량분포는 식 (1)에 따라 매 시간스텝의 “충돌”로 인하여 평행분포에 근접(완화) 한다는 방식이며, 완화계수 τ는 n 스텝 f_i^eq을 기준으로 해서 n+1 스텝에서의 f_i 값을 결정하는 역할을 한다. LBM에서 해의 안정성에 영향을 미치는 파라메터이기도 한 τ는 0.5이상의 값을 가져야만 유효성을 가진다. 따라서, τ의 값이 0.5 값에 가까워질수록 안정성에 문제를 야기 할 수 있으므로 주의해야 한다. 본 논문에서는 0.505 ≤ τ ≤ 0.953의 범위에서 시뮬레이션을 수행하였다.

본 논문에서 다루고 채택한 셀의 모양은 사각형으로서 중심에 속도와 밀도가 정의되어 있는 형태이다. 따라서 고체경계는 셀면에 위치하게 된다. Fig. 2는 직교 격자에 경사면이 있는 경우, bounce-back 방향의 밀도함수 처리방식을 설명하고 있으며, 여기서 점 P에서 방향의 8번이 진행상에 경계면이 있다고 가정할 경우를 예시하고 있다. 즉 bounce가 된 이후 6번방향의 새로운 밀도함수를 구하는 방식에 대한 내용이다. 여기서 W는 벽면, B는 고체면내, F는 P를 중심으로 6번방향의 점을 나타낸다. Fig. 2에서 나타낸바와 같이 2가지의 경우를 고려 할 수 있다. 즉, P점과 경계면과의 거리가 비교적 먼거리(A의 경우)와 가까운 경우(B의 경우)로 구분할 수 있을 것이다. Bouzidi, M. et al.[2001]은 P점과 B점사이에 W가 위치 할 때 PW의 거리가 PB의 거리와 비교 했을 때 1/2이상일 경우와 1/2미만일 때를 나누어서 반대 방향으로의 밀도 함수를 내삽하는 방식을 택하였다. 즉 로 나타내면, 각방향의 밀도 함수는 다음과 같다.

∆의 값이 1/2이라면 일반적인 bounce-back방식이 될 것이다. 여기서 는 유체가 벽면을 향하는 방향의 반대 방향을 의미한다. 또한 ∆가 1/2보다 작을 때 점 F의 값도 반영이 된다는 것을 알 수 있다. 식 (9)과 (10)의 좌변항의 밀도 함수는 충돌과 전진스텝을 모두 거친 후의 결과값이며, 우변항의 밀도함수는 충돌전의 값을 시뮬레이션에서 사용하였다.

Fig. 2 
Bounce back treatment with inclined wall.

또한, P점에서 8개의 방향마다 ∆값을 계산하는 알고리듬을 아래와 같이 나타내었다. 예를들어 Fig. 2(B)의 경우 방향1번, 4번과 8번이 해당이 될 것이며, 이 ∆값은 아래와 같이 대수학적으로 정리할 수 있다. 1번 방향은 (I×a)/(J−a), 4번 방향은 a, 8번 방향은 (a × b)/ 으로 계산할 수 있으며, 여기서 a= J−t, t=tanα, 그리고 b= sinβ/sin(135^o−β)이다. I와 J는 계산 코드상의 가로 및 세로방향의 계산배열을 나타내며, α,β는 하단과 상단의 각도를 의미한다.

본장에서는 4장의 경사면 처리방법에서 기술한 바에 따라 일반적 물체에도 유체시뮬레이션이 가능하도록 해당 알고리듬을 적용 하였으며, 경사면이 포함되어 있는 경우에 대해서 시뮬레이션을 수행하였다. Table 6에는 사용된 입력값을 나타내며. 격자수는 2000×500, ∆x는 5×10^-5m, ∆t는 0.5×10^-4s를 사용하였다. g'을 0.0005로 정한 이유는 τ값을 최대한 큰 수를 택하기 위해서이다. Fig. 10에는 댐붕괴와 같이 본 예제에서는 높이 12.5mm의 사각 물기둥이 경사면이 있는 장애물을 넘는 시뮬레이션을 수행하였다. 세 s가지 시간스텝을 보여 주고 있는데 t=0.25정도에서 장애물에 도달 후 이후 t=0.5일 때 경사면을 넘어가는 형태를 확인할 수 있다. 그 이후 수직벽에 부딪히면서 장애물 오른쪽에 유체가 담기는 형태를 볼 수 있다. 전체적인 흐름제현은 가능한 것으로 파악된다. 결과적으로 자유 수면 LBM기법을 이용하여 일반형상 주변의 시뮬레이션도 가능함을 보였다. 단, 자유수면 처리에 있어 보다 세밀한 처리가 요구되며, 추후 검증과정을 거쳐 보다 신뢰성 있는 코드개선이 요구된다.

Fig. 10 
Schemetic view for the freesurface simulation with inclined bed (units: mm).

Fig. 11 
Free surface movement under the gravitational field with inclined bed.

본 논문은 비교적 최근에 알려져 있는 Boltzmann 방정식에 기반을 둔 격자볼츠만법(Lattice Boltzmann Method)을 이용하여 자유수면 시뮬레이션 코드를 개발하였다. 격자볼츠만법은 Navire-Stokes방정식과 비교하여 Knudsen수의 적용 대상 범위가 넓은 특징이 있으므로 응용확장성이 높다. 본 연구에서는 자유수면 시뮬레이션시 격자볼츠만 고유의 특성을 비교적 상세히 고찰하였다. 본 논문은 지배방정식에 포함되어 있는 단일 밀도 함수로 인한 수면상에서의 재구성 방안, 수면의 이동 방식, 중력장내 밀도분포, 격자 의존성등에 대해서 논하였으며, 댐 붕괴 시뮬레이션 사례를 들어 계산 결과값들을 실험치등과 비교분석하였다.

특히 중력장에서의 무차원 중력세기 값인 g'의 특성을 조사하고, 댐 붕괴시뮬레이션에서 g'의 의존도도 고찰하였다. 또한 각 방향의 밀도 함수값을 정량적으로 파악하였으며, 중력세기를 적절히 선정하는 것이 요구된다. 격자볼츠만법의 시뮬레이션상의 격자 의존도도 확인하였으며, 코드의 검증을 위해 기존 댐 붕괴 시뮬레이션의 실험값 및 타 결과값들과 일치함을 보였다. 마지막으로 LBM의 응용성을 확장하기 위하여 일반형태의 장애물이 존재시에도 시뮬레이션이 가능하도록 코드를 확장하였으며, 추가적 검증단계를 거쳐 유용성 및 계산 정도에 대한 검증을 확보할 필요가 있다고 판단한다.