Development and Application of an Accurate Least Squares Method Based on the Principle of Superposition for Computational Fluid Dynamics = 전산유체역학을 위한 중첩 원리 기반의 정확한 최소제곱법 개발 및 적용

다국어 초록 (Multilingual Abstract)

Over the years, enhancing the accuracy of the gradient operator, which is regarded as a key component of differential operators, has remained a fundamental research objective in mathematics, science, and engineering. For instance, it finds wide-ranging applications in research areas such as inverse problems, optimal design, Partial Differential Equation (PDE) analysis, regression & interpolation, signal & image processing, manifold analysis, and Artificial Intelligence (AI). In particular, the Least Squares Method (LSM), a representative approach for gradient estimation, has long been utilized in the field of computational fluid dynamics (CFD). In recent years, its utility has extended beyond gradient reconstruction to the LSM-based spatial discretization in meshless method.

However, in terms of accuracy, LSM still faces limitations in its application to spatial discretization and gradient reconstruction, particularly for boundary layer problems involving complex geometries with high Aspect Ratios (AR), where it becomes difficult to accurately estimate the gradients of physical quantities. In addition, when combining upwind schemes for compressible flow analysis with LSM-based gradient operators for spatial discretization, issues related to robustness and stability are observed in the local point clouds of specific geometric configurations. These problems can lead to inaccurate prediction of physical quantities or even numerical divergence in localized regions, which ultimately propagate and deteriorate the accuracy of the global numerical solution. Finally, in terms of efficiency, high computational cost of matrix inversion, which is an essential procedure of LSM-based gradient operators for spatial discretization, remains an unresolved challenge. Therefore, the successful implementation of meshless framework in CFD fundamentally requires an accurate LSM-based gradient operator, numerical strategies for robustness and stability, and computationally efficient algorithms for fast inversion.

First, mathematical correlation between numerical accuracy of LSM-based gradient operator and geometric characteristics of two-dimensional local point clouds, including AR, curvature, and skewness, is theoretically analyzed. Based on this theoretical analysis, two gradient correction approaches are proposed to reduce numerical errors in gradient estimation across multi-dimensional spaces: (1) a multi-stage optimization procedure using the new LSM based on the principle of superposition, and (2) a Lagrange Multiplier (LM)-based LSM for improved gradient estimation in skewed point clouds. To ensure robust operation of the newly developed method-based algorithms, limiting strategies are introduced from the perspectives of boundedness and condition number. The new
gradient operator is evaluated using various test functions, local point clouds, and varying levels of point connectivity in accuracy studies. This operator is evaluated by comparison with conventional approaches, including LSM, Geometric Conservation Least Squares Method (GC-LSM), Green-Gauss theorem (GG), and Green-Gauss theorem with centroidal values using volume-weighted averaging (SGG). In addition, numerical simulations are conducted to investigate the influence of gradient operators on spatial discretization and gradient reconstruction in high-order flux schemes. Furthermore, the time convergence characteristics of the proposed operator are also examined.

Second, the robustness and stability issues arising from the combination LSM-based spatial discretization in meshless method and upwind scheme under specific local point cloud configurations are analyzed in one-dimensional space to gain numerical and physical insights. In addition, because numerical results from the two-dimensional blast wave problem reveal that the findings derived from one-dimensional analysis are insufficient to fully capture the complex behavior observed in multidimensional problems, the one-dimensional analysis was extended to multi-dimensional space, leading to the derivation of new mathematical conditions: one to ensure robustness, and another to guarantee conditional stability based on von Neumann stability analysis. To satisfy these conditions, a new algorithm is developed in this work to generate aligned meshless coefficients. Furthermore, in order to mitigate the potential loss of numerical accuracy resulting from enforcing robustness and stability, an additional indicator is introduced. These proposed methods are validated on point clouds including sharp features, demonstrating their effectiveness in terms of robustness and stability.

Third, a fast matrix inversion algorithm is developed to alleviate the high computational cost associated with local approximation in LSM-based gradient operators with LM for multi-dimensional space. To reduce the time complexity of matrix inversion in the LM-based LSM, tensor product techniques are employed to decompose high-dimensional block matrices into lower-dimensional sub-blocks and repetitive patterns. Singular Value Decomposition (SVD) is integrated into the algorithm to further improve efficiency in underdetermined systems, and block matrix inversion is adopted to enable the complete inversion process using only low-scale computations. The performance of the proposed algorithm is quantitatively evaluated against the conventional LU decomposition benchmark, demonstrating a significant improvement in computational efficiency.

All proposed methods are integrated into an in-house CFD solver implemented in Fortran 90. To evaluate the meshless framework’s applicability across a broad range of validation benchmarks, the following validation cases are considered: (1) inviscid flow (e.g., Sod shock tube), (2) viscous flow (e.g., shock wave boundary layer interaction), (3) turbulent flows (e.g., RAE2822 transonic airfoil, ONERA M6 wing), (4) hypersonic flows (e.g., frozen flow over a circular cylinder, non-equilibrium flow over a Hypersonic Glide Vehicle (HGV)), and (5) supersonic turbulent combustion flows (e.g., Burrows and Kurkov supersonic mixing/combustion). For performance benchmarking, FVM and GC-LSM are employed as reference spatial discretization approaches. The numerical results demonstrate that the developed meshless framework achieves significant improvements in accuracy, robustness, stability, and computational efficiency across a wide range of flow regimes.

국문 초록 (Abstract)

수년간, 미분 연산자의 핵심 구성 요소인 구배 연산자의 정확도 향상은 수학, 과학, 공학 전반에서 지속적인 핵심 연구 과제로 여겨져 왔으며, 역문제(inverse problem), 최적 설계(optimal design), 편미분방정식(Partial Differential Equation, PDE) 해석, 회귀 및 보간(regression & interpolation), 신호 및 영상 처리 (signal & image processing), 매니폴드 해석(manifold analysis), 인공지능(Artificial Intelligence, AI) 등 다양한 분야에 널리 적용되어 왔다. 특히, 구배 추정을 위한 대표적인 접근법인 최소제곱법(Least Squares Method, LSM)은 오랜 기간 동안 전산유체역학(Computational Fluid Dynamics, CFD) 분야에서 활용되어 왔다. 최근에는 gradient reconstruction 과정에서의 적용을 넘어, LSM 기반 공간 차분(LSM-based spatial discretization)이 무격자 기법에 적용되는 방식으로 그 활용 범위가 확장되고 있다.

그러나 정확도(accuracy) 측면에서는, 복잡한 형상과 높은 종횡비를 가지는 경계층 문제에서 물리량의 구배를 정확하게 추정하기 어렵다는 한계로 인해, 여전히 공간 차분 및 gradient reconstruction 과정에서 LSM을 활용하기에는 부족함이 있다. 또한, 압축성 유동 해석을 위한 상류차분법과 공간 차분을 위한 LSM 기반 구배 연산자를 결합할 경우, 특정 기하학적 형상을 갖는 국소 질점 분포계에서 강건성 및 안정성(robustness & stability) 문제점들이 관찰된다. 이러한 문제는 물리량의 부정확한 예측 및 국소적 발산을 초래할 수 있으며, 이는 전역 수치 해석 결과의 정확도를 저해하는 원인이 된다. 마지막으로, 공간 차분을 위한 LSM 기반 구배 연산자에서 요구되는 행렬 역산 과정은 높은 계산 비용을 수반하며, 이는 효율성 측면에서 해결되지 않은 과제로 남아 있다. 따라서, CFD 분야에서 무격자 프레임워크를 성공적으로 구현하기 위해서는 정확한 LSM 기반 구배 연산자, 강건성과 안정성을 위한 수치 전략, 고속 역산을 위한 효율적인 알고리즘이 요구된다.

본 연구에서는 첫 번째로, 2차원 국소 point cloud의 기하학적 특성(종횡비, 곡률, 왜도)과 LSM 기반 구배 연산자의 정확도 간의 수학적 상관관계를 이론적으로 분석하였다. 이론적 분석을 바탕으로, 다차원 공간에서 구배 예측 오차를 줄이기 위한 두 가지 구배 보정 기법들을 제안하였다: (1) 다단계 최적화 기법을 적용한 중첩 원리 기반의 새로운 LSM, (2) 왜도가 있는 질점 클라우드에서 구배 정확도를 향상시키기 위한 라그랑주 승수법(Lagrange Multiplier, LM) 기반의 LSM. 새롭게 개발된 알고리즘의 강건한 작동을 보장하기 위해, 유계성(boundness)과 조건수(conditioned number) 관점에서 제한 전략(limiting strategy)을 도입하였다. 제안된 구배 연산자는 다양한 시험 함수, 국소 질점 클라우드, 질점 커넥션 수에 대해 구배 예측에 대한 정확성 연구가 이루어졌으며, 기존 방법인 LSM, 기하학적 보존 최소제곱법(Geometric Conservation Least Squares Method, GC-LSM), 그린 가우스 정리(Green-Gauss theorem, GG), 체적 가중 중심값을 사용하는 그
린 가우스 정리(SGG) 방식들과 비교되었다. 또한, 구배 연산자가 공간 차분 및 gradient reconstruction 기반의 고차 플럭스 스킴에 끼치는 영향을 확인하기 위해 수치 해석을 진행하였으며, 추가로 제안된 연산자의 시간 수렴성에 대한 논의도 이루어졌다.

두 번째로, 특정 기하학적 형상을 갖는 국소 질점 클라우드에서 무격자 기법에서의 LSM 기반 공간 차분과 상류 차분법을 결합할 때 발생하는 강건성 및 안정성 문제를 1차원 공간으로 단순화하여 수치적/물리적 통찰을 얻기 위한 분석을 수행하였다. 또한, 2차원 폭풍파(Blast wave) 수치 해석 결과에서 1차원 해석을 통해 도출된 결과만으로는 다차원 문제에서 나타나는 복잡한 현상을 충분히 설명하기 어렵다는 사실을 감안하여 1차원 해석을 다차원 공간으로 확장하였으며, 이를 통해 강건성을 만족하기 위한 조건과 폰 노이만 안정성 해석(von Neumann stability analysis)을 기반으로 조건부 안정성(conditionally stable condition)을 보장하는 새로운 수학적 조건을 유도하였다. 이러한 조건을 만족시키기 위해 정렬된 무격자 계수를 생성하는 새로운 알고리즘을 개발하였으며, 강건성과 안정성을 확보하는 과정에서 발생할 수 있는 수치 정확도 손실을 최소화하기 위해 추가적인 indicator를 도입하였다. 제안된 기법은 첨점을 포함하는 질점 클라우드들을 대상으로 검증되었으며, 강건성와 안정성 측면에서 제안된 방식의 효과를 입증하였다.

세 번째로, 다차원 LM을 활용한 LSM 기반 구배 연산자의 국소 근사화 과정에서 요구되는 높은 계산 비용을 절감하기 위해 고속 행렬 역산 알고리즘(fast inversion algorithm)을 개발하였다. LM 기반 LSM에서의 행렬 역산 시간 복잡도를 줄이기 위해, 텐서 곱(tensor product)을 활용하여 고차원 행렬의 블록 행렬들을 저차원 형태의 블록 행렬과 반복되는 패턴의 텐서 곱으로 재구성하였다. 또한, 기저 함수의 수 증가에 따른 시간 복잡도 문제를 해소하기 위하여 특이값 분해 (Singular Value Decomposition, SVD)를 알고리즘에 도입하였다. 마지막으로, 블록 행렬 역산(block matrix inversion) 기법을 적용하여, 전체 행렬 역산을 저차원 연산만으로 처리할 수 있도록 구성하였다. 제안된 알고리즘의 계산 성능을 정량적으로 평가하기 위해 수치 실험을 수행하였으며, 그 결과 LU 분해(Lower–Upper decomposition, LU decomposition) 기법을 기준으로 하여 본 알고리즘이 계산 효율성을 획기적으로 향상시킴을 확인하였다.

상기 모든 기법을 Fortran 90 기반의 in-house CFD 솔버와 통합하였다. 무격자 프레임워크를 평가하기 위해 전산유체역학 benchmark 문제들에 대해 수치 해석을 수행하였으며, 다음과 같은 검증 문제들을 선정하였다: (1) 비점성 유동 (예: Sod 충격파 관 문제), (2) 점성 유동 (예: 충격파–경계층 상호작용), (3) 난류 유동 (예: RAE2822 천음속 에어포일, ONERA M6 날개), (4) 극초음속 유동 (예: 원형 실린더 주변의 동결 유동, 극초음속 활공 비행체 주변의 비평형 유동), (5) 초음속 난류 연소 유동 (예: Burrows and Kurkov 초음속 혼합/연소). 본 프레임워크의 수치적 성능을 평가하기 위하여, FVM 및 GC-LSM 공간 차분 기법을 비교군으로 설정하였으며, 다양한 유동 문제에 대한 수치 실험 결과를 통해 본 연구에서 제안한 무격자 프레임워크가 정확도, 강건성, 안정성, 계산 효율성 측면에서 유의미한 개선을 달성했음을 확인하였다.

목차 (Table of Contents)

• Chapter 1 Introduction 1
• 1.1 Significance of Least Squares Method (LSM)-Based Gradient Operators 1
• 1.2 Knowledge of Meshless Methods 5
• 1.2.1 Key Advantages of Meshless Methods 5
• 1.2.2 Historical Development of Meshless Methods 7
• 1.2.3 LSM-Based Spatial Discretization in Meshless Method for Compressible Flow: Application of LSM-Based Gradient Operator 10
• 1.3 Challenges in Applying LSM-Based Gradient Operators to Meshless Methods 12
• 1.3.1 Accuracy 12
• 1.3.2 Robustness and Stability 14
• 1.3.3 Efficiency 15
• 1.4 Outline of the Thesis 15
• Chapter 2 Overview of Gradient Operators 20
• 2.1 LSM Operator 23
• 2.1.1 LSM 23
• 2.2 LM-based LSM Operator 25
• 2.2.1 First-order Consistency 25
• 2.2.2 Geometric Conservation Law 26
• 2.2.3 GC-LSM 26
• 2.3 GG Operator 29
• 2.3.1 GG 29
• 2.3.2 SGG 30
• Chapter 3 Accurate Gradient Operator 31
• 3.1 Accuracy Limitations of LSM-Based Gradient Operator with IDW, 31
• 3.2 Principle of Superposition-Based LSM 45
• 3.2.1 Idea and Conceptualization 45
• 3.2.2 New Objective Function and Optimization Procedure 49
• 3.2.3 Generation of Basis Set for Weighting Coefficients 54
• 3.2.4 Limiting Strategy 60
• 3.2.5 General Framework of the Principle of Superposition-Based LSM 62
• 3.2.6 Algorithms 64
• 3.3 Accuracy Studies 67
• 3.3.1 Two-Dimensional Space 67
• 3.3.2 Three-Dimensional space 103
• 3.4 Numerical Experiments 108
• 3.4.1 Influence of Gradient Operator on Spatial Discretization 113
• 3.4.2 Influence of Gradient Operator on Gradient Reconstruction in High-Order Flux Schemes 118
• 3.4.3 Comparison of Time Convergence 120
• Chapter 4 Robust and Stable LSM-Based Spatial Discretization 123
• 4.1 Robustness and Stability Issues 123
• 4.2 Cure for Robustness and Stability Issues in Multi-Dimensional Space 129
• 4.2.1 General Conditions for Robustness in Multi-Dimensional Space 129
• 4.2.2 Motivation for Multi-Dimensional Stability Analysis 129
• 4.2.3 Multi-dimensional von Neumann Stability Analysis 140
• 4.2.4 Numerical Strategies for Curing Robustness and Stability Issues 143
• 4.3 Numerical Experiments 150
• Chapter 5 Fast Inversion Algorithm for LM-Based LSM 157
• 5.1 Excessive Computational Cost in LM-Based LSM 157
• 5.2 Fast Inversion Algorithm 163
• 5.2.1 Tensor Product-Based Matrix Decomposition 163
• 5.2.2 Incorporation of SVD 164
• 5.2.3 Application of Block Matrix Inversion 164
• 5.2.4 Conceptual Overview of the Proposed Algorithms 171
• 5.3 Numerical Experiments 175
• Chapter 6 Physical Models and Numerical Strategies 182
• 6.1 Physical Models 182
• 6.1.1 Governing Equations 182
• 6.1.2 Thermochemical Properties 186
• 6.1.3 Transport properties 191
• 6.1.4 Chemical Reactions 196
• 6.1.5 Turbulent Modeling 204
• 6.2 Numerical strategies 207
• 6.2.1 Spatial discretization 207
• 6.2.2 High-Order flux schemes 207
• 6.2.3 Time integration 213
• 6.2.4 Boundary condition 215
• Chapter 7 Meshless Framework-Based In-house Codes 217
• 7.1 Preprocessing Procedure 217
• 7.1.1 Point Generation Technique 218
• 7.1.2 Parallelization, 219
• 7.2 Detailed Description of Meshless Framework-based in-house CFD Solver 221
• 7.2.1 Code structure 221
• 7.2.2 Outsource 224
• Chapter 8 Numerical Experiments 225
• 8.1 Inviscid Flow 225
• 8.1.1 Sod shock tube 225
• 8.2 Viscous Flow 244
• 8.2.1 Shock Wave Boundary Layer Interaction 244
• 8.3 Turbulent Flow 251
• 8.3.1 RAE2822 Transonic Airfoil 251
• 8.3.2 ONERA M6 Wing 258
• 8.4 Hypersonic Flow 267
• 8.4.1 Frozen Flow 267
• 8.4.2 HGV with Non-equilibrium Flow 275
• 8.5 Combustion Flow 279
• 8.5.1 Burrows and Kurkov Supersonic Mixing/Combustion 279
• Chapter 9 Mathematical Proofs 286
• Chapter 10 Concluding Remarks 301
• 10.1 Summary 301
• 10.2 Originality and Contribution 307
• 10.3 Future works 310
• Reference 312
• 국문초록 335