Journal of the Korean Society for Marine Environment and Energy
[ Article ]
Journal of the Korean Society for Marine Environment and Energy - Vol. 15, No. 4, pp.273-280
ISSN: 2288-0089 (Print) 2288-081X (Online)
Print publication date Nov 2012
Received 29 Feb 2012 Revised 09 Aug 2012 Accepted 06 Sep 2012
DOI: https://doi.org/10.7846/JKOSMEE.2012.15.4.273

Feasibility Study on the Two-dimensional Free Surface Simulation Using the Lattice-Boltzmann Method

JungRho-Taek
University of Ulsan
Lattice Boltzmann Method를 이용한 2차원 자유수면 시뮬레이션 기법연구

Correspondence to: rtjung@ulsan.ac.kr

The numerical simulation using the Lattice Boltzmann Method in the field of computational fluid dynamics becomes wider in the engineering applications because of its simplicity of update rules compared to the conventional Navier-Stokes solvers. Here, a two-dimensional D2Q9 LB model is numerically tested with a few new computational treatment on the free surface. The single relaxation time is applied under the gravitational field where applied only in the higher density fluid because of its big density difference. At the free surface, the reconstruction techniques in combination with boundary conditions is adopted in order to get some distribution function coming into the fluid site from the air one, and surface tension, early stable test for the gravitional field is considered in it. With the implementation of the gravitational profile, conserving the overall mass and grid dependency are observed during the calculations and freesurface advance track is well captured with an experiment.

초록

전산유체역학의 격자볼츠만법은 Navier-Stokes방정식의 시뮬레이션 기법 보다 비교적 간략한 이산화 방식으로 인하여 공학적 응용성을 더욱 넓혀 가고 있다. 본 논문에서는 단일 완화계수 및 D2Q9 방식으로 중력장하에서 액체영역에서의 다이나믹스만 계산하는 자유수면 시뮬레이션을 수행하였으며, 그 활용성을 검증하였다. 자유표면의 재구성방법, 분포 함수의 조합으로 이루어진 경계조건, 표면장력, 중력장의 안정화, 격자의존성, 자유수면 끝단의 하단 벽면 이동 검증 등을 수행하였으며, 그 결과치가 실험치의 데이터와 일치함을 보였다.

Keywords:

Lattice-Boltzmann Method, Dam-breaking Simulation, D2Q9, 격자볼츠만방법, 댐붕괴 시뮬레이션, 차원 9방향

1. 서 론

공기와 물, 이상(two-phase)간의 자유수면 시뮬레이션은 중요한 공학적 관심 대상이다. 자유수면 시뮬레이션을 위해 여러 연구자들이 다양한 알고리즘을 개발해 왔으며, 두 상간에 급격한 밀도 차이(O(3))로 인한 자유수면 경계조건의 만족과 그 표현방법에 대해서도 여러 연구사례가 있다(Hirt et al.[1981], Sussman et al.[1994],Tryggvason et al., [2011], Yabe [2001]). 특히 물이나 연기를 중심으로 한 유체시뮬레이션은 공학자 뿐만 아니라 일반 대중들의 관심을 이끌기에 충분하므로 컴퓨터 그래픽분야에서도 다양한 시도들이 이루어지고 있다(Thurey et al. [2010], Wolfram et al. [2008],Xiaoming et al. [2003]).

공기-물의 자유수면 시뮬레이션은 공기와 물을 따로 시뮬레이션 하는 방법과 동점성계수가 공기보다 작은 물만 시뮬레이션을 하는 경우가 있다. 시뮬레이션의 목적에 따라 공기부분도 동시에 시뮬레이션 해야 하는 경우가 있을 것이나, 주로 고체면에 주어지는 유체력 추정을 대상으로 할 경우 공기보다는 물만을 대상으로 하는 것이 계산시간 및 계산메모리 할당면에서 경제적일 것이다.

한편 본 논문에서는 Navier-Stokes(N-S)방정식을 직접 이용하지 않고 유체입자를 그 위치에서의 분포함수(통계함수의 결과값)로 정의하고, 그 분포함수의 "충돌"과 "전진" 등의 조합에 의해서 속도장을 구하는 Lattice Boltzmann 방법(Lattice Boltzmann Method: LBM)을 본 논문에 도입하여 자유수면 시뮬레이션 프로그램을 개발하였다.

LBM은 입자 분포함수의 시간변화율을 나타내는 Boltzmann transport 방정식에 기반을 두고 있으며, 그 변화율은 해당 격자점로 들어오는 입자 수에 나가는 입자수의 차로 표현하게 된다. 그리고 동일면(phase)내 밀도의 변동이 무시할 수 있을 정도로 작다는 가정하에서 Chapman-Enskog전개를 거치게 되면 결국 비압축성 유동의 N-S방정식으로 변환되게 된다(He et al. [1997]). LBM은 해석 입자관점에서 보면, Continuum Mechanics와 Molecular Dynamics 의 중간 영역에서 출발되었다고 볼 수 있다. N-S방정식 또는 비선 형미분방정식은 시간 및 공간적 평균을 취하는데 비해, LBM은 동역학적 다이나믹스 문제를 직접 다룰 수 있다는 점과 연속체 개념에서 중요시되는 물리량인 밀도와 속도를 입자분포함수의 간단한 조합으로 구현할 수 있는 특징이 있다.

1980년대 개발초기에는 LBM의 공학적 응용성이 낮았으나 1990 년대에 접어들면서 다양한 분야에서 LBM이 이용되고 있다. LBM의 응용사례를 살펴보면, 원자력 플랜트내의 방사능 분열 반응시의 복잡한 다상현상을 실험적으로 구현하기 어렵기 때문에 Kato[1997] 등은 AMADEUS(Advanced Microscopic Analysis by Discrete and Emergent Computing Schemes for Complex Phenomena in Nuclear plants) 프로젝트를 통하여 LB기법을 이용해 기포의 생성과 성장을 표현하는 비등(沸)이상류 시뮬레이션을 미시적(microscopic)으로 제현하였다. 이상류에서 나타나는 접면마찰과 다른 상(phase)간에 에너지의 교환등 면을 통한 다이나믹스 문제를 본질적으로 다룰 수 있다는 것이다. 뿐만아니라, 고속도로 차량 흐름 해석(Resnics[1995], Nagel [1997]) 등과 같은 large scale 해석과 생체면역시스템 해석(Chowdhury [1990]), 다공물질 흐름(Buckels [1994]), 기포거동(Nadiga [1996]), 상분리 (Shan [1993])등의 연구가 진행되었으며, 유체진동해석과 고 레이놀즈수 시뮬레이션, 메탈 폼해석(Korner et al. [2005])에도 응용되고 있다.

1980년대 개발초기에는 LBM의 공학적 응용성이 낮았으나 1990년대에 접어들면서 다양한 분야에서 LBM이 이용되고 있다. LBM의 응용사례를 살펴보면, 원자력 플랜트내의 방사능 분열 반응시의 복잡한 다상현상을 실험적으로 구현하기 어렵기 때문에 Kato[1997] 등은 AMADEUS(Advanced Microscopic Analysis by Discrete and Emergent Computing Schemes for Complex Phenomena in Nuclear plants) 프로젝트를 통하여 LB기법을 이용해 기포의 생성과 성장을 표현하는 비등(沸)이상류 시뮬레이션을 미시적(microscopic)으로 제현하였다. 이상류에서 나타나는 접면마찰과 다른 상(phase) 간에 에너지의 교환등 면을 통한 다이나믹스 문제를 본질적으로 다룰 수 있다는 것이다. 뿐만아니라, 고속도로 차량 흐름 해석(Resnics[1995], Nagel [1997]) 등과 같은 large scale 해석과 생체면역시스템 해석(Chowdhury [1990]), 다공물질 흐름(Buckels [1994]), 기포 거동(Nadiga [1996]), 상분리 (Shan [1993])등의 연구가 진행되었으며, 유체진동해석과 고 레이놀즈수 시뮬레이션, 메탈 폼해석(Korner et al. [2005])에도 응용되고 있다.

비압축성문제를 해석시에 통상의 CFD는 압력해를 도출하기 위해서 포아슨 방정식(Poisson Equation)을 푼다든지 ACM(Artficial Compressibility Method)을 이용하지만 LBM의 경우는 주어진 문제에 따라 상태방정식에서 압력을 구하거나 운동량 변화(momentum exchange)방식에 의해서 구할 수 있다(Jung [2009], Huabing [2004]).

본 논문은 공기-물을 대상으로 LBM기법을 도입하여 유체시뮬레이션에 적용하였다. 이상류를 표현하기 위해 각 상별로 국소함수를 사용하는 기법(Orlandini et al. [1995])이 있으나 본 논문에서 다루고자 하는 공기-물은 밀도비가 크므로 물만을 대상으로 한 단일 국소함수 형태로 자유수면을 시뮬레이션 하였다. 본 논문은 지배방정식에 포함되어 있는 완화계수의 도출방법, 단일 국소 함수로 인한 수면상에서의 재구성 방안, 수면의 이동 방식, 중력장내 밀도 분포, 격자의존성 등에 대해서 논하였으며, Dam Break의 사례를 들어 실험치와 결과값들을 실험치와 비교분석하였다.


2. 자유수면 해석을 위한 격자볼츠만 방법 (Lattice-Boltzmann Method, LBM)

2.1 지배방정식

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

식 (1)은 i의 방향으로 전진하는 국소분포함수값 fi의 전진단계(좌항)와 충돌단계(우항)로 이루어져 있다. 이 충돌단계의 구성은 국소분포함수와 국소등가함수로 구성되며, 충돌항에 완화계수 가 곱해져 있는 기본적인 LBM식이다. 시간스텝동안 입자는 속도 벡터방향으로 가장 가까운 격자점까지 이동하게 된다. 이산화의 방향은 그 격자점을 포함하여 총 9개(i=0~8)로 구성되어 있으며, 계산의 단순화를 위해서 시간스텝은 1로, 공간스텝 로 두었다. 또 한 입자의 질량을 1로 설정하면 격자점 x, 시간 t에 있어 i방향의 운동량은 된다. ei는 8개방향의 격자단위속도를 나타낸다. 즉, i=0~8일 때 속도성분은 (0,0), (1,0), (0,1), (-1,0), (0,-1), (1,1), (-1,1), (-1,-1), (1,-1)으로 구성된다. 식 (1)을 N-S방정식으로 변환이 가능하며(Bhatnagar [1986]), 이때 속도에 대한 이차 식의 형태로 fieq 식이 표현된다. 본 논문에서는 비압축성유동을 다루고 있기 때문에 유체변수값에 의존하는 fieq 식은 아래와 같다(Heet al. [1997]).

wi는 계수로서 i가 가로세로 방향인 1~4일때는 1/9, 대각선방향인 5~8일때는 1/36이 쓰였다(Fig. 2 참조). 또한 0 일 때 극소등가 함수는 를 도입하였다. 비압축성유동의 범위(imcompressible limit)로서 Mach수가 낮은 경우에 식 (2)가 사용가능하며, Chapman-Enskog 가정을 적용하면 LBM은 비압축성 N-S방정식으로 유도가 가능하다. 단, 주 대상 유체가 압축성 유체 또는 상분리(phase separation)와 같은 압축성을 해석코자 할 경우 feq내의 밀도는 [·] 밖으로 나와야만 한다(Orlandini[1995]). 또한 식(2)에 나타나 있는 물리량 속도와 밀도는 다음과 같이 표현된다.

Fig. 1

Momentum distributions in two dimensional rectangular site

Fig. 2

The single-time relaxation process(a) Particle momentum that will collide at the site P(b) Equilibulium momentum distribution when τ=1(c) Outgoing momentum distribution when τ=0.5

즉, 밀도는 분포함수의 합산으로 얻을 수 있으며, 속도는 분포함수에 격자단위속도를 곱한 것을 합산하여 밀도를 나누어 계산할 수 있다. Fig. 2에 본 논문에 사용된 차분의 방향과 특정시간의 p점에서의 운동량 분포를 나타내었다.

2.2 유체중력장내에서 단일시간 완화계수 τ

단일완화시간모델(Single-time-relaxation approximation)의 운동량분포는 식 (1)에 따라 매 시간스텝의 충돌로 인하여 평행분포에 근접(완화) 한다는 방식이다. 예를들어, 완화시간 τ의 값이 1인 경우 식 (1)에 의해 충돌후의 운동량 분포는 가 되 Fig. 2의 중심 p 점으로 들어오는 모든 방향의 운동량은 평형분포가 되어 나가게 된다(Fig. 2(a)에서 (b)로의 진행). 만약 τ의 값이 0.5인 경우에는 운동량은 평형 분포값에 두배에서 충돌전의 값을 뺀 형식이 될 것이다(Fig.2(a)에서 (c)로의 진행). 중력과 같은 외력항을 다루고 있을때는 외력의 방향으로 운동량이 커지게 된다. 따라서, 완화계수 τ는 n 스텝 fieq을 기준으로 해서 n+1 스텝에서의 fi값 M을 결정하는 역할을 한다.

본 논문에서 다루고 있는 τ의 결정방법을 설명한다. 우선 셀 속도 ucell을 기본으로 한 CFL 조건과 무차원된 격자중력항(g')으로 이루어진 두 항을 비교해서 최소값인 Δt를 식 (4)와 같이 구하는 것이 기본 개념이다. 본 논문에서는 CFL 값은 0.8로 두었으며, 격자 중력항은 비압축 시뮬레이션임에 따라 0.0025이하로 제한 하였다. 계속해서 무차원화된 격자동점성계수(v')를 구하게 되면, 단일 완화계수 τ는 식 (5)에 의해서 구할 수 있게 된다. 주어진 문제에서 유체의 점성계수가 주어지면 식 (5)에 따라서 시뮬레이션에 필요한 완화시간을 계산할 수 있다. 그러나 LBM에 사용되는 τ는 식에서 알 수 있는 것과 같이 LBM에서 해의 안정성에 영향을 미치는 파라메터이기도 한 τ는 0.5이상의 값을 가져야만 유효성을 가진다. 따라서, τ의 값이 0.5 값에 가까워질수록 안정성에 문제를 야기 할 수 있으므로 주의해야 한다. 본 논문에서는 0.548≤τ≤0.596의 범위에서 시뮬레이션을 수행하였다.

2.3 자유수면에서의 분포함수 재구성

격자볼츠만방법을 이용하여 기-액 이상유체는 기체와 액체 두 개의 상을 각각 두고 격자볼츠만방정식을 정의하여 시뮬레이션하는 방법이 있으나 본 논문에서는 서술한 바와 같이 단일시간 완화개념의 격자볼츠만법을 이용하면서 경계를 포함한 액체의 경우만 시뮬레이션하는 방법을 택하였다. 액체의 동점성계수가 기체보다 O(10-2)가 작으며 액체의 경우가 레이놀즈수가 크므로 액체가 기체로부터 받는 영향은 거의 없다고 가정한다. 액체만의 단일방정식으로 해석을 할 경우 자유수면에서 공기측에서 t+1 스텝의 분포함수가 요구되는 데, 이를 위해서 아래식으로 표현하였다. 기본적인 개념은 Thurey[2003]의 논문과 동일하며, 본 장에 간략히 정리하였다. 식(1)에 τ의 값을 1, 그리고 inn으로 들어오는 분포함수의 정의 점을 자유수면 셀로 가정하게 되면 아래와 같은 식이 도출된다.

아래첨자 'inn'은 t+1스텝에서 공기측에서 내부로 들어오는 방향이 되고 'out'은 액체에서 공기측으로 향하는 방향을 나타낸다. 그리고 위첨자 'eq'는 국소등가(local equilibrium)를 나타내며, 국소등 가함수 feq는 자유수면에서의 속도와 밀도의 함수로 표현된다. 여기서 finneq 내의공기에 대한 무차원 밀도 ρ는 1의 값을 가진다. 왜냐하면 본 방법의 기본적인 개념은 공기영역은 계산에 참여하지 않으며, 자유수면셀 내의 국소등가 함수의 정의 점을 액체에 둔다는 가정이 포함되기 때문이다. 따라서 공기측에서 들어오는 분포함수는 자유수면상에서 방향이 서로 반대인 등가함수의 합과 이에 상응하는 유체의 운동으로 결정된다. 식 (6)은 2-4에서 다룰 자유수면의 이동을 위해 분포함수의 재구성에 관련한 식이라고 볼 수 있다.

2.4 자유수면의 표현과 이동

자유수면의 이동방식은 Front Tracking Method (Hirt et al. [1981],Sussman et al. [1994]등)와 Front Capturing Method (Ryskin et al. [1984], Takagi et al. [1994]등)가 있다. LBM을 이용한 자유수면의 표현은 Front Capturing Method와 유사한 방식으로 국소분포함수의 조합으로 구현할 수 있다. 즉, 충돌단계 이후의 전진단계에서 유체 스칼라양을 이동 시킨다는 개념이다(Korner et al. [2005]. 우선, 계산 격자상에서 움직이는 유체 스칼라를 ε으로 표현하면

이 된다. 여기서, 은 fluid fraction, m은 질량, ρ는 밀도, vol은 셀 부피를 나타낸다. 자유수면의 이동은 질량의 시간변화로 볼 수 있으며, 한 셀에 있어 질량의 변화량은 주위 분포함수의 관계로 표현할 수 있다. 시뮬레이션 상에 인접 셀간에 질량의 변화는 두가지의 형태로 나눌 수 있다. 첫 번째는 유체셀(100% 유체로 이루어진 셀)과 자유수면셀, 두 번째는 자유수면 셀간이다. 유체셀 상호간의 질량변화는 없는 것으로 가정한다. 분포함수로 이루어진 질량의 변화는

첨자 i와 반대되는 방향으로 식 (8)의 등호 오른쪽의 첫 번째 항은 주위의 셀에서 현재의 셀로 들어오는 방향이 되고 두번째 항은 그 주위의 셀로 나가는 방향의 국소함소임을 알 수 있다. 이 질량변화는 유체셀과 자유수면 셀에서의 질량변화를 나타내는 식이고, 자유수면간의 질량변화는 해당 셀간의 스칼라양의 변화도 포함된 일반화된 식으로 표현할 수 있다.

식 (8)과 비교하면 식 (9)에는 유체의 스칼라량을 표현하는 ε을 포함하고 있다. 즉, 주위 셀과 이산화 되는 방향을 따라 질량교환이 이루어진다. 주변의 셀로 빠져나가는 양과 들어오는 양은 같다는 질량보존법칙 따라서 최종적으로 해당시간 스텝에서의 새로운 질량은 과거의 질량에 각 방향에서 주어진 질량의 변화를 더해서 얻게 된다.

2.5 외력항

격자볼츠만법의 중력항의 취급을 위해 다양한 방식이 제시되었다(Buick et al. [2000], Wen et al. [2012]등). 비교적 간단한 방법으로서는 충돌단계에서 수정된 속도를 추가적으로 구해서 그 속도 성분을 국소등가함수내에 사용한다는 것이다. 체적력으로서 작용하는 중력의 효과를 운동량의 변화 즉, 으로 해석하여 매 스텝마다 식 (11)을 통해서 계산된다.

여기서, 는 식 (3)으로부터 구한 값이며, 는 중력의 효과가 추가된 속도이며 이 값이 국소등가함수에 사용된다. τ는 완화계수이며, wi는 식 (2)에 쓰인 값과 동일하다. 는 유체력으로 이며 는 중력가속도이다. 실제 계산에서는 무차원된 값을 사용하게 된다. 이 값은 또한 incompressible limit에 영향을 끼치는 값이다. 4,7,8 방향으로는 양의 값이 2,5,6에는 음의 값이 적용되며, 1,3방향으로는 중력가속도가 작용하지 않는다.

2.6 표면장력

격자볼츠만방법에 있어 표면장력의 표현은 2-3의 재 구성항에 표면장력의 효과를 추가하는 형식으로 표현할 수 있다. 즉

여기서 첨자 inn은 유체외부에서 내부로 향하는 방향에 해당하고, out은 내부에서 외부로 향하는 방향을 의미한다. σ는 표면장력을 그리고 k는 곡률을 나타내며, 구하는 식을 아래 식 (13)과 식 (14)에 각각 나타내었으며, VOF(Voume of Fluid)방식에 쓰이는 알고리즘을 채택하였다.

ε은 스칼라량으로 식 (7)과 동일하며, σ'은 표면장력의 무차원 수이며, 은 자유수면 셀의 normal 방향성분을 나타낸다.

2.7 경계조건

본 논문에서 채택한 셀의 모양은 사각형으로서 control 셀은 그 사각형이 되며, 중심에 속도와 밀도가 정의되어 있는 형태이다. 따라서 고체경계를 생각할 때 셀중심과 셀면사이의 1/2위치에 존재한다. 따라서 2차정도의 벽면경계조건을 적용하였다(Kim [2000], Chen et al. [1996]). 예를들어 D2Q9의 제2방향으로의 밀도함수는

으로 표현할 수 있게 된다. 왜냐하면 등가밀도함수내의 셀내의 속도와 밀도를 적극적으로 이용하기 위함이다.


3. 계산사례

3.1 중력세기

본 장에서는 식 (11)에 포함되어 있는 중력항이 LBM시뮬레이션 상에 어떻게 작용하는지 알아보았다. 초기 밀도는 1.0으로 주어졌으며, 등가상태가 되었을 때 밀도의 수직분포를 확인하기 위함이다. Fig. 3과 같은 계산영역에 100×100의 정사각형 격자를 구성하였으며, 자유수면을 가지고 있는 정수상태를 가정하고, Δx는 2.5×10-5m, Δt는 2.0×10-5sec으로 하여, 완화계수 τ는 0.596으로 설정되었으며, 무차원 중력가속도인 g'은 0.0025를 두었다.

Fig. 3

Problem setup for the density gradient

문제의 설정으로는 사각경 계면은 no-slip 경계조건으로 두고, Fig. 3과 같이 b만큼의 유체가 채워져 있을시에, 무차원시간이 경과함에 따라 밀도의 수직 분포의 확인해 보았다. Fig. 4(a)는 밀도의 연직분포를 시간에 따라 나타낸 것을 시간이 지날수록 직선에 가깝게 된다는 것을 알 수 있으며, 유체하면에서의 밀도는 최상부보다 약 3%의 밀도가 증가 되었음을 알 수 있다. g'은 중력가속도의 세기를 나타내는 상수로서 그 값이 클수록 밀도의 수직 기울기는 더욱 증가 할 것이다. 본 계산에서 예를 든 것은 dρ/db=(1.03-1.0)/50=6.0×10-4으로 이 밀도의 기울기 값은 Buick et al., [2000]의 정도(order)내에 있다. 따라서 계산 시간이 지남에 따라 연직 밀도분포가 안정화 되어 가고 있는 것을 Fig. 4(a) 에서 알 수 있다. 따라서 3-2에 댐붕괴 해석시에도 동일한 g'값을 사용하였다. Fig. 4(a)에서 화살표시는 무차원시간이 0.062일때의 안정된 밀도 분포를 지시하고 있다. 또한 Fig. 4(b)에서는 Fig. 4(a) 의 *표시(50×1) 즉, 최하단면의 밀도변화를 무차원 시간으로 나타낸 것으로 진동을 하면서 밀도의 분포값이 매우 짧은 시간(0.05) 내에 중력장이 형성됨을 알 수 있다.

Fig. 4

(a) The density as function of non-dimensional depth for a fluid with ρ=1.0 and τ=0.596. The arrow is the density profile for the steady state.(b) The density variation of a point at the bottom layer as the time is evolved.

Fig. 5

Schematic view of typical Dam-breaking

3.2 Dam breaking 적용사례

본 논문에서는 격자볼츠만법을 이용하여 자유수면 시뮬레이션을 수행하였다. 개발된 코드의 계산사례로 2차원 Dam-breaking문제를 다루었다. Fig. 5에는 계산영역을 나타내고 있으며, 유체의 초기 형태를 보여주고 있다. 또한 Table 1에는 계산시에 사용한 물리적 파라메터를 나타내었으며, Δx=2.5×10-5m를 사용하였다. 식 (4)에 의해서 Δt는 0.56×10-5s로 주어졌으며, τ=0.596를 사용하였다.

먼저 격자수에 따라 질량보존의 정도를 검증하고자 한다. 자유수면 이동은 2-4장에도 언급하였듯이 자유수면셀과 자유수면셀, 자유수면셀과 유체셀의 물질이동에 의해 표현된다. 이산화 방향을 따라 밀도함수와 유체스칼라의 조합으로 주위의 셀과 질량의 차감을 반복하며 전체 질량을 보존토록 한다.

Parameters values for the test simulation

Fig. 6

Mass conservation check in case of 100×100

Parameters values for the grid dependency test

Fig. 6는 세가지의 다른 격자수 50×50, 80×80, 100×100를 택하여 질량의 변화를 비교한 것으로 [(mini-mcal)/mini]×100 계산된 유체 전체질량과 초기 유체전체 질량과의 비율을 나타낸 것이다. 100×100의 경우에 질량 변화를 나타내었으며, 격자가 작은 두 경우도 비슷한 경향임을 확인하였다. 계산영역의 한변의 길이를 0.0025 m로 주어졌을 때에 세가지 격자수 50×50, 80×80, 100×100의 경우를 테스트해 보았다. Table 2에 각각의 파라메터를 나타내었으며, g'의 값은 동일하며, Δx가 작아질 수록 τ값이 커짐을 알 수 있다. 즉 시뮬레이션이 안정적인 경향으로 주어진다는 의미이다. Fig. 6에는 격자수에 따른 질량의 보존을 시뮬레이션 한 경우이다. 댐의 유체가 반대편 벽면에 부딪치고 반사되는 시간영역에서 최대 ±0.2%의 미소한 변화가 존재하나 그 이후에는 질량보존이 복원됨을 알 수 있다.

Dam-breaking에 관련한 기존 실험 문헌(Martin et al. [1952])과 비교를 수행하였다. 주로 비교의 대상으로는 댐이 붕괴되면서 댐 자유수면의 선단 또는 상단이 전진 또는 하강하는 위치를 시간대별로 비교하였다. 따라서 Fig. 7에 댐 자유수면의 하단부 선단이 전진되는 변화양상을 비교한 그래프를 나타내었다. 격자수 [100×100]일 때 실험값과 가장 비슷한 경향인 것을 알 수 있다.

Fig. 8에는 격자수 [100×100]일때의 무차원 시간별 유체의 속도 장을 나타내었다. 유체 속도장의 형태를 살펴보면, Fig. 8(b) 우측 아래면의 코너부위의 역류존재, 상대벽면을 부딪치고 난 뒤의 흐름이 잘 표현되고 있는 것으로 사료된다. 단, 벽면경계조건인 no-slip 에 따라 잔류유체가 벽면에 남아있는 점, 유체가 벽에 부딪치고 난 뒤 스플래쉬 된 유체 입자의 물리적인 처리가 보완되어야 하는 점이 남아 있다.

Fig. 7

Comparison of the water front position with experiments and other simulation results

Fig. 8

Free surface movement under the gravitational field


4. 결 론

본 논문은 격자볼츠만 방법(Lattice Boltzmann Method)을 이용하여 자유수면 시뮬레이션 코드를 개발하였다. LBM 방식은 Navire-Stokes방정식과 같이 대류항이 존재하지 않음으로 수치확산과 같은 부차적인 문제점이 없고, 이산화 단계가 단순함에 따라 비교적 개발이 용이하다는 점이다. 그러나 계산영역을 고려해 볼 때 Navier-Stokes Solver의 연속체적 접근방식과 비교하면 입자의 통계학적 접근방식을 취하고 있는 LBM은 Mezo 또는 Micro 스케일의 계산 영역에 적합하다고 하겠다. 그럼에도 불구하고 비교적 전산유체역학의 새로운 시뮬레이션 방식인 LBM을 사용하여 자유수면 시뮬레이션의 가능성을 연구해 보았다.

본 논문에서는 기본적인 D2Q9(2차원9개방향)의 시스템을 통하여 비압축성유동을 위한 국소등가함수가 사용되었으며, 이 등가함수내에 외력항으로서의 중력항이 포함되어 있다. 중력장과 CFL조건으로 인한 완화계수의 결정방법을 논하였다. 벽면 조건으로는 2차정도의 bounce-back 방식을 적용하였고, 스칼라로 주어지는 수면의 이동을 위해 질량 이동 방법론에 대해서도 거론하였다. 질량보존이 만족함도 시뮬레이션을 통해 확인하였다. 자유수면 계산시의 중력장의 수직분포를 위해 짧은 계산시간 내에 밀도의 수직분포가 선형적으로 안정적으로 분포됨을 확인하였다. 또한 LBM은 격자크기 의존성이 나타남을 확인하였으며, 또한 이것이 계산 영역확장을 위해 메모리의 확장과 해석시간의 단축을 위해 병렬화가 필수 불가결하다는 결론이다. 따라서 통상적인 물리적 영역을 가정하여 LBM의 결과값을 NS Solver의 값과 비교하기 위해서는 상당히 큰 메모리의 확보가 필수불가결하다. 뿐만아니라, 보다 정밀한 유체시뮬레이션을 위해서 스플레쉬 후의 입자의 처리가 물리적으로 고려되어야 할 사안이다.

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]
  • J.J Buckels, R.D Hazlett, S Chen, K.G Eggert, W Grunau, W.E Soll, “Toward Improved Prediction of Reservoir Flow Performance - Simulating Oil and Water Flows at the Pore Scale”, Los Alamos Science, (1994), 22, p112-121.
  • 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]
  • S Chen, D Martinez, R Mei, “On boundary conditions in Lattice Boltzmann methods”, Phys Fluids, (1996), 8, p2527-2536. [https://doi.org/10.1063/1.869035]
  • D Chowdhury, D Stauffer, “Systematics of the models of immune response and autoimmune disease”, Journal of Stat. Phys, (1990), 59, p1019-1042. [https://doi.org/10.1007/BF01025860]
  • C.W Hirt, B.D Nichols, “Volume of fluid(VOF) method for the dynamics of free boundaries”, Journal of Computational Physics, (1981), 39, p201-225. [https://doi.org/10.1016/0021-9991(81)90145-5]
  • 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]
  • L Huabing, L Xiaoyan, F Haiping, Q Yuehong, “Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimentions”, Physical Review E, (2004), 70, p1539-3755.
  • R.-T Jung, “Numerical simulation on phase separation by using the Lattice-Boltzmann Method”, J. the Korean Society for Marine Environmental Engineering, (2009), 12(3), p197-201, printed in Korean.
  • 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]
  • I Kim, “Second order bounce back boundary condition for the Lattice Boltzmann Fluid Simulation”, KSME International Journal, (2000), 14, p84-92.
  • C Korner, M Thies, T Hofmann, N Thurey, 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]
  • 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]
  • B.T Nadiga, S Zaleski, “Investigations of a twophase fluid model”, Jornal of Mechanics, B/Fluids, (1996), 15(n-6), p885-896.
  • K Nagel, C.L Barrett, “Using Microsimulation Feedback for Trip Adaptation for Realistic Traffic in Dallas”, Int. J. Mod. Phys. C, (1997), 8, p505-525. [https://doi.org/10.1142/S0129183197000412]
  • 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]
  • M Resnics, “Turtles, Termites, and Traffic Jams - Explorations in Massively Parallel Microworlds”, A Bradford Book, (1995), MIT Press.
  • G.R Ryskin, L.G Leal, “Numerical solution of freeboundary problems in fluid mechanics. Part 1. The finite-difference technique”, J. Fluid Mech, (1984), 148, p1-17. [https://doi.org/10.1017/S0022112084002214]
  • X Shan, H Chen, “Lattice boltzmann model for simulating flows with multiple phased and components”, Physical Review E, (1993), 47(3), p1825-1819.
  • M Sussman, P Smereka, S Osher, “A level set approach for computating solutions to incompressible two-phase flow”, Journal of Computational physics, (1994), 114, p272-280. [https://doi.org/10.1006/jcph.1994.1155]
  • S Takagi, A Prosperetti, Y Matsumoto, “Drag coefficient of a gas bubble in an asymmetric shear flow”, Physics of Fluids, (1994), 6, p3186-3188. [https://doi.org/10.1063/1.868097]
  • N Thurey, “A single-phase free-surface Lattice Boltzmann Method”, Master thesis, INSTITUT FUR INFORMATIK (MATHEMATISCHE MASCHINEN UND DATENVERARBEITUNG), (2003).
  • N Thurey, C Wojtan, M Gross, G Turk, “A multiscale approach to mesh-based surface tension flows”, ACM SIGGRAPH 210 papers, (2010).
  • G Tryggvason, , “Direct numerical simulations of Gas-Liquid Multiphase flows”, Cambridge University Press, (2011). [https://doi.org/10.1016/j.fluiddyn.2005.08.006]
  • B Wen, H Li, C Zhang, H Fang, “Lattice-typedependent momentum-exchange method for moving boundaries”, Physical Review E, (2012), 85.
  • F Wolfram, W Tino, T Holger, S Hans-Peter, “Smoke Surfaces: An Interactive Flow Visualization Technique Inspired by Real-World Flow Experiments”, IEEE Transactions on Visualization and Computer Graphics (Proceedings Visualization 2008), (2008), 14(6).
  • W Xiaoming, Z Ye, F Zhe, L Wei, Y.-S Suzanne, K Arie, “Blowing in the Wind”, Eurographics/SIGGRAPH Sympposium on Computer Animation, The Eurographics Association, (2003), p75-85.
  • T Yabe, F Xiao, T Utsumi, “The Constrained Interpolation Profile Method for Multiphase Analysis”, JCP, (2001), 169, p556-593. [https://doi.org/10.1006/jcph.2000.6625]

Fig. 1

Fig. 1
Momentum distributions in two dimensional rectangular site

Fig. 2

Fig. 2
The single-time relaxation process(a) Particle momentum that will collide at the site P(b) Equilibulium momentum distribution when τ=1(c) Outgoing momentum distribution when τ=0.5

Fig. 3

Fig. 3
Problem setup for the density gradient

Fig. 4

Fig. 4
(a) The density as function of non-dimensional depth for a fluid with ρ=1.0 and τ=0.596. The arrow is the density profile for the steady state.(b) The density variation of a point at the bottom layer as the time is evolved.

Fig. 5

Fig. 5
Schematic view of typical Dam-breaking

Fig. 6

Fig. 6
Mass conservation check in case of 100×100

Fig. 7

Fig. 7
Comparison of the water front position with experiments and other simulation results

Fig. 8

Fig. 8
Free surface movement under the gravitational field

Table 1

Parameters values for the test simulation

Table 2

Parameters values for the grid dependency test