DEVELOPMENT AND VALIDATION OF A FLOW SOLVER BASED ON THE ENTROPIC LATTICE BOLTZMANN METHOD

J. Kim; Y. Jo

doi:10.6112/kscfe.2024.29.4.233

Preview

Original Article

Journal of Computational Fluids Engineering. 31 December 2024. 233-246
https://doi.org/10.6112/kscfe.2024.29.4.233

DEVELOPMENT AND VALIDATION OF A FLOW SOLVER BASED ON THE ENTROPIC LATTICE BOLTZMANN METHOD

엔트로피 격자 볼츠만 기법을 기반으로 한 유동해석 코드 개발 및 검증

J. Kim¹

Y. Jo²^†

김 준재¹

조 영민²^†

¹Dept. of Aeronautical System Engineering, Hanseo University

²Dept. of Unmanned Aircraft Systems, Hanseo University

¹한서대학교 항공시스템공학과

²한서대학교 무인항공기학과

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted noncommercial use, distribution, and reproduction in anymedium, provided the original work is properly cited.

ABSTRACT

The lattice Boltzmann method(LBM) has been widely used to model fluid dynamics by reproducing the Navier-Stokes equations, and successfully handling complex phenomena such as multi-phase, multi-component, and microscale flows. Despite its versatility, the standard LBM faces significant limitations when applied to high Mach number and low-viscosity flows, where numerical inaccuracies become prominent. To overcome these challenges, Karlin et al.[Europhys. Lett., 47, pp. 182–188(1999)] introduced the entropic lattice Boltzmann method(ELBM), a significant advancement that employs Boltzmann’s H-theorem and the maximum entropy principle(MEP) to ensure compliance with the second law of thermodynamics, which is often violated by the standard LBM. In this study, we present the development and validation of a flow solver based on the ELBM, implemented with GPU parallel computing. The solver’s performance is demonstrated through simulations of flow over a 2D circular cylinder and the NACA0012 airfoil, validating its accuracy.

Keywords

Computational Fluid Dynamics, CFD

Lattice Boltzmann Method, LBM

Entropic Lattice Boltzmann Method, ELBM

D2Q21 lattice

2D circular cylinder

2D NACA0012 airfoil

Grad’s boundary condition

키워드

전산유체역학

격자볼츠만기법

엔트로피격자볼츠만기법

D2Q21격자

2차원 원형 실린더

2차원 NACA0012

Grad경계조건

MAIN

1. 서 론
2. 이론적 배경
2.1 격자 볼츠만 기법(LBM)
2.2 엔트로피 격자 볼츠만 기법
2.3 Grad’s moment approximation을 이용한 No-slip 경계 조건
2.4 공력 계산
3. 코드 개발 및 검증
3.1 2차원 원형 실린더에 대한 유동해석
3.2 2D NACA0012 Airfoil에 대한 유동해석
4. 결 론

1. 서 론

격자 볼츠만 기법(Lattice Boltzmann Method, LBM)은 Navier-Stokes 방정식을 풀기 위한 대안으로 제시된 방법으로, 특히 다상 유동(multi-phase flow), 다성분 유동(multi-component flow), 다공성 매질(porous media) 등 복잡한 유동 해석 문제에서 성공적으로 적용되고 있다[1, 2]. LBM은 단순한 지배 방정식과 해석 알고리즘을 바탕으로 하며, 병렬 계산에 적합한 특성 덕분에 최근에는 GPU를 활용한 유동 해석 소프트웨어로 활발히 개발되고 있다.

LBM 알고리즘은 크게 두 가지 단계로 구성된다. 첫 번째는 입자 간 충돌 과정(collision process)이며, 두 번째는 충돌 후 입자들이 사전에 정의된 경로를 따라 인접 지점으로 이동하는 과정(propagation process)이다. 충돌 과정은 전적으로 지역화(localized)되어 있어 각 격자 노드에서 독립적으로 계산될 수 있으며, 이동 과정에서는 인접한 값만 참조하기 때문에 GPU를 활용한 병렬 계산에 최적화되어 있다[3].

이러한 LBM의 강력한 장점에도 불구하고, 이 기법의 적용은 여전히 높은 점성과 낮은 마하 수에 대한 유동해석으로 제한되고 있다. 이는 운동 점성(kinematic viscosity)의 함수인 relaxation factor가 높아짐에 따라 Fig. 1과 같이 시뮬레이션의 안정성이 급격하게 저하되기 때문이다[4]. Relaxation factor란 입자들이 평형 상태로 돌아가고자 하는 경향성을 나타내며, relaxation factor가 0.5일 때 입자는 충돌 과정 후 즉시 평형 상태가 되고 1.0에 가까워질수록 천천히 평형 상태에 도달한다. 또한, 한 지점에서 특정 속도를 가질 확률을 나타내는 입자 분포 함수는 Navier-Stokes 방정식에서 사용되는 모멘트보다 더 고차원의 모멘트(high-order moments)들에 대한 정보를 포함하고 있으며, LBM 알고리즘 내에서 이러한 고차 모멘트의 정보가 적절히 계산되지 못하면(under-resolved) 비물리적 현상이 발생할 수 있다[5]. 가장 단순한 Bhatnagar- Gross-Krook(BGK) 충돌 모델을 사용하는 LBM의 경우, 항상 양수인 입자 분포 함수가 특정 상황에서 음의 값을 가지거나, 갈릴레오 불변성(Galilean invariance)을 충족하지 못하거나, 입자들의 충돌 과정 전후로 계의 엔트로피가 감소하는 현상이 발생하여 열역학 제2 법칙을 위배하는 수치적인 문제가 발생하기도 한다[6]. 이러한 문제로 인하여 기존 LBM은 relaxation factor가 0.9 이하에서만 안정적으로 유동 해석이 가능했다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F1.jpg

Fig. 1.

Stability regions depending on relaxation factor(𝛽)

이러한 LBM의 근본적인 문제를 해결하기 위해 많은 연구가 진행되어 왔다. 특히, 각각의 모멘트를 다른 속도로 완화시키는 다중 완화 시간(Multi-Relaxation Time, MRT) 모델이 대표적인 예시 중 하나로, LBM 시뮬레이션의 정확성과 안정성을 획기적으로 개선하는 데 기여했다[7]. 그러나 다수의 relaxation factor들을 free parameter로서 설정해야 하기 때문에, 최적의 relaxation factor들을 찾기 위해서는 많은 경험적 테스트가 동반되어야 했다. 한편, BGK-LBM의 엔트로피 감소 문제를 해결하기 위해 Karlin et al.[8]은 볼츠만의 H-정리를 활용하여 계의 엔트로피가 감소하지 않도록 하고 매 time-step간의 엔트로피의 변화를 최소화함으로써 시뮬레이션의 안정성을 향상시킨 ELBM(Entropic LBM)을 제안하였다. 그리고 Frapolli et al.[9]은 고 마하수 유동에 대한 시뮬레이션의 안정성과 정확성을 향상시키고 갈릴레오 불변성을 충족하기 위해 d-차원의 유동에 대해 최소 $D d Q 7^{d}$ 이상의 격자 모델을 사용할 것을 제안하였다. 여기서 $Q$ 는 입자들의 이산화된 속도의 수를 의미한다. 최근에는 Karlin과 그의 연구팀이 극초음속(Ma > 5) 유동 해석을 위해 Particle on Demand(PonD) 모델을 적용한 LBM 모델을 제안하기도 하였다[10, 11].

가장 진보된 LBM 모델 중 하나로 평가받는 ELBM은 상대적으로 높은 레이놀즈 수의 유동에 대한 해석이 가능하며, 여러 연구에서 그 유효성이 검증되어 왔다[12, 13]. 본 연구에서는 ELBM을 활용한 in-house 코드를 소개하고, 높은 relaxation factor(𝛽≥0.90)에서의 시뮬레이션 정확도와 안정성에 대한 검증 결과를 분석하며, 유동 수치해석 기법으로서 LBM의 가능성을 고찰하였다. 본 논문의 Section 2에서는 ELBM의 이론적 배경을 소개하고, Section 3에서는 in-house 코드를 이용한 2차원 원형 실린더 및 2차원 NACA0012 에어포일에 대한 해석 결과를 다른 수치해석 결과와 비교하였다.

2. 이론적 배경

2.1 격자 볼츠만 기법(LBM)

LBM은 평형 상태일 때 특정 속도를 가지는 입자들의 분포를 나타내는 맥스웰-볼츠만 분포 함수와 이러한 입자 분포의 변화를 계산하는 볼츠만 운동 방정식을 활용한다.

(1)

\frac{d}{d t} f (x_{β}, ξ_{β}, t) = \frac{\partial f}{\partial t} + ξ_{β} \frac{\partial f}{\partial x_{β}} + \frac{\vec{F_{β}}}{ρ} \frac{\partial f}{\partial ξ_{β}} = Ω (f)

위 방정식에서 시간, 위치, 입자 속도(phase space)에 따른 입자 분포 함수 $f$ 의 변화는 입자 간 충돌 𝛺에 의해 설명된다. 이를 이산화하면 아래와 같이 표현할 수 있다.

(2)

f_{i} (x + c_{i} Δ t, t + Δ t) - f_{i} (x, t) = Ω_{i} (x, t)

여기서 $\vec{c_{i}}$ 는 이산화된 입자들의 속도를 의미하며, $i$ 는 입자들의 이동 방향 번호를 의미한다(D2Q9 격자 모델의 경우 $i ∊ 0, 1, 2, . . ., 8$ ). 본 연구에서는 입자 간 충돌을 일대일 충돌로 가정하며(dilute gas), 모든 미시적 물리현상을 무시하는 single-relaxation factor BGK 충돌 모델을 사용하였으며, 이는 식 (3)과 같이 표현될 수 있다.

(3)

Ω^{B G K} = 2 β (f^{e q} - f)

$f^{e q}$ 는 평형 상태에서의 입자 분포 함수이며, 𝛽는 relaxation factor이다. 이 relaxation factor는 운동 점성 𝜈와 온도 $T$ 의 함수로, 식 (4)과 같이 정의될 수 있다.

(4)

β = \frac{T}{2 ν + T}

위의 방정식들을 적용하여, 충돌 과정 후 입자 분포 함수를 식 (5)와 같이 표현할 수 있다. 여기에서 𝛼는 ELBM에서 사용되는 over-relaxation factor이며, 국소적인 점성(local viscosity)을 제어하는 핵심적인 요소이다. 입자 분포 함수가 평형 상태에 근접하다면 over-relaxation factor는 2가 되며, 이 경우 일반적인 BGK 모델과 동일하다.

(5)

f^{*} = f + α β (f^{e q} - f), α = 2 (B G K)

충돌 과정 후, 입자들은 식 (6)와 같이 지정된 격자 속도(lattice velocity)를 따라 이동하게 되며 이를 전파 과정이라 한다.

(6)

f (\vec{x} + {\vec{c}}_{i} Δ t, t + Δ t) = f^{*} (\vec{x}, t)

충돌과 전파 과정을 수행한 후 입자 분포 함수를 이용하여 유체의 밀도(𝜌), 운동량( $ρ \vec{u}$ ), 총 에너지( $E$ )와 같은 거시적인 물리량을 식 (7)과 같이 구할 수 있다. 추후 고차 모멘트를 고려하기 위해, 이를 일반적인 형태로 식 (8)과 같이 표현할 수 있다.

(7)

ρ = \sum_{i}^{Q} f_{i} ρ u_{α} = \sum_{i}^{Q} f_{i} c_{i α} 2 ρ E = \sum_{i}^{Q} f_{i} c_{i}^{2}

(8)

\sum_{i}^{Q} f c_{i x}^{p} c_{i y}^{q} c_{i z}^{r} = M_{p q r}

2.2 엔트로피 격자 볼츠만 기법

2.2.1 지수함수 형태의 평형 분포 함수(Exponential Equilibrium Distribution Function)

일반적인 LBM에서는 평형 분포 함수를 구하기 위해 이산화된 맥스웰-볼츠만 평형 분포 함수의 이차 다항식 형태의 방정식을 사용한다[4]. 이는 가장 간단하고 효율적인 방법이지만, 삼차 이상의 항에 대해서는 고려하지 않으므로, 마하수가 높거나 점성이 낮을 때 무시할 수 없는 오류가 발생한다( $O (u^{3})$ ). 이러한 경우 평형 분포 함수가 양의 값을 잃는 문제가 발생할 수 있다[14]. 이를 해결하기 위해 Karlin, Succi, 그리고 그 연구팀은 최대 엔트로피 원리(Maximum Entropy Principle, MEP) 및 볼츠만의 H-정리를 기반으로 하여 매 time-step마다 엔트로피의 변화가 최소화되는 방향으로 지수함수 형태의 평형 분포 함수를 계산하는 방법을 제안하였다[15]. 이를 통해 어떤 조건에서도 평형 분포 함수가 양의 값을 가질 수 있으며, 따라서 보다 안정적이고 정확한 시뮬레이션이 가능하게 되었다. 이산화된 H-정리를 이용한 엔트로피 계산은 다음과 같다.

(9)

H = \sum_{i}^{Q} f_{i} [\ln (\frac{f_{i}}{a})]

여기서 $a$ 는 일반적으로 분포 함수의 가중치(weight)를 사용한다. 이와 같은 H-정리와 라그랑주 승수법(Lagrange multiplier), 그리고 뉴턴-랩슨(Newton-Raphson) 반복법을 사용하여 다음과 같은 방법으로 평형 분포 함수를 계산할 수 있다[16]. 한 위치에서 입자들의 밀도, 운동량, 그리고 에너지와 같은 물리량은 입자들이 평형 상태에 도달하더라도 보존 법칙에 따라 그 값이 변하지 않아야 한다. 따라서, 이러한 보존 법칙은 식 (10)과 같이 뉴턴-랩슨 과정에서의 제약 조건(constraint)으로 활용될 수 있다.

(10)

G_{ρ} = \sum_{i} f_{i}^{e q} - ρ = 0, G_{ρ u_{α}} = \sum_{i} f_{i}^{e q} c_{i α} - ρ u_{α} = 0, G_{ρ E} = \sum_{i} f_{i}^{e q} c_{i}^{2} - 2 ρ E = 0

최대 엔트로피 원리에 따라 평형 상태에서는 엔트로피의 변화가 없어야 하므로 $\nabla_{f_{i}} H = \frac{\partial H}{\partial f_{i}} = 0$ 이 되도록 설정하고 라그랑주 승수를 이용하여 식 (11)과 같이 정리할 수 있다.

(11)

- \frac{\partial H}{\partial f_{i}} = λ_{ρ} \frac{\partial G_{ρ}}{\partial f_{i}} + λ_{ρ u_{α}} \frac{\partial G_{ρ u_{α}}}{\partial f_{i}} + λ_{ρ E} \frac{\partial G_{ρ E}}{\partial f_{i}}

다음으로, 볼츠만 H-정리를 $f_{i}$ 에 대해 편미분을 취하면 $\nabla_{f_{i}} H = \ln f_{i} + 1$ 이 되므로 평형 분포 함수는 최종적으로 식 (12)과 같이 표현될 수 있다.

(12)

f_{i}^{e q} = a \exp [- (1 + \sum_{p, q, r} λ_{M_{p q r t}^{e q}} c_{i x}^{p} c_{i y}^{q} c_{i z}^{r})]

볼츠만 H-정리에 기반한 지수함수 형태의 평형 분포 함수를 사용하면 높은 마하수 및 낮은 점성 유동에 대해서도 양의 값이 보장될 수 있다. 그러나, 평형 분포 함수를 구하는 방정식이 닫힌 형태(closed-form)로 존재하지 않기 때문에, 매 time-step마다 뉴턴-랩슨 과정을 반복해야 하는 번거로움이 있다. 또한, 적용된 격자 모델에 따라 사용해야 하는 제약 조건의 수가 더 많아질 수 있다. 본 연구에서는 Fig. 2과 같이 3-level $D 2 Q 7^{2}$ 격자 모델에서 간소화된 $D 2 Q 2^{1}$ 격자 모델을 사용하였으며, 이에 따라 식 (13)과 같이 3차 모멘트를 포함하여 총 8개의 제약 조건을 적용하였다[17]. 여기서 0차 모멘트는 밀도, 1차 모멘트는 운동량, 2차 모멘트는 운동량 유속(momentum flux), 그리고 3차 모멘트는 열 유속(heat flux)에 대한 값을 나타낸다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F2.jpg

Fig. 2.

D2Q21 3-level lattice velocity model

(13)

M_{00} = \sum_{i}^{Q} f_{i} = ρ M_{10} = \sum_{i}^{Q} f_{i} c_{i x} = ρ u_{x} M_{01} = \sum_{i}^{Q} f_{i} c_{i y} = ρ u_{y} M_{20} = \sum_{i}^{Q} f_{i} c_{i x}^{2} = ρ (u_{x}^{2} + T) M_{11} = \sum_{i}^{Q} f_{i} c_{i x} c_{i y} = ρ u_{x} u_{y} M_{02} = \sum_{i}^{Q} f_{i} c_{i y}^{2} = ρ (u_{y}^{2} + T) M_{30} = \sum_{i}^{Q} f_{i} c_{i x} c_{i}^{2} = 2 ρ u_{x} (\frac{1}{2} (u_{x}^{2} + u_{y}^{2} + 2 T) + T) = 2 ρ u_{x} (E + T) M_{03} = \sum_{i}^{Q} f_{i} c_{i y} c_{i}^{2} = 2 ρ u_{y} (\frac{1}{2} (u_{x}^{2} + u_{y}^{2} + 2 T) + T) = 2 ρ u_{y} (E + T)

2.2.2 Over-relaxation factor를 이용한 충돌 과정

ELBM의 핵심은 볼츠만 H-정리를 이용한 평형 분포 함수의 도출뿐만 아니라, over-relaxation factor를 정의하는 데 있다. 입자 충돌 과정에서 엔트로피 변화가 최소화되는 방향으로 over-relaxation factor, 𝛼를 설정하는 것이 중요한 요소이다. Over-relaxation factor는 각 격자 노드에서 개별적으로 정의되며, 이를 통해 충돌 과정에서의 안정성과 정확성을 높일 수 있다. 국소적인 𝛼를 계산하기 위해서 마찬가지로 H-theorem를 이용한다[8].

(14)

H (f + α (f^{e q} - f)) = H (f)

위의 엔트로피 조건을 만족하는 양의 해(positive root)를 𝛼로 설정한다. 이를 통해 점성이 0에 근접해도 (𝛽→1) 충돌 과정 후 엔트로피가 감소하지 않도록 할 수 있다[18]. 그러나, 엔트로피 조건을 만족하는 𝛼를 구하기 위해서는 또다른 뉴턴-랩슨 반복 과정을 사용해야 하기 때문에 연산비용이 증가된다. 따라서, 분포 함수가 평형 상태에 근접해 있다고 가정할 수 있다면, 식 (15)과 같은 명시적(explicit) 형태로 over-relaxation factor를 대략적으로 계산하여 추가적인 뉴턴-랩슨 반복과정을 줄일 수 있다[19].

(15)

α = 2 + \frac{2 Δ t}{3 β} (\frac{S_{i j} S_{i k} S_{j k}}{S_{l m} S_{l m}})

여기서 $S_{α β} = \frac{1}{2} (\partial_{α} u_{β} + \partial_{β} u_{α})$ 이다.

2.3 Grad’s moment approximation을 이용한 No-slip 경계 조건

일반적으로 LBM에서는 no-slip 경계 조건을 적용하기 위해 half-way bounce-back(HBB) 기법이 사용된다. 이 방법은 물체(obstacle)의 경계면(boundary)으로 향하는 입자들이 동일한 경로를 통해 반사된다는 개념에 기반하여, 가장 간단하게 질량 보존을 유지하면서 no-slip 조건을 적용할 수 있는 방법이다[4]. 그러나 이 방법은 경계면을 격자를 따라 계단식으로 표현하게 되어 경계면이 왜곡된다는 단점이 있다. 이를 보완하기 위해, 경계면에 가장 인접한 경계 노드(boundary node)에서 Grad의 모멘트 근사(Grad’s moments approximation)를 사용하여 비평형 함수 $(f^{neq} = f - f^{eq})$ 를 계산하는 Grad 경계 기법을 적용하였다[20]. 이를 위해 우선적으로 경계 노드에서 필요한 밀도와 속도를 계산해야 한다. 질량 보존을 위해 밀도는 식 (16)과 같이 HBB를 이용해 계산할 수 있으며, 속도의 경우 경계면에서 속도를 0으로 설정하고(no-slip 경계 조건), 입자가 이동하는 방향과 거리에 따라 보간법(interpolation)을 통해 경계 노드에서 필요한 속도를 식 (17)과 같이 계산할 수 있다.

(16)

ρ^{target} = \sum_{i \in \bar{D}} f_{i}^{H B B} + \sum_{i \notin \bar{D}} f_{i}

(17)

{\vec{u}}^{target} = \frac{1}{n_{\bar{D}}} \sum_{i \in \bar{D}} \frac{q_{i} {\vec{u}}_{f, i}}{1 + q_{i}}

여기서 $\bar{D}$ 는 missing population( $f_{i} (i \in \bar{D}) = f_{i}^{m i s s}$ ), 즉, 유체 영역(fluid domian) 밖에서 내부로 들어오는 입자를 의미한다. $n_{\bar{D}}$ 는 각 노드에서 missing population의 수이다. ${\vec{u}}_{f, i}$ 는 경계 노드에 인접한 유체 노드(fluid node)에서의 속도를 의미하며, $q_{i}$ 는 유체 노드로부터 물체의 경계면까지 정규화(normalized) 된 거리를 의미한다.

(18)

q_{i} = \frac{||\vec{x_{b}} - \vec{x_{w, i}}||}{||\vec{c_{i}}||}

계산된 목표 밀도와 속도를 이용하여 경계 노드에서의 비평형 입자 분포 함수를 식 (19)과 같이 계산할 수 있다.

(19)

f_{i}^{n e q, Grad} = W_{i} (\frac{P_{α β}^{(1)} (c_{i α} c_{i β} - T_{0} δ_{α β})}{2 T_{0}^{2}})

여기서 $W_{i}$ 는 각 입자 분포 함수가 갖는 가중치(weight)이며, $δ_{α β}$ 는 Kronecker delta이다. 그리고 $T_{0}$ 는 기준 온도(reference temperature)이다. LBM에서 기준 온도는 사용하는 격자 모델에 따라 달라진다. 예를 들어, D2Q9의 경우 기준 온도는 $\sqrt{1 / 3}$ , D2Q49의 경우 기준 온도는 0.697954이다[6]. 본 연구에서 사용한 D2Q21 격자 모델에 대해서는 기준 온도를 0.7로 설정하였다[16]. 비평형 압력 텐서 $P_{α β}^{(1)}$ 는 식 (20)을 통해 구할 수 있다.

(20)

P_{α β}^{(1)} = - \frac{1}{2 β} ρ T_{0} S_{α β}

마지막으로 목표 밀도와 속도를 이용하여 경계 노드에서의 평형 분포 함수를 계산하고 이를 Grad의 모멘트 근사를 이용하여 계산한 비평형 분포 함수와 더하면 식 (21)과 같이 경계 노드에서의 missing population을 구할 수 있다.

(21)

f_{i}^{m i s s} = f_{i}^{e q} + f_{i}^{neq, Grad} (i \in \bar{D})

2.4 공력 계산

LBM에서 물체 벽면에 작용하는 힘을 평가하는 방법으로는 두 가지가 널리 사용된다: 응력 적분(stress integration, SI)과 운동량 교환(momentum exchange analysis, MEA) 방법이다. 운동량 교환 방법은 특히 HBB 경계 조건을 사용하는 LBM 모델에 적합하다. 이는 HBB을 적용하는 과정에서 간단하게 운동량 변화를 계산할 수 있기 때문이다. 운동량 교환 방법의 유용성에도 불구하고, 시뮬레이션의 격자 해상도가 충분하지 않을 경우 왜곡되는 결과를 초래할 수 있다. 따라서, 낮은 해상도의 시뮬레이션에서도 그 정확성을 높이기 위해 Kaushal et al.[22]은 이산 레이놀즈 수송 정리(Discrete Reynolds Transport Theorem, DRTT)를 적용한 새로운 운동량 변화 계산 방법을 제안했다. LBM 프레임워크에서 DRTT는 일반적인 유한체적법(Finite Volume Method, FVM)과 같이 control volume 안의 운동량 변화를 계산하고, control surface를 통해 전달되는 운동량 교환을 계산한다. 이러한 운동량의 변화를 합산하여 공력 $F$ 를 계산할 수 있으며, 그 계산식은 식 (22)과 같다:

(22)

F_{α} (t) = J_{α}^{C V} (t) - J_{α}^{C V} (t + Δ t) + \sum_{x \in C S_{out}} \sum_{i \in D_{in}} f_{i} (x, t) c_{i α} - \sum_{x \in C S_{in} i \in D_{out}} f_{i} (x, t) c_{i α}

여기서 $J_{α}^{C V}$ 는 control volume 내의 운동량을 나타낸다. $x ∊ C S_{o u t}$ 은 control surface의 바깥쪽에 있는 노드를 나타내며, $i ∊ D_{in}$ 은 control volume 밖에서 내부로 향하는 입자 분포를 의미한다. 마찬가지로, $x ∊ C S_{in}$ 은 control surface 안쪽에 있는 노드를 나타내며, $i \in D_{o u t}$ 은 control volume 내에서 외부로 향하는 입자 분포를 의미한다.

응력 적분과 운동량 교환 방법의 경우 물체의 실제 표면이 아닌 가장 인접한 경계 노드에서의 운동량 변화를 계산하기 때문에 물체 경계면에서의 운동량 변화와 다소 차이가 있을 수 있다. 따라서 격자의 해상도가 충분하지 않을 때 경계 노드와 물체 경계면 간의 거리가 증가함에 따라 운동량의 변화량을 더욱 심하게 왜곡할 수 있다. 반면에 DRTT를 사용하면 격자의 해상도, 또는 실제 물체의 형상과 상관없이 control volume 내에서 물체에 의해 발생하는 운동량 변화를 계산할 수 있기 때문에 격자의 해상도가 낮더라도 더 정확하게 운동량 변화를 계산할 수 있다[22].

3. 코드 개발 및 검증

본 연구에서는 ELBM과 Python 기반의 CuPy 라이브러리를 이용하여 GPU 컴퓨팅을 활용한 in-house 코드를 개발하였다. LBM의 특성상 각 격자 노드가 독립성을 가지므로 GPU를 이용한 병렬 계산에 적합하다[4]. 따라서 ELBM 기반의 유동해석 알고리즘을 Fig. 3과 같이 구성하였다. 일반적인 LBM 알고리즘에서는 초기화 다음 충돌 과정으로 시작하지만, 본 연구에서는 전파 과정를 먼저 수행하였다. 이로써 충돌 후 입자 분포 함수(post-collision)에 대한 변수를 별도로 생성할 필요가 없으므로 GPU 메모리 효율성을 높일 수 있다[23]. 또한, 평형 분포 함수를 계산하는 과정 중 라그랑주 승수에 대한 값을 구하기 위하여 뉴턴-랩슨 반복법을 사용하였으며, 식 (23)과 같이 Jacobian 행렬을 구성하고 그 역행렬을 이용함으로써 라그랑주 승수의 값을 결정할 수 있다[16].

(23)

J = (\begin{matrix} \sum_{i} f_{i}^{e q} & \sum_{i} f_{i}^{e q} c_{i α} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} & \sum_{i} f_{i}^{e q} c_{i α} c_{i}^{2} \\ \sum_{i} f_{i}^{e q} c_{i α} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i γ} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i}^{2} \\ \sum_{i} f_{i}^{e q} c_{i α} c_{i β} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i γ} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i γ} c_{i κ} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i γ} c_{i}^{2} \\ \sum_{i} f_{i}^{e q} c_{i α} c_{i}^{2} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i}^{2} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i γ} c_{i}^{2} & \sum_{i} f_{i}^{e q} c_{i α} c_{i β} c_{i}^{2} c_{i}^{2} \end{matrix})

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F3.jpg

Fig. 3.

Simulation algorithm of in-house code based on the entropic lattice Boltzmann method

D2Q21 격자 모델과 8개의 구속 조건을 사용하는 경우 평형 분포 함수를 계산하는 데 있어 각 격자 노드에서 21×8×8형태의 텐서 변수를 생성하여야 하므로 이는 상당한 GPU 메모리를 차지하는 요소이다. 따라서 본 연구에서는 GPU가 가진 thread를 최대한 활용하되, 최대 memory-pool의 크기를 줄이기 위하여 시뮬레이션 도메인을 여러 서브 도메인으로 나누어 순차적으로 계산을 진행하였다.

3.1 2차원 원형 실린더에 대한 유동해석

엔트로피 격자 볼츠만 기법 기반의 유동해석 코드의 정확성과 성능을 검증하기 위해 우선 2차원 원형 실린더에 대한 시뮬레이션을 수행하였다. 원형 실린더에 대한 해석은 수치해석 기법의 검증에 널리 사용되어 왔으며, 이에 대한 풍부한 실험적 자료가 존재한다. 그러나 2차원 해석은 3차원에 대한 실험결과와 차이가 발생하므로 2차원 유동에 대한 다른 수치해석 결과와 비교를 진행하였다. 비교에 사용된 Henderson의 시뮬레이션에서는 unstructured spectral element method를 사용하였다[24, [25]. 유동해석은 레이놀즈수 10부터 100까지는 10 단위로, 100부터 1,000까지는 100 단위로 해석을 수행하여 비교 자료와 동일한 조건을 유지하였다. 시뮬레이션 도메인은 Fig. 4와 같이 구성하였으며, 노드 수는 코드 길이 기준으로 모든 레이놀즈수에 대해 C = 60으로 설정하였다. Inflow boundary 조건은 밀도, 속도, 그리고 온도를 할당하여 분포 함수가 식 (24)과 같이 평형 분포 함수와 일치하도록 설정하였다.

(24)

f_{i} (x_{in}, t) = f_{i}^{eq} (x_{in}, t) (ρ^{set}, u_{α}^{set}, T^{set})

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F4.jpg

Fig. 4.

Simulation domain for a flow past 2D circular cylinder

Outflow boundary는 missing population에 대해서만 이전 time-step의 값과 같도록 식 (25)과 같이 설정하였다.

(25)

f_{i}^{miss} (x_{out}, t) = f_{i}^{miss} (x_{out}, t - Δ t)

Open boundary의 경우, missing population에 대해서만 $f_{i}^{miss} (x_{open}, t) = f_{i}^{e q} (x_{open}, t)$ 로 설정하였는데, 이는 가장 간단한 방법이지만 비평형 분포 함수에 대한 값이 소실되므로 시뮬레이션의 정확도에 영향을 줄 수 있다. 따라서 이러한 방식으로 경계 조건을 설정하기 위해서는 충분한 크기의 도메인이 보장되어야 한다. 위와 같은 조건을 사용하여 2차원 원형 실런더에 대한 항력 계수를 계산하고, 비교 자료의 결과와 비교하였다. 비교 자료에서는 점성에 의한 항력과 압력항력을 별개로 계산하였으나, DRTT는 이를 구분하여 계산할 수 없으므로 전체 항력 값만 비교하여 Fig. 5(a)에 제시하였다. 결과적으로 본 연구에서 개발한 in-house 코드를 이용하여 해석을 수행한 모든 레이놀즈수 구간에서 비교 자료의 수치해석 결과와 잘 일치함을 볼 수 있다. 또한, 각 레이놀즈수에 대한 시뮬레이션에서 적용된 relaxation factor를 Fig. 5(b)에 제시하였으며, 이는 ELBM을 사용하였을 때 relaxation factor가 0.9 이상에서도 상대적으로 시뮬레이션이 안정적임을 보여준다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F5.jpg

Fig. 5.

(a) Drag coefficients and (b) relaxation factor for 2D circular cylinder at various Reynolds number

3.2 2D NACA0012 Airfoil에 대한 유동해석

다음으로 레이놀즈수 1,000에서 2차원 NACA0012 에어포일에 대한 유동해석을 진행하였으며, 받음각 0도에서 29도까지 1도 간격으로 해석을 진행하였다. 시뮬레이션 도메인을 Fig. 6와 같이 구성하였으며, 경계 조건은 2차원 원형 실린더 시뮬레이션의 조건과 동일하게 설정하였다. 추가적으로 받음각이 증가함에 따라 난류의 영향이 중요해지므로 이를 잘 해석할 수 있도록 노드 수를 C = 200으로 설정하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F6.jpg

Fig. 6.

Simulation domain for a flow past 2D NACA0012 airfoil

해석 결과, 받음각 8도부터 Kármán vortex street가 관찰되기 시작하였으며, 이는 비교 자료인 Kurtulus의 시뮬레이션 결과와 일치하였다[26]. 추가적으로 비교 자료의 시뮬레이션에서는 받음각 22도부터 periodic doubling bifurcation으로 인하여 vortex가 두 갈래로 갈라지는 현상이 관찰되었으나, 본 연구에서 개발한 코드에서는 Fig. 7와 같이 받음각 24도부터 관찰되었다. 이는 받음각 22도에서 에어포일의 상부에서 발생하는 주요 와류의 강도 차이에 기인한 것으로 판단된다. 받음각 28도 이후로는 chaotic wake가 발생함을 확인할 수 있었다. 추가적인 비교를 위해 DRTT를 이용하여 양력과 항력 계수를 계산하였으며, 이 결과를 유한체적법과 유한요소법(Finite Element Method, FEM) 기반의 수치해석 결과와 비교하여 Fig. 8에 제시하였다[26, 27, 28]. Fig. 8에서 보이는 바와 같이 유동이 비교적 안정적인 받음각 19도 이하에서는 다른 수치해석 결과들과 매우 유사한 값이 도출되었으며, 유동박리가 심한 받음각 영역에서도 전반적인 경향성이 유사함을 확인할 수 있었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F7.jpg

Fig. 7.

Velocity field for NACA0012 airfoil at AoA = 24°

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290419/images/jkscfe_2024_294_233_F8.jpg

Fig. 8.

(a) Lift and (b) drag coefficients for NACA0012 airfoil at various angle of attack

4. 결 론

ELBM의 핵심은 평형 분포 함수를 정의하는 과정과 충돌 과정에서 발생하는 엔트로피의 변화를 최소화하여 시뮬레이션의 안정성과 정확도를 향상시키는 데 있다. 이로써 LBM이 갖는 해석 능력이 낮은 relaxation factor에 국한되는 한계를 극복할 수 있었다. 따라서 본 연구에서는 ELBM 기반의 코드를 개발함으로써 ELBM이 실용적인 유동해석 기법으로 사용될 가능성을 확인하고자 하였다. 그리고 검증을 위하여 2차원 원형 실린더와 2차원 NACA0012 에어포일에 대한 해석을 수행하고 그 결과를 다른 수치해석 결과와 비교하였다. 결론적으로 다른 수치해석 기법과 매우 유사한 공력 값이 도출되었으며, 상대적으로 높은 relaxation factor 영역대에서도 기존 LBM 대비 그 안정성을 확인할 수 있었다. 따라서 ELBM이 신뢰할 수 있는 해석 기법임을 보여줄 수 있었다.

ELBM은 Smargorinsky subgrid scale model을 사용하는 Large Eddy Simulation(LES)와 유사하게 최대 엔트로피 원리(MEP)와 볼츠만 H-정리를 이용하여 over-relaxation factor를 결정하고 이를 통해 국소 점성을 조절함으로써, 시뮬레이션 안정성을 높일 수 있었다[18]. 그러나 ELBM은 여전히 Direct Numerical Simulation(DNS)에 가까운 시뮬레이션 기법으로서 기존 상용 유동해석 소프트웨어와 비교하였을 때 상당한 연산비용을 요구한다. 따라서 ELBM이 보다 실용적이고 범용적인 유동해석 기법으로 사용되기 위해서는 여전히 시뮬레이션 안정성과 효율성 개선이 필요하다. 이에 후속 연구에서는 Karlin et al.[29]이 제안한 Entropic Multi-Relaxation Time(EMRT) 기법을 적용한 ELBM을 테스트하고자 한다. EMRT에서는 입자 분포 함수를 운동항(kinetic part), 전단항(shear part), 그리고 고차 모멘트항(high-order moments)으로 분할하여 계산하고 기존의 MRT 모델과 유사하게 각 파트 별로 개별적으로 정의되는 relaxation factor를 적용한다. 이 과정이 또한 많은 연산비용을 요구하므로 이 과정을 간소화할 수 있는 추가적인 연구가 필요할 것으로 보인다. 또한, 압축성 유동에 대한 열 해석을 고려하기 위하여 유체의 비열(specific heat)을 제어할 수 있도록 에너지에 대한 입자 분포 함수를 별개로 할당하는 double distribution function(DDF) model[30]과 유체의 Prandtl number를 제어할 수 있도록 quasi-equilibrium collision model을 적용하여 테스트하고자 한다[31].

Acknowledgements

이 논문은 2024년도 정부(방위사업청)의 재원으로 국방기술진흥연구소의 지원을 받아 수행된 연구임(KRIT-CT-23-01, 국방수직이착륙기특화연구센터).

References

2020, Sahu, A. and Bhowmick, S., "Applications of lattice Boltzmann method in multi-component and multi-phase flow: A review," AIP Conf. Proc., Vol.2273, p.030007.

10.1063/5.0024322

2012, Nor Azwadi, C.S. and Afiq Witrib, M.Y., "Simulation of multicomponent multiphase flow using lattice Boltzmann method," AIP Conf. Proc., Vol.1440, pp.1025-1039.

10.1063/1.4704318

2023, Plekhanov, M., Ivashchenko, V., Karpenko, A. and Mullyadzhanov, R., "Numerical simulation of a turbulent pipe flow: FluidX3D LBM validation," E3S Web Conf., Vol.459, p.03010.

10.1051/e3sconf/202345903010

2017, Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G. and Viggen, E.M., "The Lattice Boltzmann Method: Principles and Practice," Springer International Publishing.

2001, Dellar, P.J., "Bulk and shear viscosities in lattice Boltzmann equations," Phys. Rev. E., Vol.64, p.031203.

10.1103/PhysRevE.64.03120311580323

2017, Frapolli, N., "Entropic lattice Boltzmann models for thermal and compressible flows," Ph.D. thesis, ETH Zurich.

1992, d'Humières, D., "Generalized lattice-Boltzmann equations," in Rarefied Gas Dynamics: Theory and Simulations, Prog. Astronaut. Aeronaut., Vol.159, pp.450-458.

1999, Karlin, I.V., Ferrante, A. and Öttinger, H.C., "Perfect entropy functions of the lattice Boltzmann method," Europhys. Lett., Vol.47, pp.182-188.

10.1209/epl/i1999-00370-1

2016, Frapolli, N., Chikatamarla, S.S. and Karlin, I.V., "Entropic lattice Boltzmann model for gas dynamics: Theory, boundary conditions, and implementation," Phys. Rev. E, Vol.93, p.063302.

10.1103/PhysRevE.93.06330227415382

2024, Kallikounis, N.G. and Karlin, I.V., "Particles on demand method: Theoretical analysis, simplification techniques, and model extensions," Phys. Rev. E, Vol.109, p.015304.

10.1103/PhysRevE.109.01530438366517

2024, Bhadauria, A. and Karlin, I.V., "High-speed flows with particles on demand: Boundary conditions," Comput. Fluids, Vol.271, p.106152

10.1016/j.compfluid.2023.106152

2015, Frapolli, N., Chikatamarla, S.S. and Karlin, I.V., "Entropic lattice Boltzmann model for compressible flows," Phys. Rev. E, Vol.92, p.061301(R).

10.1103/PhysRevE.92.06130126764625

2017, Dorschner, B., Chikatamarla, S.S. and Karlin, I.V., "Transitional flows with the entropic lattice Boltzmann method," J. Fluid Mech., Vol.824, pp.388-412.

10.1017/jfm.2017.356

2022, Tran, S.B.Q., Leong, F.Y., Le, Q.T. and Le, D.V., "Lattice Boltzmann method for high Reynolds number compressible flow," Comput. Fluids, Vol.249, p.105701.

10.1016/j.compfluid.2022.105701

1998, Karlin, I.V. and Succi, S., "Equilibria for discrete kinetic equations," Phys. Rev. E, Vol.58, p.R4053-R4056.

10.1103/PhysRevE.58.R4053

2020, Latt, J., Coreixas, C., Beny, J. and Parmigiani, A., "Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria," Philos. Trans. R. Soc. A, Vol.378, p.20190559.

10.1098/rsta.2019.055932833583PMC7333948

2023, Thyagarajan, K., Coreixas, C. and Latt, J., "Exponential distribution functions for positivity-preserving lattice Boltzmann schemes: Application to 2D compressible flow simulations," Phys. Fluids, Vol.35, p.126115.

10.1063/5.0175908

2023, Hosseini, S.A., Atif, M., Ansumali, S. and Karlin, I.V., "Entropic lattice Boltzmann methods: A review," Comput. Fluids, Vol.259, p.105884.

10.1016/j.compfluid.2023.105884

2008, Malaspinas, O. and Deville, M., "Towards a physical interpretation of the entropic lattice Boltzmann method," Phys. Rev. E, Vol.78, p.066705.

10.1103/PhysRevE.78.06670519256979

2015, Dorschner, B., Chikatamarla, S.S., Bösch, F. and Karlin, I.V., "Grad's approximation for moving and stationary walls in entropic lattice Boltzmann simulations," J. Comput. Phys., Vol.295, pp.340-354.

10.1016/j.jcp.2015.04.017

2023, Kaushal, S., Succi, S. and Ansumali, S., "Flow force calculation in the lattice Boltzmann method," Phys. Rev. E,Vol.108, p.045304.

10.1103/PhysRevE.108.04530437978617

2018, Succi, S., The Lattice Boltzmann Equation: For Complex States of Flowing Matter, Oxford University Press.

1995, Henderson, R.D., "Details of the drag curve near the onset of vortex shedding," Phys. Fluids, Vol.7, pp.2102-2104.

10.1063/1.868459

1994, Henderson, R.D., "Unstructured spectral element methods: Parallel algorithms and simulations," Ph.D. thesis, Princeton University.

2015, Kurtulus, D.F., "On the unsteady behavior of the flow around NACA 0012 airfoil with steady external conditions at Re=1000," Int. J. Micro Air Vehicles, Vol.7, pp.301-326.

10.1260/1756-8293.7.3.301

2012, Liu, Y., Li, K., Zhang, J., Wang, H. and Liu, L., "Numerical bifurcation analysis of static stall of airfoil and dynamic stall under unsteady perturbation," Commun. Nonlinear Sci. Numer. Simul., Vol.17, pp.3427-3434.

10.1016/j.cnsns.2011.12.007

2012, Khalid, M. and Akhtar, I., "Characteristics of flow past a symmetric airfoil at low Reynolds number: A nonlinear perspective," In Proceedings of the ASME 2012 International Mechanical Engineering Congress and Exposition, Vol.7, pp.167-175.

10.1115/IMECE2012-87389

2014, Karlin, I.V., Bösch, F. and Chikatamarla, S.S., "Gibbs' principle for the lattice-kinetic theory of fluid dynamics," Phys. Rev. E, Vol.90, p.031302(R).

10.1103/PhysRevE.90.03130225314388

1998, He, X., Chen, S. and Doolen, G.D., "A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit," J. Comput. Phys., Vol.146, pp.282-300.

10.1006/jcph.1998.6057

2007, Ansumali, S., Arcidiacono, S., Chikatamarla, S.S., Prasianakis, N.I., Gorban, A.N. and Karlin, I.V., "Quasi-equilibrium lattice Boltzmann method," Eur. Phys. J. B, Vol.56, No.2, pp.135-139.

10.1140/epjb/e2007-00100-1

Journal of Computational Fluids Engineering ISSN:1598-6071(Print) 3022-6252(Online) 한국전산유체공학회지

Preview

DEVELOPMENT AND VALIDATION OF A FLOW SOLVER BASED ON THE ENTROPIC LATTICE BOLTZMANN METHOD

ABSTRACT

MAIN

Fig. 1.

Stability regions depending on relaxation factor(𝛽)

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

Fig. 2.

D2Q21 3-level lattice velocity model

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

Fig. 3.

Simulation algorithm of in-house code based on the entropic lattice Boltzmann method

(24)

Fig. 4.

Simulation domain for a flow past 2D circular cylinder

(25)

Fig. 5.

(a) Drag coefficients and (b) relaxation factor for 2D circular cylinder at various Reynolds number

Fig. 6.

Simulation domain for a flow past 2D NACA0012 airfoil

Fig. 7.

Velocity field for NACA0012 airfoil at AoA = 24°

Fig. 8.

(a) Lift and (b) drag coefficients for NACA0012 airfoil at various angle of attack

Acknowledgements

References