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

Ships, Splashes, and Waves on a Vast Ocean 논문을 이해하는 것을 목표로 Robert Bridson 의 Fluid Simulation for Computer Graphics 등을 기반으로 정리한 내용입니다.

1. Navier-Stokes 방정식

유체의 움직임을 제어하는 기본방정식은 나비에-스토크스 방정식이다.

위 식에서 u가 유체의 속도장을 의미한다. 참고로 마지막 항은 점성에 대한 것인데, 점성이 없는 물에 대해 다룰 것이므로 해당 항은 고려하지 않는다.

위 식이 나비에-스토크스 방정식이고, 아래 식은 유체가 비압축성이 되도록 하는 식이다. 즉 두 식을 이용하여 비압축성 유체의 움직임을 표현할 수 있다.

• 비압축성 유체란 ? 압축을 가하였을 때 부피가 변하지 않는 유체, 다른 물리량(압력, 온도)등의 변화가 밀도에 영향을 미치지 않는 유체를 말한다. 대부분의 액체는 비압축성이다.
  물이 꽉 차있는 페트병에 물을 압축시켜 더 채워넣을 수는 없다. 이게 비압축성을 의미한다. 그리고 기체는 압축성이다.

2. 오일러(Eulerian) 관점과 라그랑지안 (Lagrangian) 관점

연속체의 움직임에 대해 생각할 때, 운동을 추적하는 데는 두 가지 접근법, 라그랑지안 관점과 오일러 관점이 있다.

라그랑지안은 입자 하나하나를 따라가면서 이 입자가 어떻게 변하는지를 살펴보는 것이고, 오일러 관점에서는 입자를 추적하지 않고, 어느 고정된 '위치'를 기준으로 그 위치의 물리량이 어떻게 변하는지를 살펴본다.

예를 들어 온도 10도인 입자가 지나가고 온도 20도인 입자가 지나간다고 했을 때, 오일러 관점에서는 그 온도 10인 입자가 온도가 변하는지 안변하는지에는 관심이 없고 단지 여기 우리가 보고있는 위치의 온도가 10도에서 20도로 변했다 이런 관점으로 말하는 것이다.

라그랑지안 관점은 입자 시스템에, 오일러 관점은 고정된 그리드 시스템에 해당한다고 하기도 한다.

유체역학에서는 보통 오일러관점을 이용한다.

그런데 오일러 관점에서는 입자를 추적하지 않으니까 주어진 입자의 물리량이 얼마나 빠르게 변하는지, 예를 들어 가속도 같은 것들을 판단할 수 가 없다.

그래서 오일러 관점에서 물질을 추적하면서 미분하기 위해 등장한게 물질 미분(material derivate)이다.

어떤 물리량에 대한 수식(ex. q(t,x) - 시간 t, 위치 x에서 물리량 q의 값)를 그냥 통으로 미분해버리면 물질 미분식이 나온다.

그리고 이 물질 미분을 이용하여 아까 말한 나비에-스토크스 방정식을 표현하면, 오일러 관점에서의 비압축성 나비에-스토크스 방정식이 된다. 이게 논문에서의 주된 방정식이다.

3. Numerical Simulation - Split

이제 저 방정식으로 표현되는 유체를 컴퓨터로 시뮬레이션 하려면 이산화 과정이 필요하다. 이산화는 연속적인 수치를 이산적인, 불연속적인 수치로 변화시키는 것이다.

그러기 위해 먼저 복잡한 나비에-스토크스 방정식을 분할한다.

방정식을 3개로 나눈 다음, 어떤 항을 계산하고, 그 효과를 포함시키면서 갱신하고 이런식으로 진행된다.

나비에-스토크스 방정식을 아래처럼 3개로 쪼갠다.

1. 이류(advection) : 이류란 물리량(농도, 온도 등)의 차이로 인해 유체가 흐르는 것을 말한다. 우리는 속도에 관심이 있으므로 물리량 q를 아래처럼 u로 놓았다.
이류 항에서는 아래 간단한 이류 방정식을 풀면 된다.

논문에서는 초기 속도를 u-, 이류 후 속도를 u, 압력 투영 후 속도를 u+라고 보았다.

즉, 위 식을 풀어 우리는 u-에서 u로 값을 업데이트 한다.

그리고 여기다가 중력을 고려하면 된다. 그게 두번째 항이다.

2. body force : 부피 전체에 걸쳐 작용하는 힘이다. 여기서는 중력만을 의미하는 것으로 보인다.

3. 압력 투영 및 비압축성 (pressure projection + incompressibility) : 우리는 비압축성 유체를 다루기 때문에 비압축성 제약을 강제해야 한다. 그리고 그걸 하게 되는 항이 이 3번째 항이다.
압력이 부피가 축소되지 않도록 하는 일을 한다고 보기 때문에 여기서 이걸 고려한다.

여기서 일단 앞 식, 속도장을 시간에 대해 편미분 한 식을 이산화해보자.

시간에 대한 미분식을 이렇게 이산화 하면 3번째 항은 아래와 같이 변한다.

논문에서는 아래와 같이 pressure projection 후 속도 u+에 대한 식으로 표현하였다.

여기서 우리가 모르는 것 즉 구해야되는 것은 압력 P이다. 논문에서는 이를 구하기 위해 포아송 방정식(Poisson's equation)을 이용한다.

이 식을 푸는 방법이 BEM과 FLIP에서 다르다. BEM은 아직 공부하지 않았으니 우선 FLIP에서 이를 푸는 방법만 알아보자.

4. How to get Pressure

우선 기본적으로 포아송 방정식을 풀어 압력을 구하기 위해 유한차분(finite difference)을 이용한다.

유한차분이란 함숫값의 차를 이용해 미분계수를 근사하는 방법이다. 이렇게 말하면 잘 안와닿을 수 있는데,

아래와 같은 함수 f(x)가 있다고 가정해보자.

그럼 이 함수에서 도함수 f'(x)는 아래와 같이 근사할 수 있다.

이게 유한차분이다.

이걸 하는 이유는 포아송 방정식을 이산화하기 위해서이다.

그리고 2D 이상에서 이 유한차분을 실행하는 방법이 '그리드(Grid)'이다. 공간을 그리드로 차분화(Discretization)하여, 아래 그림 처럼 x의 성분, y의 성분, 압력 등을 다 다른 위치에 저장한다.

압력을 그리드의 중앙에, x성분을 세로선의 중앙에 u로, y성분을 가로선의 중앙에 v로 저장한다.

그럼 이걸 이용하여 아까 포아송 방정식의 첫번째 식의 우변(상수 제외)을 이산화해보자.

첫번째 식은 그냥 공식이다.

그리고 이걸 그리드 격자를 기반으로 이산화 한 것이 두번째 줄의 식이다.

그리고 이 식은 행렬의 곱으로 표현할 수 있다.

이 행렬곱을 보면 좌변이 div에 대한 식, 우변이 속도에 대한 식이라는 걸 알 수 있다.
즉 div를 이산화한게 저 좌변의 식이라는 거 같은데, 이건 좀 더 알아봐야 될 거 같다.

이렇게 해서 우변을 이산화 했다!

그럼 이제 같은 방법으로 좌변에서 divP를 이산화해보자.

마찬가지로 행렬 곱으로 나타내면 아래와 같다.

이렇게 하여 압력과 관련된 포아송 방정식을 모두 이산화 하였고, 행렬 곱으로까지 나타내었다 !

정리해보면 아래와 같다.

이렇게 해서 나비에-스토크스의 세번째 식, pressure projection이 얼추 해결되어 보인다.

그럼 이제 첫번째 식 advection을 이산화하는 방법에 대해 알아보자.

5. Semi-Lagrangian Advection

advection을 위에서 사용한 방법처럼 이산화하기는 어렵다.

대신 여기에서는 semi-Lagrangian method 라고 불리는 방법을 이용한다.

공간의 어떤 점 x에서 q의 새로운 값을 얻기 위해서는 개념적으로 x에 도달하는 입자를 찾아 q의 값을 구할 수 있다.

이 논리를 그리드에 적용한게 semi-Lagrangian method이다.

예시를 들어 자세히 살펴보면 다음과 같다.

그렇기 때문에 현재 어떤 그리드 포인트에서의 물리량을 알고 싶으면 속도장을 따라 역추적하면 된다.

그렇다면 Xp의 위치는 어떻게 알아낼까?

입자는 아래 방정식에 따라 ∆t의 시간동안 Xg까지 움직인다.

따라서 Xg의 위치로 Xp의 위치를 구하면 아래와 같다.

단순하게, u(Xg)는 속도에 관한 식이고, ∆t는 시간이니까 시간*속도해서 ∆t 시간동안 이동한 거리만큼을 빼준 것이다.

이런식으로 어떤 입자가 시작한 위치인 Xp를 구하고 나면, 이 Xp에서의 물리량 q가 현재 Xg에서의 물리량이 된다.

---

이렇게 해서 유체 시뮬레이션의 기본이 되는 나비에-스토크스 방정식을 3개로 분할하고, 각각을 이산화하는 방법에 대해 알아보았다.

다음 글에서는 FLIP에 대해 더 알아보도록 하자! 🙂