1. 서 론

1970년대부터 항공우주공학 분야를 중심으로 전산유체역학(computational fluid dynamics, CFD) 기법들이 지속적으로 개발 및 확장되어 왔으며, 광범위한 산업분야에서 다양한 CFD 기법들이 활용되고 있다. 유체 흐름을 모델링하는 방법에 따른 전형적인 유동 모델링 체계는 크게 두 가지로 분류될 수 있다. 첫 번째는 유체 흐름을 입자들의 조합으로 고려하는 Molecular 접근법과 두 번째는 전통적 격자 기반의 Continuum 접근법에 기반한 수치해석 기법이다. Fig. 1은 Knudsen number(Kn)에 기반한 유동 모델링 접근법의 분류와 이들의 상호 관계를 보여주고 있다[1]. 여기서, Kn는 유체역학 및 열역학에서 유체의 연속체 거동과 분자적 거동 간의 상대적 중요성을 평가하기 위해 사용되는 무차원 값이며, 일반적으로 유체가 연속체로 고려될 수 있는지, 아니면 분자적 효과를 고려해야 하는지를 결정하는 지표로 활용된다. 물리적으로는 분자 간 충돌 없이 입자가 이동할 수 있는 평균 거리와 특성길이(물체의 크기, 유동 통로의 폭 등)의 비에 해당한다. Lattice Boltzmann Method(LBM)은 Lattice Gas Automata(LGA)의 통계적 잡음을 제거하고 더 우수한 갈릴레이 불변성(Galilean invariance)을 달성하기 위해 개선된 수치 유동해석 기법으로, 이미 오래전부터 제안된 개념이다[2, 3]. 또한, Kinetic theory와의 밀접한 연관성을 통한 높은 유연성 덕분에, LBM은 여러 물리적 현상을 모델링할 수 있도록 적합하게 조정될 수 있는 특성을 가진다. 입자 운동 Lagrangian 차분화 기법의 LBM은 유체 유동 방정식을 해결하는 데 있어 점차 중요한 방법으로 고려되고 있는데, 이 기법은 기체 운동론(Kinetic Theory of Gases)과 직접적으로 연결되어 있어 물리적 현상에 대한 심층적 이해와 더욱 세밀한 모델링을 가능하게 한다[4, 5]. 이러한 특징은 다상/다성분 유동(multiphase/multicomponent flow) 연구에서 사실상의 표준 도구로 많은 연구 그룹에 의해 채택되는 추세로 이어지고 있으며[6], LBM을 압축성 유동(compressible flow) 해석에 적용하는 연구도 수행되었다[7].

Fig. 1.

Classification of flow regimes and interrelationships of various flow models[1]

통상적인 CFD 기법과 비교했을 때, LBM 기법의 가장 큰 장점은 첫 번째로 복잡한 3차원 형상에 대해 전처리(pre-processing) 과정에서 해석자가 body-fitted 개념의 세밀한 격자를 직접 생성하는 과정이 필요하지 않다는 점과(CAD 형상 파일을 입력받아 필요한 3차원 Octree-based lattice structure(OLS)를 자동으로 생성하는 개념 활용), 두 번째로 기반 이론적 특성상 뭉툭한 물체나 고받음각 형상에서 발생하는 비정상(unsteady) 특성이 강한 박리 유동 현상에 대한 신뢰성이 높으며, 세 번째로 격자의 꼬임 현상이 거의 없이 물체의 동적 움직임을 고려한 유동 해석을 안정적으로 수행할 수 있다는 점으로 판단된다. 반면, LBM의 주요 단점으로는 외부 유동의 경우 정확도를 보장하기 위해 요구되는 particle resolution 수준이 기반 연구 사례의 부족으로 아직 명확하지 않아 보이며, 전형적인 Octree-based lattice structure 생성 요구도에 따라 물체 표면 부근에 배치되는 lattice 조밀도 증가를 위해서는 급격한 자유도 증가로 수치해석 소요 시간이 기존 연속체 기반 CFD 해석 대비 대폭 증가할 수 있다는 점이다.

압축성 고속 외부 유동 해석의 경우, LBM 기반 XFlow 프로그램을 활용한 국내외 연구 사례는 아직 드문 실정이며, 이를 활용한 주요 연구 사례들은 주로 저속 유동 영역에서 특정 모델이나 현상에 대한 해석 및 분석에 한정되어 있는 경향이다[8, 9, 10, 11]. 연속체 기반 CFD 기법은 과거 수십 년간 방대한 연구 사례와 경험 등을 통해 다양한 모델에 대한 유동 해석을 상당히 정확하게 수행할 수 있는 정보가 충분한 상황이라고 볼 수 있다. 반면, LBM 기반의 압축성 유동 해석은 전술한 장점이 있음에도 불구하고 아직은 활용 및 검증 사례가 충분하지 않아 더 다양한 모델에 대한 응용해석과 검토가 필요해 보인다. 이후 본문에 상세히 기술하였지만, LBM 기법은 meshless 유동해석 기법으로 분류되나 내부적으로는 수치해석 효율성 때문에 Octree lattice structure(OLS)에 기반하여 필요한 lattice structure를 생성하여 활용한다. 따라서 LBM 기법을 활용한 유동의 경우 가장 우선적이면서 민감한 사안 중 하나는 실용적이면서도 적절한 lattice resolution 수준의 결정이다. 이 논문에서는 전통적인 연속체 기반 CFD 기법과 비교했을 때, 일반적인 수준의 병렬처리 해석 서버에서 실용적인 신뢰도 수준의 결과를 도출하기 위해 요구되는 particle 분포 및 변수 설정 방법론과 압축성 유동 해석 검증 연구를 수행하였다. 이러한 연구는 항공우주공학 분야에서 LBM 기법의 활용성을 확보하기 위해 선행적으로 중요한 부분으로 고려되었다. 이에 LBM 기반 XFlow 소프트웨어를 활용하여 다양한 풍동시험 모델들에 대한 유동 해석을 수행하고 비교 및 분석 결과를 제시하고자 한다. 참고로, 다양한 슬로싱(sloshing) 현상에 대한 유동 해석과 비교 검증을 수행한 결과 LBM 기법의 효율성과 신뢰성이 기존 연속체 기반 CFD 기법에 비해 더 우수한 경향성을 확인하였으며, 이에 대해서는 향후 별도 논문을 작성할 계획이다.

2. 수치해석 기법

2.1 Lattice Boltzmann Method(LBM)

LBM에는 다양한 접근 방식이 존재하는데, 이는 nonlinear hyperbolic conservation law에 기반한 방정식 시스템을 해결하기 위한 일련의 방법론으로 간주할 수 있다. 모든 LBM 모델의 공통된 특징은 lattice discretization과 propagate- collide scheme에 기반한 time-stepping 모델을 채택하고 있다는 점이다. 격자에서 수행되는 입자의 전파 운동은 시간 간격(dt)과 이산적인 속도 집합(e→i,i=1,2,...,b)에 기반하여 계산되지만 각 입자 운동의 위치는 lattice site로 제한되게 된다. 이러한 속도 집합은 lattice를 생성하며, 각각의 lattice site에는 b개의 fi(r→,t)가 저장되는데, fi(r→,t)는 확률분포함수(probability distribution function, PDF)를 의미한다. XFlow의 3차원 lattice structure는 27개의 속도 방향(D3Q27)을 포함하는데 이는 통상적인 LBM 기법에 비해 더 높은 차원의 공간 이산화 기법에 해당한다[10]. LBM은 확률 분포 함수를 사용하여 유체 흐름을 mesoscopic scale로 고려하여 볼츠만 수송 방정식(Boltzmann transport equation)을 푼다.

이산화된 속도를 가지는 연속체 공간(continuum space)에서 볼츠만 수송 방정식은 다음과 같이 나타낼 수 있다[12].

여기서, e→i는 이산속도(discrete velocity)를 의미하며, Ωi는 질량과 선형 운동량을 보존하면서 충돌 이후의 상태를 계산하는 collision operator를 나타내며, Ωi = 0인 경우(입자 간 충돌이 없는 경우)는 단지 streaming 연산만 수행하게 된다. 볼츠만 수송 방정식은 collision-streaming scheme을 적용하여 lattice에 대해 다음과 같이 이산화(discretization)될 수 있다[11].

위 식에서 여기서 fi는 𝑖 방향의 확률분포함수, 𝑄는 속도 방향의 개수, x→는 위치벡터, e→i=(eix,eiy,eiz)는 이산 속도 방향벡터, 그리고 Ωi는 collision operator를 나타낸다. 첫 번째와 두 번째 항은 스트리밍 단계(streaming steps)라고 불리며, 우변은 collision operator을 포함한 충돌 단계에 해당한다. LBM의 stream-and-collide scheme은 연속적인 볼츠만 방정식의 이산 근사(discrete approximation)로 해석될 수 있다. Stream(또는 propagation) 단계는 이산 방향에 따른 입자 분포 함수의 advection 현상을 모델링하는 반면 대부분의 물리적 현상은 collision operator에 의해 모사 되는데, collision operator는 scheme의 수치적 안정성에도 큰 영향을 미친다. Stream-and-collide 접근법은 이산적인 속도 방향 집합을 갖는 Caretsian 점들로 구성된 lattice에서 정의된다. 이러한 lattice scheme은 일반적으로 DnQm으로 표기하는데 n은 문제의 차원을 m은 속도 방향의 개수를 의미한다. XFlow의 LBM solver는 D3Q27 lattice scheme과 central moment collision operator[13]를 기반으로 하며, 이는 표준 LBM 기법에 비해 수치적 안정성을 증대시킬 수 있는 것으로 알려져 있다. 또한 XFlow에서는 volumetric-based 접근법을 적용하기 때문에 확률분포함수 fi가 격자의 중심(barycenter)에 위치한다고 가정한다[14, 15].

LBM에서는 스트리밍(streaming) 단계 이후 생성된 미지의 확률분포함수(PDF) 값을 경계조건을 통해 재계산하며, 벽(wall)의 위치와 방향에 따라 이러한 미지 PDF 값이 달라질 수 있다. LBM에서 bounce-back 경계조건은 기하학적 벽에서 no-slip condition을 구현하는 데 사용되며, 벽 노드로 들어오는 PDF는 initial fluid node로 반사되어 미지 PDF 값에 할당된다. CFD 해석에서 주로 사용되는 주요 입구 및 출구 경계조건은 속도와 압력이다. 이는 일반적으로 알려져 있거나 쉽게 측정 가능한 데이터이기 때문이다. 이러한 접근법은 선형 시스템을 구성하며, 이는 첫 두 모멘트(질량과 운동량)의 보존에 기반한다. 이 선형 시스템은 스트리밍 단계 이후 생성된 미지의 PDF 값뿐만 아니라, 속도 또는 압력 경계조건이 부여되었을 때는 그 값을 직접적으로 할당한다. 또한, 경계에 수직인 입자 분포의 비평형(non-equilibrium) 부분에 대해서는 반사 규칙(bounce-back rule)을 적용한다. 이 접근법은 비교적 안정적이고 단순하지만, LBM의 비정상적이고 압축성 특성으로 인해 출구와 입구에서 수치적 반사 현상이 발생할 수 있다. 이러한 문제를 해결하기 위해 non-reflecting 경계조건이 도입되기도 한다[12].

2.2 Turbulence Modeling

난류는 유체역학 분야에서 가장 복잡한 메커니즘 중 하나로, 항상 일관된 결과를 제공하는 난류 모델링 기법의 개발은 CFD 분야에서 지속적인 도전 과제였다. 대부분의 유동 해석이 난류 조건에서 이루어지기 때문에, 정확한 난류 모델링 기법의 적용은 매우 중요한 문제로 여겨진다. Navier-Stokes 방정식은 강한 단순화와 가정을 고려하지 않는 한 해석적인 해를 가질 수 없으며, Direct Numerical Simulation(DNS)은 계산 자원의 한계로 인해 실제 공학적 사례에 대한 적용이 아직 비실용적이고 대부분 학문적 응용에 국한하여 사용되는 실정이다. 따라서, 산업계에서 요구되는 빠른 해석 시간을 충족시키기 위해 난류 효과를 적절한 정확도로 반영할 수 있는 실용적인 모델이 필요하다. 난류 모델링은 Joseph Boussinesq가 난류 응력을 시간 평균 유동으로 기술하기 위해 난류 점성(turbulent viscosity) 개념을 도입하면서 시작되었다. 이 난류 점성은 난류 와류에 의해 유발되는 운동량 전달을 모델링하며, 이후 Prandtl의 혼합 길이 개념이 보완되어 벽면 효과를 효율적으로 고려할 수 있게 되었다. 이러한 가정을 바탕으로, Reynolds-Averaged Navier-Stokes(RANS) 방정식을 해결하기 위한 다양한 난류 모델들이 개발되었다. 가장 널리 사용되는 모델로는 Spalart-Allmaras, 𝑘−ϵ, 𝑘−𝜔 및 SST(Shear Stress Transport) 모델 등이 있다. RANS 모델은 낮은 계산 비용과 합리적인 정확성 덕분에 현재 산업계에서 폭넓게 적용되고 있다. 그러나, RANS 모델은 경험적 모델이기 때문에 많은 상수를 포함하며, 유동 상황에 따라 난류를 올바르게 모델링하기 위해 이러한 상수를 신중하게 조정해야 하는 단점이 있다. 각 모델은 서로 다른 상수 집합을 포함할 수 있으며, 난류에 따른 분리 예측이 크게 달라질 수 있으므로 문제에 따라 완전히 다른 해를 예측할 가능성도 있다. 따라서, 적절한 모델의 선택은 고려하는 응용 분야에 크게 의존하며, 제품 설계에 큰 변화가 있는 경우 이미 설정 및 검증한 난류 모델링 기법에 불확실성이 초래될 가능성도 있다. 또한, RANS 난류 모델은 모든 와류 규모를 구별 없이 모델링하며, 시간 평균 방정식에 기반하고 있어 난류에 의해 유발되는 순간적인 변화를 모두 필터링한다. 따라서 이는 주로 정상상태(steady state) 유동 해석에 한정되며, Kolmogorov가 설명한 실제 난류 에너지 스펙트럼을 완전하고 정확히 포착할 수 없다. 이러한 제약은 특히 aeroacoustics 분야에서 중요한 단점으로 작용하게 된다.

난류 모델의 중간 단계로 분류되는 모델은 1964년 Joseph Smagorinsky에 의해 처음 제안되었다. 이는 해석되지 않은 서브 격자(sub-grid) 규모에서 난류 점성을 적용하는 방법을 일반화한 것이다. 이러한 난류 모델링 접근법은 Large Eddy Simulation(LES)이라고 불리며, RANS와 DNS 사이에서 좋은 중간 해결 방안이 된다. LES는 대부분의 난류를 직접적으로 해석하면서도 더 작은 규모의 난류를 모델링하여, Fig. 2에 나타난 것처럼 계산 시간을 현실적인 수준으로 유지할 수 있도록 한다. LES 난류 모델은 Navier-Stokes solver와 함께 사용될 때 계산 비용이 RANS 모델에 비해 상대적으로 높기 때문에 널리 활용되지는 않고 있다. 그러나 LES 모델은 고급 응용 분야나, RANS 모델이 정확도가 현저하게 떨어지는 고난류 및 유동박리 현상이 심한 transient flow에 대한 정교한 해석이 요구되는 경우 한정적으로 사용되는 경향이 있다. 따라서 난류의 시간 변동성이 심한 상황에 대한 연속체 기반 CFD 해석은 RANS와 LES 모델의 장점을 결합하여 실용성이 우수한 Detached Eddy Simulation(DES) 기법이 많이 활용되고 있다.

Fig. 2.

Comparison of turbulence modeling approaches

유체의 비선형적인 특성 때문에 Navier-Stokes 방정식 시스템의 해를 구하려면 원래 시스템의 이산화된 형식에 적합한 선형화가 필요하다. 유한요소법과 유한체적법에서는 Eulerian 수학적 정식화가 일반적인 접근법이지만, 이는 운동량 방정식의 non-linear advection term의 수치해석에 많은 어려움을 유발하는 특성이 있다. XFlow의 난류 모델링 접근법은 기본적으로 LES에 기반하며, Wall-Adapting Local Eddy(WALE) 점성 모델을 도입하여 일관된 local eddy-viscosity와 near-wall behavior를 제공할 수 있도록 한다[16]. 또한 XFlow에서 사용하는 Cartesian lattice structure는 LES 난류 모델에 매우 적합한데, 난류가 경계층 외부에서 등방성(isotropic) 특성을 가지며 LES 난류 모델은 비율이 균형 잡힌 셀 형상을 기본적으로 요구하기 때문이다. 하지만, lattice structure의 등방성은 경계층을 해석하기 위해 비합리적으로 많은 요소 수가 필요할 수 있으므로 XFlow에서는 이러한 상황을 대비하기 위해 Wall-Modeled LES(WMLES) 벽함수(wall function) 모델을 포함하고 있다[12]. 이는 LES의 효율성과 계산 비용을 최적화하기 위해 벽면 근처의 난류를 벽함수로 처리하는 기법이다. LES는 전체 난류 스펙트럼을 직접 계산해야 하므로 경계층에서의 작은 와류까지 해석하려면 과도한 계산 자원이 요구되는 문제가 있는데, WMLES 기법은 벽함수 또는 근사 모델을 사용해 경계층 영역을 처리하며 벽에서 멀리 떨어진 영역만 LES 기법을 활용하므로 난류 모델 수치해석의 효율성을 상당히 증가시킬 수 있다. 이 기법의 단점으로는 벽 모델의 정확도에 따라 해석 결과의 신뢰도가 달라질 수 있으며, LES 영영과 벽모델 간의 상호 작용이 정확히 구현되려면 여전히 상당한 수준의 격자 조밀도가 요구되며, 경계층 유동이 매우 얇거나 불안정한 경우 벽모델이 난류를 제대로 표현하지 못할 가능성이 있다.

2.3 Octree Lattice Structure(OLS)

XFlow에서 사용되는 격자 구조는 D3Q27 형식으로, 옥트리(octree) 구조로 구성되어 있다. 이 격자 구조는 유한 체적법 또는 일반적인 LBM 코드보다 개별 요소당 자유도가 더 많으며, 4-th order 공간차분 정확도를 가진다. Octree lattice structure(OLS)는 유체 영역의 서로 다른 위치에서 서로 다른 공간 해상도를 가지는 비균일(non-uniform) 분포로 생성될 수 있다. 공간상의 격자 요소의 크기는 ∆x=H/2L의 세분화(refinmement) 관계식을 따르며, H는 계산영역(computational domain)의 최대 길이를 L은 세분화 수준(refinement level)을 의미한다. XFlow의 전처리기는 입력된 기하학적 형상과 각 형상에 대해 사용자가 지정한 격자 해상도 및 far-field 해상도를 기반으로 초기 lattice structure를 생성한다. 사용자는 구, 상자, 원기둥 등과 같은 별도 정의 영역을 설정하여 다른 lattice 조밀도 임의의 영역을 세밀화할 수 있는데, 사용된 서로 다른 공간 스케일은 계층적으로 배열된다. 각 레벨은 이전 레벨보다 공간 및 시간 스케일이 절반인 문제를 해결할 수 있도록 하는데, 이는 매우 효율적인 방식으로 유체 영역 전체에서 ∆x/∆t 비율을 유지하므로 모든 위치에서 일정한 CFL 조건과 음속이 보장된다. 따라서 local time stepping 접근법에 매우 유리하고 유체 영역의 각 lattice 크기에 항상 적합한 시간 스텝을 적용할 수 있는 장점이 있다[12].

XFlow 전처리기에서 생성된 초기 lattice structure는 여러 설정 기준에 따라 고전적인 adaptive mesh CFD 기법과 유사하게 해석 과정에 따라 동적으로 변경될 수 있다. 움직이는 기하학적 형상으로 인해 계산 영역이 변경되는 경우, 격자는 각 시간 스텝마다 기하학적 형상의 새로운 위치를 따라가도록 동적으로 세밀화될 수 있다. 또한, 물리적인 유동 특성 변화를 반영한 적응적 세밀화 기준도 반영할 수 있다. 와도(vorticity) 수준에 기반한 세밀화 알고리즘은 높은 와도 변동성을 가지는 후류 영역을 동적으로 세밀화하는 데 효과적으로 활용될 수 있다. Wake refinement 기능은 local dimensionless vorticity가 임계 ω*=(ω/(∆x/∆t))2 값보다 클 때 활성화되도록 설정할 수 있다. 또한 wake refinement가 적용되는 거리를 별도로 설정할 수 있는데, 물체로부터 세밀화가 적용될 lattice node까지의 범위가 한정되게 된다. 이러한 조건을 적용하기 위해 물체로부터 격자 노드까지의 거리 계산은 Manhattan distance 알고리즘[17]이 적용된다. 기타 자유표면유동(free surface flow) 및 다상유동(multi-phase flow) 현상의 경우, 자유표면 또는 상 경계면 등에 대해서도 dynamic lattice refinement 적용이 가능한 특성을 가진다[11].

3. 해석결과 및 검토

LBM 유동해석 기법은 meshless 기법으로 분류되지만, 전통적인 CFD 모델링과 마찬가지로 물체 면 부근에서 Octree lattice structure(OLS)의 조밀도 수준에 따라 결과에 차이가 발생할 수 있으며, 정확한 해석을 위해서는 다양한 모델에 대한 유동 해석과 정량적인 평가를 통한 경험 축적이 요구된다. 따라서 본 연구에서는 LBM 기법을 활용한 유동 해석 시 가장 우선적인 민감도 사항이 되는 lattice structure의 조밀도 수준에 대한 보조지표로 Element Resolution Ratio(ERR=ReferenceLength/TargetResolvedScale)를 계산하여 활용하였다. 여기서, 기준길이(reference length)는 해석 대상의 크기와 관계된 값이며, Target Resolved Scale(TRS)은 물체 표면에 분포되는 lattice의 최소 크기 설정에 관계되는 값이다.

3.1 F-5 Wing Model

F-5 주익 축소 풍동시험 모델의 날개 단면은 NACA 64A004.8 에어포일이고 두께비는 4.8% 이다[18]. 풍동시험 모델 날개의 참조면적(reference area, Sref)은 0.260621 m²이며, 날개 뿌리에서 코드 길이(c_r)는 0.6396 m 이다. Fig. 3은 LBM 유동해석을 위해 생성한 octree lattice 분포 예를 보여주고 있는데, near wall refinement 옵션을 적용하여 생성한 경우로 해석 영역 크기는 32 m(x) × 4 m(y) × 32 m(z) 이다. 이 모델에 대한 유동해석 및 평가를 위해 다양한 lattice 크기에 대한 해석이 수행되었으나, 지면 관계상 대표적인 lattice 분포만 예시로 제시하였다. F-5 날개에 대한 LBM 유동해석은 물체 표면에서 lattice structure 조밀도 크기에 대한 다양한 테스트 결과 중 shape refinement의 TRS 크기가 4 mm인 경우(ERR = 639.6 mm/4 mm = 160)와 3 mm인 경우(ERR = 639.6 mm/3 mm= 213)에 대한 풍동시험 비교 결과만 제시하였다. 참고로 정육면체 요소에 기반하는 Octree lattice 생성 알고리즘 특성상 ERR 값이 160인 경우 날개 위/아래 면에서 코드 방향으로 배치되는 평균 요소 개수는 160개, 스팬 방향으로는 약 155(622.6 mm/4 mm)개로 이론적으로는 날개의 윗면에만 약 24,800개의 요소가 배치된다. 이는 연속체 기반 CFD 해석 격자에 비해 상대적으로 매우 조밀한 격자 배치에 해당한다. 참고로 Octree lattice 생성 알고리즘이 사용자에게는 편리하지만, 큰 단점 중 하나는 span 방향 또한 등방성 육면제 요소 배치 특성 때문에 연속체 기반의 CFD 격자 생성과 같이 조밀도 함수를 별도로 적용하여 조정할 수 없다는 점이다. 이는 세장비가 큰 형상의 경우, 기본적인 정확도 유지를 위해 코드 방향으로 요구되는 lattice의 개수가 많은 경우 표면에 과도한 수의 lattice가 표면에 배치되는 문제를 유발하게 된다. LBM 해석에서 주로 사용하는 OLS의 경우 표면 격자 크기를 1/2로 축소하고 shape refinement 옵션을 적용하는 경우 생성되는 요소 개수와 해석 시간이 약 4배 이상 증가하는 특성을 가지기 때문에 참고문헌[19]의 로켓 형상과 같이 세장비가 큰 모델에 대해 동적 pitching motion을 고려한 해석까지 고려하는 경우 해석모델의 실용성이 급격하게 저하되는 단점이 있다.

Fig. 3.

Geometric configuration of F-5 wing wind-tunnel test model with NACA 64A004.8 airfoil

F-5 날개 유동해석을 위한 경계조건은 유입 경계면은 inlet velocity로 유출 경계면은 gauge pressure outlet으로 설정하였으며, 대칭면과 측면 원방 경계면은 free slip 경계조건을 상하 원방 경계면은 periodic 조건을 설정하였다. 날개면 경계조건은 비점성 유동 해석의 경우는 free-slip 조건을 점성 유동 해석의 경우는 non-equilibrium enhanced wall function 조건을 부가하였다. 마하 0.8 및 마하 1.35 조건에 대해서는 비점성(inviscid) 유동해석 결과를 제시하였으며, 날개 면에서 수직 충격파가 강하게 발생하는 마하 0.95 조건에 대해서는 Wall-Modeled Large-Eddy Simulation(WMLES) 기법을 적용한 해석 결과를 제시하였다. Fig. 4는 F-5 날개 풍동시험 모델에 대해 유입 마하수 0.6, 받음각 0° 조건에서 비점성 압축성 LBM 유동해석으로 구한 압력계수를 풍동시험과 비교한 결과를 보여주고 있으며, 전반적으로 LBM 해석 결과가 풍동시험 결과와 잘 일치하는 경향을 나타내고 있다. 이 해석의 경우 LBM 모델의 물체 표면 요소 크기는 TRS 값이 4 mm인 경우로 ERR 값은 약 160에 해당하며, 계산 영역 내의 전체 요소 개수는 약 42만 개이다. Fig. 5는 유입 마하수 0.95, 받음각 0° 조건에서 점성 압축성 유동해석을 수행하여 구해진 날개면 압력계수를 풍동시험과 비교한 결과이다. 이 해석의 경우 물체 표면에서 TRS 값은 3.0 mm인 경우로 ERR 값은 약 210에 해당하여, 계산 영역 내의 전체 요소 개수는 약 72만 개다. LBM 유동해석은 기본적으로 비정상(unsteady) 해석에 기반하기 때문에, 이 경우 양력계수가 수렴된 결과를 도출하는데 물리적으로 약 0.6초 동안의 해석이 소요되었다. XFlow의 automatic time-stepping 옵션을 적용하여 수행한 비정상 유동해석은 Intel Xeon CPU E5-2687W v2 3.40 GHz(16 cores, 128 GB RAM) 서버로 약 4시간이 소요되었다. LBM 해석 결과가 풍동시험 결과를 전반적으로 잘 일치하나, 아래면 앞전 부근과 87.5% 스팬 위치에서 50~70% 코드 영역에서는 다소 차이를 나타내고 있다. 이러한 차이는 연속체 기반 CFD 해석에서도 초기 생성한 격자의 품질에 따라 발생할 수 있는 경우로 LBM 해석의 경우도 물체 표면 및 인근의 요소 분포 수준에 따라 차이가 발생할 수 있는 부분에 해당하며, 이는 수직 충격파 발생 유동조건에서 요소 분포가 다소 부족하기 때문으로 판단된다. Fig. 6은 유입 마하수 1.35, 받음각 0° 초음속 조건에서 비점성 LBM 유동해석과 풍동시험 압력계수 비교 결과를 나타내고 있다. 마하 0.8의 경우와 마찬가지로 물체 표면에서 TRS 값이 4 mm인 경우로 ERR 값은 160에 해당하는 경우로 초음속 유동해석의 경우 마하 0.95 조건보다 오히려 유동의 선형 특성이 강하기 때문에 상대적으로 낮은 조밀도의 격자 분포로도 풍동시험 결과와 전반적으로 잘 일치하는 결과를 얻을 수 있음을 확인할 수 있다.

Fig. 4.

Comparison of sectional pressure coefficient for the F-5 wing model at M = 0.8 and α = 0°

Fig. 5.

Comparison of sectional pressure coefficient for the F-5 wing model at M = 0.95 and α = 0°

Fig. 6.

Comparison of sectional pressure coefficient for the F-5 wing model at M = 1.35 and α = 0°

3.2 Standard Dynamic Model(SDM)

본 절에서는 LBM 기법의 받음각 변화에 따른 양력 및 모멘트 계수 예측 성능을 평가하기 위해 풍동시험 결과가 제시되어 있는 Standard Dynamic Model(SDM) 형상[20]에 대한 3차원 CAD 모델링과 유동해석을 수행하였다. SDM 모델은 Leading-Edge Extension(LEX)을 가지는 주익, 동체, 수평 및 수직미익 및 ventral fin을 포함한 단순화된 전기체 F-16 전투기 형상에 해당하며, 풍동시험용 축소 모델의 형상은 Fig. 7과 같다. 이 모델의 경우 주익, 수평 및 수직미익의 단면이 wedge-flat-wedge 형상으로 되어 있다. 시험 모델의 폭은 609.6 mm(24 in), 동체 길이는 942.9 mm, 동체 중심에서 수직미익 끝단까지 높이는 214.5 mm(8.446 in) 이다. 그리고 주익과 수평미익 앞전의 후퇴각은 40°, 수평미익 하반각(anhedral angle)은 10° 이고, 수직미익 앞전의 후퇴각은 47.5°, 동체 중심 수직면에 대한 ventral fin 경사 각도는 25° 이다. Fig. 8은 SDM 풍동시험 모델 사진과 유동해석을 위해 생성한 3차원 CAD 모델 형상을 나타내고 있다.

Fig. 7.

Standard Dynamic Model(SDM) layout[19]

Fig. 8.

The SDM wind-tunnel support system(left) and generated CAD model(right)

Fig. 9는 SDM 형상에 대한 LBM 유동해석을 위해 생성한 OLS 격자를 보여주고 있는데, 해석 영역 크기는 24 m × 5 m × 24 m, 최대 resolved scale은 0.256 m, 물체 표면에서 목표 최소 요소 크기인 TRS 값으로 4 mm(ERR = 60.96 mm/4 mm = 152)와 1 mm(ERR = 60.96 mm/1 mm = 610)인 두 가지 조건에 대해 shape refinement 옵션을 추가로 적용하여 생성하였다. 이 경우, 생성되는 총 요소 개수는 TRS가 4 mm 일 때 1,024,116개이고, TRS가 1 mm인 경우는 9,640,185개로 요소 수가 약 8배 이상 증가하였다. 유동 흐름 방향은 x-축 방향이고, 유입 경계면은 inlet velocity 및 유출 경계면은 gauge pressure outlet으로 설정하였다. 해석모델의 측면은 y-축 방향으로, 측면의 경계조건은 free slip 경계조건을 부가하였다. 이 모델의 경우 공력계수 산정을 위한 참조면적은 0.1238 m², 기준길이는 0.6096 m이며, 피칭 모멘트 중심은 주익 root chord의 앞전에서 후방으로 0.3944 m인 지점에 위치한다. 풍동시험 유입 마하수가 0.6인 조건으로 강한 충격파 효과는 기대되지 않으므로 XFlow의 thermal model 옵션은 isothermal 모델을 적용하였고, 풍동시험 결과가 받음각 범위에서 선형 특성을 나타내고 있어 참고문헌[21]과 비교평가를 위해 비점성 유동해석을 수행하였다.

Fig. 9.

Generated lattice structures with near-wall refinement for the SDM shape(Total elements: 1 M(Left, ERRB=152, TRS=4 mm), 9.6 M(Right, ERRB=610, TRS= 1 mm)

Fig. 10은 유입 마하수 0.6에서 받음각 –2.5°에서 12°까지 범위에 대해 LBM 유동해석을 수행하고 풍동시험 데이터와 참고문헌[21]의 해석결과로 제시되어 있는 양력 및 피치 모멘트 계수에 대한 상호 비교 결과를 보여주고 있다. LBM 해석을 위한 lattice structure의 요소 수가 100만 개(TRS=4 mm, ERR=152) 수준에 해당하는 해석 결과의 경우 받음각 6° 이상에서 풍동시험 양력계수보다 더 큰 값을 예측하는 경향이 있으며, 피치 모멘트 계수의 경우는 전반적으로 약 0.02 정도의 편향된 차이를 보이고 있다. Lattice structure의 요소 수가 960만 개(TRS=4 mm, ERR=152) 수준의 해석 결과는 풍동시험 양력계수와 더 잘 일치하는 결과를 보이고 있다. 피치 모멘트 계수의 경우는 받음각 2°~ 6° 범위에서는 풍동시험 결과와 잘 일치하고 있으나, 2° 이하의 받음각에서는 다소 편차가 나타나고 있다. 참고문헌[20]에서는 SDM 형상에 대해 연속체 CFD 기반의 Euler solver를 활용하여 세 가지의 서로 다른 조밀도 격자(HC: 502,497 cells, HM: 1,144,774 cells, HF: 1,579,152 cells)에 대한 유동해석을 수행하고 비교하였는데, 피치 모멘트 계수의 경우 참고문헌의 HC 격자는 비교적 큰 차이를 보였고, 모든 격자에 대해 2° 이하의 받음각에서는 본 연구에서와 유사한 경향의 차이가 나타났다.

Fig. 10.

Comparison of lift(left) and pitching moment coefficients(right) for the F-16 SDM model at various angles- of-attack conditions(M = 0.6)

3.3 Generic Fighter Model

Fig. 11은 Generic Fighter 풍동시험 모델 형상[22]을 보여주고 있다. 그림의 크기 치수는 실제 대상 항공기 크기와 유사하게 10배 확대된 값이며, 실제 풍동시험은 Fig. 11의 좌측과 같이 1/10 scale 모델을 대상으로 수행되었다. 주익은 NACA 64A204 에어포일(4% 두께비)로 구성되며, 날개 뿌리(wing root) 각도는 –1°, 날개 끝단 각도는 –2.4° 이다. 풍동시험 결과와 비교하기 위한 유동해석 조건은 마하 0.85 및 받음각 2.12° 이다. Fig. 12는 1/10 scale 형상에 대한 생성된 OLS 격자를 동체 및 날개 주변만 나타낸 것으로, 물체면에서 TRS를 2 mm(ERR = 192) 설정하고 near wall refinement 옵션을 추가 적용하여 생성한 경우로 총 요소 수는 약 425만 개인 경우이다. 이 모델의 경우 lattice 조밀도에 대한 수렴성을 평가하기 위해 물체 표면 요소 크기를 1/2로 축소한(TRS = 1 mm, ERR = 385) 경우에 대한 해석도 수행하였는데, 이 경우 총 요소 수는 약 1,430만 개이다. 전체 계산 영역 원방경계 크기는 20 m(x) × 4 m(y) × 15 m(z) 이고, 유동 흐름은 x-축 방향이다. 유입 경계면은 inlet velocity 및 유출 경계면은 gauge pressure outlet 조건을 부가하였고, 상하 및 좌우 원방 경계에는 free slip wall 조건을 부여하였다.

Fig. 11.

Geometric configuration of a generic fighter wind tunnel test model(10X scale-up size)[22]

Fig. 12.

Generated lattice structures with near-wall refinement for the generic fighter model

이 해석의 경우 유입 마하수가 0.8 이상인 0.85이므로 비선형 압축성 효과를 더 정확하게 고려할 수 있도록 XFlow의 thermal model은 supersonic 옵션을 적용하였다. 또한 난류 점성 효과를 고려하기 위해 WALE 기법과 경계층에는 non-equilibrium enhanced wall function 모델을 적용하였다. 또한 해석모델에 받음각 조건을 부여하기 위해 참고문헌에서 주어진 무게중심 좌표를 기준으로 CAD 형상을 받음각만큼 회전하여 OLS 격자를 생성하였다. Fig. 13은 유입 마하수 0.85, 받음각은 2.12°인 조건에서 날개 스팬 방향으로 59% 및 85% 지점에서 단면 압력계수를 비교한 결과이다. 그림에서 연속체 기반 CFD 해석 결과는 참고문헌[22]에서 인용한 것이며, LBM 해석모델에 받음각 조건을 부여하기 위해 참고문헌에서 주어진 무게중심 좌표를 기준으로 형상을 받음각만큼 회전하여 해석을 수행하였다. Fig. 13의 우측 결과는 85% 스팬 지점에서 단면 압력계수 비교 결과를 나타낸다. LBM 해석 결과 ERR 190(4.25 M elements) 및 385(14.3 M elements)인 경우 모두 풍동시험 결과를 전반적으로 잘 예측하고 있음을 확인할 수 있다. 하지만 참고문헌의 CFD 결과는 날개 아래면 20% 코드 지점 부근에서의 압력계수 분포에 다소 차이를 나타내고 있다. Fig. 13의 좌측 결과는 59% 지점에서 단면 압력계수 비교 결과를 나타낸다. LBM 해석의 경우 날개 윗면은 풍동시험 결과와 거의 일치하였으나, 날개 아래면에서는 ERR 385 조건에서도 다소 차이를 보이고 있다. 아직 차이의 명확한 원인은 파악하지 못했으나, 참고문헌에서 풍동시험 모델과의 일치를 위해 날개의 CAD 형상을 약간 변경하였다는 언급을 바탕으로 유추해 볼 때, 본 연구에 사용된 Generic Fighter CAD 모델의 날개 형상에 미소한 편차가 존재하거나, 날개 면에서 lattice structure의 조밀도가 아직 충분하지 않은 원인일 가능성이 있다. 참고로 ERR 385 조건에 대한 XFlow 해석은 AMD Ryzen Threadripper PRO 5995WX CPU(2.7 GHz, 64 cores)와 512 GB RAM을 갖춘 병렬 연산 서버에서 수행되었으며, 64코어를 사용하여 수렴된 결과를 도출하는 데 약 26시간이 소요되었다.

Fig. 13.

Comparison of sectional pressure coefficient for Generic fighter figuration model at M = 0.85 and α = 2.12°

3.4 Wing/Pylon/Finned Store Model

Fig. 14는 날개 앞전 후퇴각이 45°이고, 단면 형상이 NACA 64A010 에어포일(10% 두께비)인 날개에 파일런 및 장착물이 추가된 형상을 나타내며, 실제 풍동시험 모델은 그림의 치수 대비 1/10 scale에 해당한다[23]. Fig. 15에는 이 모델에 대해 near wall refinement 옵션을 부가하여 생성한 OLS 격자를 제시하였으며, 물체 표면에서 TRS는 2 mm 이다. 계산 영역의 전체 크기는 25 m(x) × 3 m(y) × 15 m(z) 이고, 초기 생성된 총 요소 수는 약 252만 개로 날개 뿌리의 코드 길이를 참조 길이로 고려했을 때 ERR 값은 약 380에 해당한다. 점성 압축성 유동 해석을 위해 WALE 기법을 적용하였으며, 날개면 경계조건은 non-equilibrium enhanced wall function 조건을 부가하였고 XFlow solver의 thermal model은 비선형 충격파 효과를 고려하기 위해 supersonic으로 설정하였다.

Fig. 14.

Geometric configuration of a generic fighter wind tunnel test model(10 times scale-up size)

Fig. 15.

Generated lattice structures with near-wall refinement for a wing/pylon/finned store model

Fig. 16은 날개 및 장착물 단면에서 압력계수를 풍동시험 및 참고문헌[24]의 USM3D CFD solver(비정렬 overset 격자, Cell 수 680만 개) 해석 결과와 비교한 것이다. 참고문헌에서는 장착물 후방 끝단에 지지봉(sting) 형상을 추가로 고려하였으나, 본 해석에서는 Fig. 14와 같이 이를 고려하지 않고 해석을 수행하였다. Fig. 16의 좌측 결과는 날개 40.8% 스팬 위치에서 압력계수를 비교한 결과로 LBM 해석 결과가 풍동시험 결과와 유사한 경향을 보이나 파일런이 위치한 날개 아래면에서 충격파 발생 이후 중간 영역에서 압력계수(절대값)가 풍동시험에 비해 다소 낮은 수준을 나타내고 있다. 파일런과 날개가 접하는 부분에서 복잡한 유동현상이 발생하는 이유는 Fig. 17에 제시한 경계층 외부 표면 속도분포를 보면 파악할 수 있다. 날개 아래면에서는 파일런 앞전 부근에서 발생한 1차 충격파가 날개와 간섭을 유발하고 있으며, 이후 파일런 뒷전 부근에서 날개 스팬방향으로 뻣어 나가는 형태로 발생하는 강한 수직충격파 전후로 압력계수가 급격히 상승하게 된다(그래프에서는 y-축이 –Cp이므로 하락으로 나타남). Fig. 16의 우측 결과는 장착물 단면에서 압력계수를 비교한 결과로 LBM 해석 결과가 풍동시험과 전반적으로 잘 일치하고 있음을 확인할 수 있다.

Fig. 16.

Comparison of sectional pressure coefficient for wing/pylon/finned store model at M = 0.95 and α = 0°

Fig. 17.

Velocity distribution contour for the wing/pylon/finned store model at M = 0.95 and α = 0°

이 논문에서는 총 4종의 다른 풍동시험 날개 모델에 대한 압축성 LBM 유동해석을 수행하고 이를 평가하였고, 1,000만개 이하의 격자 수준에서 일부 특이 사항을 제외하면 전반적으로 풍동시험과 비교적 일치하는 결과를 도출할 수 있음을 보였다. 연구팀 내부적으로는 다양한 다른 형상의 날개에 대한 해석 및 평가가 있었는데 결론에 앞서 LBM을 활용한 유동해석 및 검증과 관련된 교신저자의 경험을 기술하고자 한다. LES 모델에 기반한 LBM 기법이 연속체 기반 CFD 해석에 비해 정성적으로나 정량적으로 강점을 가지는 대표적인 조건은 박리가 심한 비압축성 고받음각 유동의 경우이며 이는 기존 참고문헌[10]에도 나타나 있다. 특이 사항은 압축성 유동의 경우 아직 항공우주분야 응용이 활발하지 않아서 3차원 날개에 대한 LBM 공력해석 검증 연구 사례가 매우 드문 편이다. 교신저자의 연구실에서는 과거 LBM 기반 XFlow를 활용한 NACA 0012 에어포일(12% 두께비), RAE Tailplane(두께비 10%, 앞전 후퇴각 50.2°) 및 ONERA M6 Wing(두께비 10%, 앞전 후퇴각 30°) 모델과 등과 같이 두께비가 큰 에어포일을 가지는 형상에 대해 강한 수직 충격파가 발생하는(천음속 플러터 발생 위험성이 높은) 마하수 영역(M = 0.85~0.95)에 대한 압축성 유동해석을 수행하고 기존 풍동시험 결과와 압력계수를 비교 검증하는 연구를 수행하였다. 유동 가시화를 통한 정성적 특성은 대체로 유사한 경향을 보였으나 연속체 기반 CFD 해석 격자에 비해 상당히 많은 요소를 생성하고도 정량적인 압력계수에는 여전히 큰 편차가 존재할 수 있음을 경험하였다(단, 두께비가 큰 경우도 수직 충격파가 강하게 발생하지 않는 마하수 조건에서는 풍동시험 결과와 비교적 잘 일치함). 이러한 특성이 LBM 기반의 XFlow solver에 한정된 특성인지 아니면 LBM 기법의 공통된 특성인지는 아직 확실하게 파악하지는 못하였다. 다만, 과거 XFlow을 활용한 관련 해석 및 비교 경험에 기반해 볼 때 연구실 내부적으로는 양력면(lifting surface) 두께비가 5% 이하인 경우만 실용적인 OLS 격자 수준에서 정량적으로 유용한 해석 결과를 기대할 수 있을 것으로만 판단하고 있다. 향후 천음속 영역에서는 정적으로 고정된 물체 조건 이외에도 동적 움직임을 고려한 다양한 조건에 대한 비교 검증 연구가 추가로 수행될 필요성이 크다는 판단이다. 그리고 LBM 기법은 기존 연속체 기반의 CFD 기법과는 solver의 특성이 다르고 OLS 개념의 격자를 내부적으로 자동 생성 및 적용하기 때문에 전통적인 CFD 유동해석과 마찬가지로 편차를 최소로 하기 위해서는 검증 가능한 다양한 모델에 대한 OLS 격자 민감도 분석이 선행적으로 중요한 사항이다.

4. 결 론

본 논문에서는 Lattice-Boltzmann Method(LBM) 기반의 XFlow 프로그램을 활용하여 기존 풍동시험 결과가 제시되어 있는 다양한 날개 모델에 대한 압축성 유동 해석과 비교평가를 수행하였고 이를 통해 다음과 같은 결론을 얻을 수 있었다. 첫째, 격자 생성 측면에서 LBM 유동해석 기법에서 주로 활용하는 Octree-based lattice structure(OLS) 기법은 단순히 물체 CAD 형상을 입력하는 것만으로 복잡한 3차원 형상에 대한 고품질 격자를 효율적으로 생성할 수 있음을 확인하였다. 둘째, 압축성 유동 해석 측면에서, 주로 전투기 날개 설계에 활용되는 두께비 5% 이하의 날개 형상에 대해 강한 천음속 및 초음속 유동 조건에서도 기존 풍동시험과 대체로 일치하는 결과를 얻을 수 있음을 확인하였다. 셋째, 3차원 날개에 대한 압축성 유동 해석 시 본 논문에서 정의한 ERR 기준 200~600 범위의 particle resolution으로 실용적 신뢰도 수준의 공력해석 결과를 도출할 수 있음을 확인하였다. 다만, 날개의 두께비가 5% 이상인 날개 형상이나 천음속 영역에서 강한 수직 충격파가 발생하는 압축성 유동 조건에 대해서는 향후 다양한 추가 검증 연구가 필요하다. 본 연구의 결과는 향후 LBM 기반 압축성 유동 해석 코드 개발 및 관련 응용연구에 유용한 자료로 활용될 수 있을 것으로 판단된다.