이전 글에서는 1차원 열방정식을 유한차분법(FDM)으로 근사 풀이하는 과정을 살펴보았습니다. 유한차분법은 격자점을 기반으로 단순하고 직관적인 접근법이지만, 복잡한 기하 형태나 경계조건을 다루기에 어려울 수 있습니다. 이번 글에서는 유한요소법(Finite Element Method, FEM)의 기초를 소개하고, 이를 통해 2차원 PDE 문제에 접근하는 방법에 대해 살펴보겠습니다.
유한요소법은 연속 영역을 삼각형, 사각형, 사면체, 육면체 등 요소(Element)로 분할하고, 각 요소 내에서 특정 형태 함수(Basis Function)를 사용해 PDE 해를 근사하는 강력하고 유연한 기법입니다. FEM은 기하학적 복잡성을 다루는 데 뛰어나며, 구조해석, 유체역학, 전자기학 등 광범위한 분야에서 사용됩니다.
FEM 기초 개념 정리
- 도메인 분할(Mesh Generation): PDE를 정의한 2D 영역(예: Ω)을 삼각형 요소들로 분할(mesh)합니다.
- 약형(Weak Form) 도출: PDE 문제를 적분 형태로 변환(weak form)한 뒤, 시험 함수(Test Function)와 기준 함수(Basis Function)를 사용해 이산화합니다.
- 어셈블(Assemble): 요소별로 국소적 stiffness matrix와 load vector를 계산하고, 이를 전역 행렬과 벡터로 어셈블합니다.
- 선형시스템 풀이: Ax=b 형태의 선형시스템을 해(A에 해당하는 것은 stiffness matrix, b는 load vector), 경계조건 적용 후 해를 구합니다.
- 해석 및 후처리: 계산된 해를 기반으로 물리량 해석, 시각화 수행.
이번 글에서는 전체 FEM 구현을 하기에는 매우 방대하므로, 핵심 개념에 초점을 맞추고, 단순화된 2D PDE 예제를 통해 FEM 접근 과정을 스케치해보겠습니다.
예제 PDE: 2D 푸아송 방정식(Poisson Equation)
푸아송 방정식:
-∇²u = f(x,y) in Ω
u = 0 on ∂Ω (Dirichlet boundary condition)
여기서 ∇²는 라플라시안 오퍼레이터, f(x,y)는 소스항.
FEM 접근:
- 도메인을 예: 단위사각형 Ω=[0,1]x[0,1]로 설정
- 이 도메인을 NxN 그리드 점을 바탕으로 삼각형 요소들로 분할 (각 사각형을 2개 삼각형으로 나눔)
- 각 요소 내에서 선형 기준 함수(Linear Basis Function) 사용
- 약형:
문제 -∇²u = f의 약형 형태는 ∫_Ω ∇u·∇v dx = ∫_Ω f v dx, (v는 시험함수)
경계조건 u=0 적용 시 ∂Ω에서 테두리 조건 반영 - 요소별로 stiffness matrix Ke와 load vector fe 계산 후 전역 시스템 어셈블
- 경계조건 처리 후 Ax=b 형태의 시스템 풀이
단순화된 FEM 구현 흐름 예시 (의사코드)
// src/fem.rs (개념적 스케치, 실제 동작 코드 아님)
use crate::matrix::Matrix;
use crate::vector::Vector;
pub struct Mesh {
pub coords: Vec<(f64,f64)>, // 노드 좌표
pub elements: Vec<[usize;3]>, // 삼각형 요소(노드 인덱스 3개)
}
pub fn generate_simple_mesh(N: usize) -> Mesh {
// 0..N 사이를 균등분할한 그리드 점을 만들고,
// 각 사각형 셀을 두 삼각형으로 나누어 mesh.elements 구성
// coords에 (x,y) 저장
// 간략화 위해 구현 생략
unimplemented!()
}
pub fn assemble_system(mesh: &Mesh, f: &dyn Fn(f64,f64)->f64) -> (Matrix<f64>, Vector<f64>) {
let nnodes = mesh.coords.len();
let mut A = Matrix::new(nnodes, nnodes, vec![0.0; nnodes*nnodes]);
let mut b = Vector::new(vec![0.0; nnodes]);
for elem in &mesh.elements {
// elem: [i,j,k] 삼각형 요소
// 노드 좌표 (x_i,y_i), (x_j,y_j), (x_k,y_k)
// 요소 stiffness matrix Ke, load vector fe 계산
// Ke_ij = ∫_elem ∇φ_i·∇φ_j dxdy
// fe_i = ∫_elem f(x,y)*φ_i dxdy
// 여기서 φ_i는 i번째 노드에 대응하는 선형 기준함수
// 실제 계산: Jacobian, 영역 적분, barycentric coords 이용
// 간단화: 구현 생략, 실제론 수치 적분(Gaussian Quadrature) 필요
// A와 b에 Ke, fe를 어셈블
unimplemented!()
}
// Dirichlet 경계조건 적용: 경계 노드 u=0
// A의 해당 행을 identity, b=0 세팅
unimplemented!();
(A, b)
}
// 이후 (A,b) 시스템을 LU 분해 등으로 해풀기
위 코드는 개념적 설명이며, 실제 FEM 구현은 훨씬 방대합니다.
- Mesh 생성기,
- 요소 적분(Gaussian Quadrature),
- Basis Function 정의,
- 경계조건 처리,
- 대규모 sparse 시스템 해석 등 고려 필요.
Sparse Matrix 고려
FEM 시스템은 대부분 희소행렬(Sparse Matrix)로 나타납니다. 대규모 PDE 문제에서는 희소성 활용이 필수적입니다. Rust에서는 sprs 등 희소행렬 크레이트나 C/C++ 라이브러리 연동을 통해 효율적인 선형 시스템 해법을 구현할 수 있습니다. 여기서는 개념만 다루고 구체적 구현은 생략합니다.
향후 개선점
- 고급 메셔(Mesh Generator): Gmsh나 mesh 러스트 라이브러리와 연계해 복잡한 영역 메쉬 생성
- 다양한 Basis Function: 선형, 2차 근사 등 고차 정확성 요소
- 비정형 메쉬: 삼각형/사면체 요소를 임의 기하에 맞춰 사용
- PDE 솔버 고급 기법: FEM-FVM 혼용, 적응적 메쉬 리파인먼트, 병렬 해석(Rayon, MPI)
- 비선형 PDE, 시간 의존 PDE: Newton-Krylov 방법, implicit time stepping
결론
이번 글에서는 FEM의 개념과 2D PDE 문제 접근법을 간략히 소개했습니다. FEM은 유한차분법보다 복잡하지만, 기하적 유연성과 수학적 정교함 덕분에 실무 엔지니어링 문제에서 널리 쓰입니다. 러스트로 FEM 기반 PDE 솔버를 구현하려면 메쉬, 적분, 경계조건 처리, 희소행렬 해법 등 다양한 영역을 결합해야 하지만, 언어의 안전성과 생태계 확장이 이러한 도전을 지원합니다.
유용한 링크와 리소스
- FEM 기초: Zienkiewicz & Taylor "The Finite Element Method"
- Gmsh: https://gmsh.info/ (메쉬 생성기)
- Rust Sparse Matrix: sprs crate
- Rust FEM 예제: rfenics(개발 중), ndarray-linalg 연계