📘 유체 시뮬레이션(2) - FLIP

이전 포스트에서 나비에-스토크스 방정식을 이산화 하는 방법에 대해 알아보았다.

📘 유체 시뮬레이션(1) - Navier-Stokes 방정식

압력 투영(pressure projection)은 그리드를 이용하여 문제를 해결하지만, 이류(advection)는 압력을 이용한 해결이 어려워 그리드 내 입자를 이용하여 풀어야 했다.

그렇기 때문에 파티클에 대한 물리량을 그리드로 바꾸고 저장하는 과정과, 반대로 그리드의 물리량을 파티클에 저장하는 과정이 필요하다.

이 과정은 크게 PIC(Particle-In-Cell)과, FLIP(Fluid-Implicit-Particle) 두가지 방법으로 나뉜다.

이 중 특히 FLIP에 대하여 자세히 알아보고자 한다.

Ships, Splashes, and Waves on a Vast Ocean 논문에서는 전반적인 FLIP과정에서 Animating Sand as a Fluid 논문에 제시된 방법을 따랐다고 되어 있다.

따라서 Animating Sand as a Fluid 에서 사용한 방법을 기반으로 알아보았다.

FLIP pseudocode

1) 입자 값 q_p를 그리드 값 q_i,j,k로 바꾼다.
2) 그리드 값 q_i,j,k를 그리드에 저장한다.
3) 압력 투영, 외부 힘 등을 계산하여 그 결과인 q(new)_i,j,k를 얻는다.
4) 각각의 입자에 대해 그리드에서 변경된 값 (delta)q_i,j,k = q(new)_i,j,k - q_i,j,k를 보간하여 입자의 값에 추가한다.
5) 그리드 속도장 내 입자들에 이류를 적용한다.

(q: 물리량. 속도에 대해 다루고 있으므로 아래에서는 q를 u로 바꿔서 서술하였다)

각 단계들을 좀 더 자세히 살펴보자.

1,2) Particle to Grid - 입자의 값 q_p를 그리드 값 q_i,j,k로 바꾸고 저장

논문에서는 2x2x2 크기의 그리드 하나에 8개의 입자를 담았다.

그런다음 각 그리드 포인트는 근처 입자들의 가중 평균 값을 가지게 된다.
여기서 '근처'의 기준은 그리드 점을 중심으로 한 그리드 셀 폭의 2배인 큐브에 포함되는 것으로 정의하였다.

이걸 수식으로 나타내면 아래와 같을 것이다.

(점 A에 저장되는 속도의 값)

가중치는 standard trilinear weighting을 이용하였다고 되어 있다.

3) 압력 투영을 계산하여 그 결과인 q(new)_i,j,k를 얻는다.

먼저 중력에 대한 힘을 계산해야 한다.
이건 간단하다. y성분의 속도를 g*timestep 만큼 증가시키면 된다.

그런 다음 압력 투영을 적용해야 한다. 이전 글에서 압력을 구하는 방정식인 포아송 방정식을 이산화 하고 이를 행렬 곱으로 표현하였다.

이렇게 표현된 식을 더 단순화 하여 Ap = d의 형태로 표현할 수 있다. 따라서 Ap = d로 부터 압력을 얻어야 하며 이 때 conjugate gradient method를 이용한다.

그리고 여기서 구해진 압력을 나비에-스토크스 방정식의 pressure projection과 관련된 수식에 넣으면 업데이트 된 물리량(여기선 속도)를 얻을 수 있다.

이렇게 해서 중력과 압력 투영이 적용된 새로운 물리량을 얻었다. PIC에서는 이 값을 그대로 기존 값과 바꾸지만 FLIP에서는 기존 값과 변경된 값의 차이를 이용한다.

4) 각각의 입자에 대해 그리드에서 변경된 값 (delta)q_i,j,k = q(new)_i,j,k - q_i,j,k를 보간하여 입자의 값에 추가한다.

이렇게 새로운 속도에서 기존 속도를 뺀 값을 구한 다음, 그 값을 보간하여 기존 입자의 값에 추가해준다.

이를 수식으로 나타내면 아래와 같다.

위 수식에서 2D의 상황이라 가정하고 A,B,C,D 4개만 적었는데 실제로 우리는 3D를 다루니 4개가 아닌 8개의 그리드 지점에서 기록된 속도 변화를 삼선형 보간(trilinear interpolate)해야 된다.

5) 그리드 속도장 내 입자들에 이류를 적용한다.

이제 중력과 압력 투영이 적용된 값을 입자에 업데이트 했으니, 입자에 이류를 적용하면 된다.

해당 논문에서는 CFL 조건에 의해 제한된 5개의 서브스텝들이 있는 RK2 ODE 솔버를 이용하였다고 되어 있다.

CFL은 적절한 time step을 구하는 것에 대한 조건으로 이를 적용하여 각 서브스텝에서 입자가 1개의 그리드 셀 이상을 움직이지 못하도록 한다.

RK2 ODE 솔버는 2차 Runge-Kutta 방법이다.

이전 포스트에서 이류를 해결하기 위해 semi-lagrangian 방법을 이용해서 입자의 현재 위치 Xg에서 과거 위치 Xp를 구할 때 아래와 같은 간단한 방법을 이용했었다.

여기서 좀 더 발전한 방법이 2차 Runge-Kutta 방법이다.

이런식으로 하나의 time step에서 두번에 걸쳐 Xp 값을 구하는 것이다. 1/2 timestep에서의 중간 위치를 구하고, 그 위치를 이용하여 최종적으로 Xp를 얻는 방법이다.

이제 정말로 Xp를 구했다. 입자가 시작된 위치 Xp를 알았으니 그 위치에서의 물리량만 알면 이류 문제가 해결된다.

그런데 과거의 위치 Xp가 현재 표현되고 있는 그리드에 없을 가능성이 있다. 그럴 경우 인근 그리드 포인트에서 보간하여 적절한 근사치를 얻을 수 있다.

여기서도 마찬가지로 삼선형 보간을 이용할 수 있다.

이렇게 하여 전반적인 FLIP이 적용되는 과정에 대해 알아보았다.