Study on the Free Surface Behavior Using the Lattice Boltzmann Method

1. 서론

공간과 시간에 따라 변화하는 유체현상을 시뮬레이션 할 경우, 관심영역(공간스케일, 시간스케일)을 정의하게 된다. 주로 유체 시뮬레이션시 활용되는 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)]. 본 연구에서는 고전적 댐 붕괴시뮬레이션으로 검증단계를 거쳐, 자유수면 격자볼츠만법의 활용도를 넓히기 위해 경사면 형태의 물체가 존재할 경우 자유수면 시뮬레이션을 수행하였다.

2. Lattice-Boltzmann Method(격자볼츠만법, LBM)

마크로(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의 범위에서 시뮬레이션을 수행하였다.

3. LBM 자유수면 시뮬레이션

4. 경사면 처리 방법

본 논문에서 다루고 채택한 셀의 모양은 사각형으로서 중심에 속도와 밀도가 정의되어 있는 형태이다. 따라서 고체경계는 셀면에 위치하게 된다. 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는 계산 코드상의 가로 및 세로방향의 계산배열을 나타내며, α,β는 하단과 상단의 각도를 의미한다.

5. 검증계산 및 고찰

5.1 격자볼츠만의 중력효과

중력장이 작용하는 자유수면 문제를 시뮬레이션하기 위해서는 중력항이 시뮬레이션 상에 작용하는 효과를 우선적으로 파악해야 한다. Fig. 3에는 정사각형 계산영역(2.5×10^-3m, 2.5×10^-3m)에 높이 H만큼 물이 차있다고 가정한다. 100×100의 정사각형 격자를 구성하였으며, 초기 정수상태에서의 물의 무차원 밀도(ρ'=ρ∆x³/∆m>)를 1로 가정하였다. pt1은 자유수면의 셀, pt2는 H만큼 떨어진 바닥위 셀을 나타낸다. 아래 Table 1에는 계산시 사용된 파라메터이다.

Fig. 3

Problem setup for the density strength along vertical direction.

Table 1

Parameters of simulation for gravitational effect

Fig. 4

Time history of density values at pt1 and pt2.

무차원 중력값인 g'은 g(∆t)²/∆x로 정의하였다. g는 중력가속도 값을 사용하였다. 시뮬레이션 과정에서 LBM은 중력의 효과를 안정화시키는 방향으로 수렴해 나간다. 이것은 Fig. 4에 나타낸것과 같이 pt1 점근방은 자유수면 인접점으로 중력항의 영향을 거의 받지 않으므로 무차원 밀도의 기본값인 1을 유지하고 있는 반면, 바닥부근의 pt2은 중력항의 영향을 가장 많이 받는 지점으로 초기의 불안정화 단계를 거쳐 안정화되어 가는 것을 알 수 있다.

밀도의 수직분포를 살펴보았다. 무차원 중력값 g'의 변화에 따른 밀도 구배를 Table 2의 값을 사용하여 시뮬레이션을 수행했다. g'값이 크면 밀도구배가 커진다는 것을 알 수 있으며 계산된 dρ/dH값 즉 O(10^-4)은 Buick et al.[2000] 결과값의 범위에 포함된다(Fig. 5). 또한 Fig. 5에서 해당 정수압력값(ρgH)에 대한 LBM의 계산결과 압력값은 g' 따라 약 2%~4.5%의 범위내에 있다는 것을 알 수 있다.

Table 2

Parameters of simulation for gravitational effect

Fig. 6는 계산영역상에 밀도강도(g'=0.0025)가 안정화 된 상태가 될 때까지 각 방향(0~8)에 따른 밀도 함수값을 나타낸 그림이다. 밀도값을 구성하는 요소중 f0의 값이 약 44%, 십자방향의 f1, f2, f3 그리고 f4가 약 10%씩, 대각방향의 f5, f6, f7 그리고 f8이 약 3%를 점유하는 것을 알 수 있다. 또한 중력의 효과로 인해 pt2에서의 밀도함수들의 값이 pt1보다 더 커진 것을 알 수 있으며, 그 밀도의 양은 약 3% 정도 차이가 난다는 것을 Fig. 5에서 파악할 수 있다. Table 3에는 Fig. 6의 막대그래프의 값을 나타낸 것으로 밀도의 평균값을 비교해보면 pt2의 값이 더 큰 것을 알 수 있다.

Fig. 5

Density gradient with three different gravity strengths.

Fig. 6

Values of density functions at certain point pt 1 and pt 2 with g'=0.0025.

5.2 LBM 댐 붕괴 시뮬레이션 사례 및 특성

자유수면 시뮬레이션 코드 개발시의 대표적 검증사례인 댐 붕괴 시물레이션을 수행하였다. 중력의 영향으로 유체가 하강하면서 이동하는 시뮬레이션을 수행하여 기존의 계산 및 실험값과 비교 검토하였다. 댐 자유수면 하단부 선단이 오른쪽으로 전진하는 전진거리를 기존 실험 문헌(Martin et al.[1952])과 비교 하였다. 또한 격자볼츠만법의 밀도함수 f의 특성에 대해서도 기술하였다.

Table 3

Density function break-down at pt1 and pt2 with =0.0025

Fig. 7

Comparison of the water front position by present methodwith those of experiments and other simulation results.

Table 4

Grid dependency test with four different sizes in the same g'

먼저 Fig. 7은 유체 하단부의 선단부의 진행거리를 무차원시간에 대하여 나타난 그래프이다. 다른 실험결과와 VOF시뮬레이션 결과와 함께 나타내었다. 격자 사이즈에 대한 해의 의존도를 파악하기 위해서 쓰인 격자수는 50×50, 80×80, 100×100과 200×200등이다(Table 4). 50×50의 경우를 제외하고는 실험치와 근접함과 동시에 격자의존성이 두드러지게 나타나지는 않는 양상을 보이며 수렴되는 경향을 타나낸다. ∆x값 5.0×10^-5 정도의 범위에서 수렴되었다고 판단된다. 본 수렴정보를 바탕으로 이후 시뮬레이션시 ∆x값을 참고 하였다. 무차원 ∆t의 범위는 3.57×10^-5s~7.14×10^-5s 이며, 완화계수 τ의 범위는 0.505~0.543이며, 무차원 격자중력값 g'은 0.00025로 동일하다.

또한 시뮬레이션 시 질량보존의 정도를 검증 한 결과 시뮬레이션 초기에 약 0.15%의 차이가 발생하였으나 계산이 진행되면서 0으로 가까워져 가는 경향을 보였다(Jung[2012]). 질량변화의 정의는 [(m_ini-m_cal]/m_ini]×100로 나타내었으며, 계산된 유체 전체질량(m_cal)과 초기 유체전체질량(m_ini)과의 비율을 나타낸 것이다.

Table 5는 무차원 격자중력값의 영향을 확인하기 위하여 시뮬레이션을 수행하였다. 즉, 검증된 격자크기 ∆x를 고정하고, g'값을 0.0001~0.0005의 범위내이며, τ는 0.5096~0.5214이다. g'값의 선택은 시뮬레이션 하고자 하는 대상에 따라 달라지며, 댐 붕괴와 관련해서는 0.00025를 택하였다. g' 값과 ∆t 값은 비압축성 계산시 연계가 되어있으므로 주의해서 선정해야 한다.

Table 5

Dependency of gravity strength with a specific Δx

Fig. 8

The effect of non-dimensional gravitation strength with Dam-Breaking simulation.

Fig. 8의 시뮬레이션 결과를 보면 작은 g' 값이 높은 g' 값보다 유체이동량에 있어 미세하게 진동이 있다. g' 이 0.0005일때는 유체선두가 선형적으로 전진하는 경향이 있다. 즉, 큰 g' 값은 시뮬레이션 전반에 걸쳐 안정화에 기여하는 것으로 판단한다.

Fig. 9는 100×100의 경우 셀 위치 [50,5]에서의 밀도 함수 자체값(f₁~f₈)의 값과 속도 u와 v를 시간 스텝에 대해 나타낸 것이다. 수평방향의 속도성분이 지배적인 위치(Fig. 7의 사각내 별표위치)를 선정했다. 여기서, u에 관여하는 f 값은 격자속도성분 주축성분중 수평성분 1번과 3번, 그리고 대각방향으로 5,6,7과 8번이다. 그 값을 살펴보면 주축성분은 평균 0.11정도의 값을 가지며, 대각방향의 성분(중간 표기)은 약 0.03의 값을 가지고 있다. Fig. 9에서 보듯이 f 값은 정과 골로 이루어지는 형상을 보이고 있다. u의 속도 형태를 보면 첫 번째 정은 사각형의 유체가 아래로 내려오면서 큰 모멘텀을 받았기 때문이고 두 번째 정은 첫 번째 보다는 작지만 어느 정도의 유체량이 흘러간 것을 의미하며 그 다음 골은 오른쪽 벽면을 터치하고 돌아오는 영역이고 다음의 정은 왼쪽벽면을 치고 오른쪽으로 유체가 흐르고 있다는 것을 짐작할 수 있다. 특히 f₁과 f₃은 위상이 반대임을 비교적 명확히 알 수 있다. 이와 같이, 시뮬레이션중에 모든 f 값은 그들 고유의 방향성분 값을 가지며, 결국 속도 값으로 잘 표현되고 있음을 알 수 있다.

Fig. 9

Time variation with density functions and velocity at a certain point [50,5] (Fig. 7의 사각내 별표위치).

6. LBM 경사면 시뮬레이션

본장에서는 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기법을 이용하여 일반형상 주변의 시뮬레이션도 가능함을 보였다. 단, 자유수면 처리에 있어 보다 세밀한 처리가 요구되며, 추후 검증과정을 거쳐 보다 신뢰성 있는 코드개선이 요구된다.

Table 6

Parameters values for the test simulation of inclined walls

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.

7. 결 론

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

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

Acknowledgments

이논문은 울산대학교 교내연구비 지원으로 수행되었습니다.

References

• P.L Bhatnagar, E.P Gross, M Krook, “A model for collision processes in gases. I : small amplitude processes in charged and neutral one-component system”, Phys. Rev, (1954), 94, p511-525. [https://doi.org/10.1103/PhysRev.94.511]
• M Bouzidi, M Firdaouss, P Lallemand, “Momentum transfer of a Boltzmann-lattice fluid with boundaries”, Phys. Fluids, (2001), 13, p3452-3459. [https://doi.org/10.1063/1.1399290]
• J.M Buick, C. A Greated, “Gravity in a lattice Boltzmann model”, Physical Review E, (2000), 61(5), p5307-5320. [https://doi.org/10.1103/PhysRevE.61.5307]
• H Chen, S.A Orszag, I Staroselsky, “Macroscopic description of arbitrary Knudsen number flow using Boltzmannn-BGK kinetic theory”, J. Fluid Mech, (2007), 574, p495-505. [https://doi.org/10.1017/S0022112006004241]
• X He, L.-S Luo, “Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation”, Journal of Statistical Physics, (1997), 88(Nos.3/4), p927-944. [https://doi.org/10.1023/B:JOSS.0000015179.12689.e4]
• C.W Hirt, B.D Nichols, “Volume of fluid (VOF) method for for the dynamics of free boundaries”, Journal of Computational Physics, (1981), 93, p201-225. [https://doi.org/10.1016/0021-9991(81)90145-5]
• Rho-Taek Jung, “Feasibility study on the two-dimensional free surface simulation using the Lattice-Boltzmann Method”, J. of the Korean Society for Marine Environmental Engineering, (2012), 15(4), p273-280. [https://doi.org/10.7846/JKOSMEE.2012.15.4.273]
• Y Kato, K Kono, T Seta, D Martinez, S Chen, “AMADEUS Project and Microscopic Simulation of Boiling Two-Phase Flow by the Lattice-Boltzmann Method”, International J. of Modern Physics C, (1997), 8, p843-858. [https://doi.org/10.1142/S0129183197000722]
• M.-S Kim, F Perot, M Meskine, “Aerodynamics and acoustics predictions of the 2-blade NREL wind turbin using Lattice Boltzmann Method”, 14th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, ISROMAC-14, Honolulu, USA, (2012).
• C Korner, M Thies, T Hofmann, N Thürey, U Rude, “Lattice Boltzmann Model for Free Surface Flow for Modeling Foaming”, Journal of Statistical Physics, (2005), 121(1/2), p179-196. [https://doi.org/10.1007/s10955-005-8879-8]
• G Lu, D.J Depaolo, Q Kang, D Zhang, “Lattice Boltzmann simulation of snow crystal growth in clouds”, J. of Geophysical Research, (2009), 114(D07305). [https://doi.org/10.1029/2008JD011087]
• J.C Martin, W.J Moyce, “Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane”, Phil. Trans. R. Soc. Lond. A, (1952), 244, p312-324. [https://doi.org/10.1098/rsta.1952.0006]
• R Mei, D Yu, Wei Shyy, “Force evaluation in the lattice Boltzmann method involving curved geometry”, Physical Review E, (2002), 65(041203). [https://doi.org/10.1103/PhysRevE.65.041203]
• E Orlandini, M.R Swift, J.M Yeomans, “A lattice Boltzmann model of binary-fluid mixtures”, Europhysics Leters, (1995), 32(6), p463-468. [https://doi.org/10.1209/0295-5075/32/6/001]
• N Thürey, “A single-phase free-surface Lattice Boltzmann Method”, Master thesis, INSTITUT FÄUR INFORMATIK (MATHEMATISCHE MASCHINEN UND DATENVERARBEITUNG), (2003).
• B Wen, H Li, C Zhang, H Fang, “Lattice-typedependent momentum-exchange method for moving boundaries”, Physical Review E 85, (2012), 016704. [https://doi.org/10.1103/PhysRevE.85.016704]