Smoothed Particle Hydrodynamics(SPH) Fluid
카테고리: Graphics
태그: DirectX
홍정모님의 그래픽스 새싹코스 강의를 듣고 정리한 내용입니다.
Paper:
- SPH Fluids in Computer Graphics
Smoothed Particle Hydrodynamics
Smoothed Particle Hydrodynamics(SPH)의 기본적인 컨셉은 Incompressibility(비압축성) 이다. 고무 공의 한쪽 면을 누르면 반대쪽 면이 튀어나오는 것처럼 유체를 이루는 파티클들은 자신이 가지고 있는 물리량(ex. 밀도, 압력)을 균일하게 유지하려는 성질을 있다. 따라서 물리량이 높은 입자는 낮아지는 방향으로 속도를 생성하고, 물리량이 낮은 입자는 높아지는 방향으로 속도를 지정하는 방식이다.
먼저 임의의 position i에서의 물리량 A_i는 위와 같이 정의한다.
- 여기서 말하는 물리량에는 다양한 값들이 들어갈 수 있는데, SPH에서는 압력 p, 속도 v가 들어가게 된다.
- m_j : mass of particle j
- rho_j : density of particle j
- W_ij : Kernel Weight (주변 입자의 가중치)
- A_i의 Gradient, Laplacian은 수치적으로 위와 같이 표현
- 결국은 위 Navier-Stokes 방정식을 풀어야 함
- P term : 압력에 의한 변화량
- v term : 점성에 의한 변화량
- 점성(viscosity)은 주변 입자들의 속도 차이를 줄여주는 역할
- F term : 일반적으로 중력
- 식에 등장하는 Gradient, Laplacian은 위 A에 관한 식을 이용해 계산
- 압력 p의 경우 밀도 rho에 의해 결정되는 위 모델이 사용 (A 버전의 간소화 버전..?)
- rho_0는 유체가 평형 상태일 때의 밀도
- Kernel Function의 경우 위와 같은 Cubic-Spline 함수 사용
- Domain 범위가 0 ~ 2 이기 때문에 input은
distance/radius * 2.0f
SPH Algorithm
- radius 내에 포함된 주변 입자 찾기
- 주변 입자들을 이용해 density, pressure 업데이트
- Gradient, Laplacian 식을 이용해 Force 업데이트
- Navier-Stokes으로 구한 dv를 이용.
- v 업데이트
- 업데이트 된 v로 position 업데이트
Undate Density & Pressure
\[\rho_{i} = \sum_{j}^{} m_{j} W_{ij}\] \[p_{i} = k \left ( \left ( \frac{\rho_{i}}{\rho_{0}} \right )^{7} - 1 \right )\]void SphSimulation::UpdateDensity() {
#pragma omp parallel for
for (int i = 0; i < m_particlesCpu.size(); i++) {
m_particlesCpu[i].density = 0.0f;
for (size_t j = 0; j < m_particlesCpu.size(); j++) {
if (m_particlesCpu[j].life < 0.0f)
continue;
const float dist =
(m_particlesCpu[i].position - m_particlesCpu[j].position)
.Length();
if (dist >= m_radius)
continue;
m_particlesCpu[i].density +=
m_mass * SphKernels::CubicSpline(dist * 2.0f / m_radius);
}
m_particlesCpu[i].pressure =
m_pressureCoeff *
(pow(m_particlesCpu[i].density / m_density0, 7.0f) - 1.0f);
}
}
- 각 입자별로 주변 입자들을 찾아야 하기 때문에 이중 for문 형태
- 성능을 위해 OpenMP로 병렬 처리
if (dist >= m_radius)문으로 이웃 찾기- Density 먼저 계산 -> density를 이용해 Pressure 계산
Update Forces
\[\bigtriangledown A_{i} = \rho_{i} \sum_{j}^{} \left ( \frac{A_{i}}{\rho_{i}^{2}} + \frac{A_{j}}{\rho_{j}^{2}} \right ) \bigtriangledown W_{ij}\] \[\bigtriangledown ^{2} A_{i} = 2 \sum_{j}^{} \frac{m_{j}}{\rho_{j}} A_{ij} \frac{X_{ij} \cdot \bigtriangledown W_{ij}}{X_{ij} \cdot X_{ij} + 0.01 h^{2}}\]void SphSimulation::UpdateForces() {
#pragma omp parallel for
for (int i = 0; i < m_particlesCpu.size(); i++) {
// ...
for (size_t j = 0; j < m_particlesCpu.size(); j++) {
// ...
if (dist >= m_radius)
continue;
const Vector3 gradPressure =
rho_i * m_mass *
(p_i / (rho_i * rho_i) + p_j / (rho_j * rho_j)) *
SphKernels::CubicSplineGrad(dist * 2.0f / m_radius) *
(x_i - x_j) / dist;
const Vector3 laplacianVelocity =
2.0f * m_mass / rho_j * (v_i - v_j) /
(x_ij.LengthSquared() + 0.01f * m_radius * m_radius) *
SphKernels::CubicSplineGrad(dist * 2.0f / m_radius) *
x_ij.Dot(x_ij / dist);
pressureForce -= m_mass / rho_i * gradPressure;
viscosityForce += m_mass * m_viscosity * laplacianVelocity;
}
m_particlesCpu[i].force = pressureForce + viscosityForce;
}
}
- 식에 등장하는
(x_i - x_j) / dist은 방향을 나타내는 단위 벡터 - 코드의 재활용을 위해 중력 계산은 render() 부분에서..
Update Velocity & Position
// ...
for (int i = 0; i < m_particlesCpu.size(); i++) {
if (m_particlesCpu[i].life < 0.0f)
continue;
m_particlesCpu[i].velocity += m_particlesCpu[i].force * dt / m_mass;
m_particlesCpu[i].position += m_particlesCpu[i].velocity * dt;
}
Results
댓글 남기기