[NDC2018] 유체역학 엔진이 직면한 문제와 미래

▲ 넥슨코리아 옥찬호

유체역학은 게임 프로그래머들에게 다소 낯선 주제다. 유체역학은 유체를 다루는 학문인데, 게임에서 유체역학 엔진을 사용하면 물이나 눈 등 유체를 좀 더 현실감 있게 표현할 수 있다. 문제는 유체역학을 처리하는 계산량이 너무 많다는 점이다.

하지만 최근에는 딥러닝을 융합해 계산 속도를 비약적으로 개선하려는 연구가 진행 중이며 게임에서 실제로 유체 역학을 사용한 사례도 한둘씩 생기고 있다. 넥슨코리아의 옥찬호는 오픈소스 유체역학 엔진을 만들면서 어떤 기술을 적용했는지 경험을 공유했다.

---

유체역학?

---

옥찬호는 겨울왕국이 인기를 끌 무렵 겨울왕국 속 눈 움직임에 관심을 가졌다. 실사와 같이 표현한 눈, 파라미터에 따라 눈이 구체화되는 모습에 관심을 가졌다. 심심풀이로 시작하려 했지만, 쉽지 않았다. 관련 서적을 보았지만, 하얀 건 종이고 검은 건 글씨였다. 그러다 어떤 서적을 계기로 다시 유체 역학에 관심을 가지기 시작했다. 목표는 게임에서 유체역학을 쓸 수 있게 보고자 함이었다.

사실 유체 역학을 활용해 만든 게임이나 애니메이션이 없는 건 아니다. 2011년 작 'Where's My Water?'나 'Sprinkle', 'Watee' 그리고 2012년도에 나온 'Splash!!' 같은 게임이 존재한다. 애니메이션에서는 종종 사용되고는 한다. 특히 디즈니는 애니메이션 한 편당 논문을 하나씩 낼 정도로 유체역학을 활용한다.

▲ 석사에게도 검은건 글씨고 흰건 종이였다.

유체역학은 유체의 운동에 관해서 연구하는 학문이다. 유체는 고체보다 형상이 일정하지 않아 변형이 쉽고 자유로이 흐를 수 있는 액체와 기체 그리고 플라즈마를 총칭하는 말이다. 유체는 흐르는 성질을 가지고 있어 모양이 정해져 있지 않다. 그래서 담는 모양에 따라 모양이 결정된다는 특징이 있다.

또한, 유체는 점성과 압축성이라는 특징도 가지는데, 점성은 유체의 서로 붙어있는 부분이 떨어지지 않으려는 성질을 뜻하며, 압축성은 압력이나 온도가 변하면 밀도가 변하는 성질을 말한다.

옥찬호는 본격적으로 유체역학 엔진을 만지기 전에 간단한 유체 시뮬레이터 '1차원 물결 시뮬레이션'을 콘솔로 구현했다. 서로 부딪힌 후 반대 방향으로 퍼져나가는 걸 구현하기 위해서 상태를 정의하고, 움직임을 계산하고, 충돌을 처리하고, 시각화하는 과정을 거쳤다. 다음은 해당 코드다.

#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif

#include
#include
#include
#include
#include
#include
#include

const size_t BUFFER_SIZE = 80;
const char* GRAY_SCALE_TABLE = " .:-=+*#%@";
const size_t GRAY_SCALE_TABLE_SIZE = sizeof(GRAY_SCALE_TABLE) / sizeof(char);

void UpdateWave(const double timeInterval, double* x, double* speed)
{
(*x) += timeInterval * (*speed);

// Boundary reflection
if ((*x) > 1.0)
{
(*speed) *= -1.0;
(*x) = 1.0 + timeInterval * (*speed);
}
else if ((*x) < 0.0)
{
(*speed) *= -1.0;
(*x) = timeInterval * (*speed);
}
}

void AccumulateWaveToHeightField(const double x, const double waveLength, const double maxHeight, std::array* heightField)
{
const double quarterWaveLength = 0.25 * waveLength;
const int start = static_cast((x - quarterWaveLength) * BUFFER_SIZE);
const int end = static_cast((x + quarterWaveLength) * BUFFER_SIZE);

for (int i = start; i < end; ++i)
{
int iNew = i;
if (i < 0)
{
iNew = -i - 1;
}
else if (i >= static_cast(BUFFER_SIZE))
{
iNew = 2 * BUFFER_SIZE - i - 1;
}

const double distance = fabs((i + 0.5) / BUFFER_SIZE - x);
const double height = maxHeight * 0.5 * (cos(std::min(distance * M_PI / quarterWaveLength, M_PI)) + 1.0);
(*heightField)[iNew] += height;
}
}

void Draw(const std::array& heightField)
{
std::string buffer(BUFFER_SIZE, ' ');

// Convert height field to grayscale
for (size_t i = 0; i < BUFFER_SIZE; ++i)
{
const double height = heightField[i];
const size_t tableIndex = std::min(static_cast(floor(GRAY_SCALE_TABLE_SIZE * height)), GRAY_SCALE_TABLE_SIZE - 1);
buffer[i] = GRAY_SCALE_TABLE[tableIndex];
}

// Clear old prints
for (size_t i = 0; i < BUFFER_SIZE; ++i)
{
printf("b");
}

// Draw new buffer
printf("%s", buffer.c_str());
fflush(stdout);
}

int main()
{
const double waveLengthX = 0.8;
const double waveLengthY = 1.2;

const double maxHeightX = 0.5;
const double maxHeightY = 0.4;

double x = 0.0;
double y = 1.0;
double speedX = 1.0;
double speedY = -0.5;

const int fps = 100;
const double timeInterval = 1.0 / fps;

std::array heightField;

for (int i = 0; i < 1000; ++i)
{
// March through time
UpdateWave(timeInterval, &x, &speedX);
UpdateWave(timeInterval, &y, &speedY);

// Clear height field
for (double& height : heightField)
{
height = 0.0;
}

// Accumulate waves for each center point
AccumulateWaveToHeightField(x, waveLengthX, maxHeightX, &heightField);
AccumulateWaveToHeightField(y, waveLengthY, maxHeightY, &heightField);

// Draw height field
Draw(heightField);

// Wait
std::this_thread::sleep_for(std::chrono::milliseconds(1000 / fps));
}

printf("n");
fflush(stdout);

return 0;
}

---

유체역학의 핵심 나비에-스톡스 방정식

---

▲ 나비에-스톡스 방정식

옥찬호는 "어려워 보이지만, 쪼개서 보면 어렵지 않다"라고 설명했다. 나비에-스톡스 방정식은 점성을 가진 유체의 운동을 기술하는 비선형 편미분방정식으로 뉴턴의 운동 제2 법칙(F=ma)을 유체역학에서 사용하기 쉽게 바꾼 것이다.

첫 번째 항은 속도를 시간으로 편미분 해서 얻는 가속도다. 두 번째 항은 압축성을 가지느냐, 가지지 않느냐에 따라 변화하는데 보통 전산 유체역학 에서는 비압축성 유체를 고려한다. 압축성을 다루면 계산량이 너무 많아지기 때문이다. 그래서 계산할 필요가 없다. 세 번째 항은 외력을 뜻한다. 중력, 마찰력 등이 속한다. 네 번째 항은 단위 면적당 받는 힘, 즉 압력이다. 마지막은 점성에 관련한 항으로 점성에 따라 유체를 판단하는 값으로 사용한다.

유체를 계산할 때 크게 두 가지 방식을 사용한다 첫 번째는 입자 방식(Particle-based)이며, 두 번째는 격자방식(Grid-based)이다. 입자 방식은 라그랑주(Lagrangian) 관점이라고도 하며 입자 하나하나에 초점을 맞춰 각 입자를 따라가면서 물리량을 나타내는 기술법이다. 시간이 지난 후의 위치를 매핑 함수를 이용해 표기한다.

격자 방식은 오일러(Eulerian) 관점이라고도 한다. 관찰할 지점을 고정하여 점을 지나는 유체의 물리량을 표현한다. 유체가 포함된 영역의 각 점에 속성값을 저장한다. 격자방식도 격자를 표현하고, 형태 저장하는 방식에 따라 많은 방식이 존재한다.

옥찬호는 컴퓨터 게임을 위한 복셀 기반의 유체엔진 CubbyFlow를 사용했다. CubbyFlow는 MIT 라이선스라 무료로 사용할 수 있다. 다양한 수학, 기하, 애니메이션, 충돌, 장, 입자, 격자, 하이브리드, 유틸 클래스를 이 Jet Framework 기반 유체엔진은 애니메이션과 프레임을 기반으로 피직스 애니메이션을 상속하는 방식으로 솔버(Solver)를 제공한다.

▲ 예시

#include
#include
#include
#include

using namespace CubbyFlow;

class SimpleMassSpringAnimation : public PhysicsAnimation
{
public:
struct Edge
{
size_t first;
size_t second;
};

struct Constraint
{
size_t pointIndex;
Vector3D fixedPosition;
Vector3D fixedVelocity;
};

std::vector positions;
std::vector velocities;
std::vector forces;
std::vector edges;

double mass = 1.0;
Vector3D gravity = Vector3D(0.0, -9.8, 0.0);
double stiffness = 500.0;
double restLength = 1.0;
double dampingCoefficient = 1.0;
double dragCoefficient = 0.1;

double floorPositionY = -7.0;
double restitutionCoefficient = 0.3;

VectorField3Ptr wind;

std::vector constraints;

SimpleMassSpringAnimation() {}

void MakeChain(size_t numberOfPoints)
{
if (numberOfPoints == 0)
{
return;
}

size_t numberOfEdges = numberOfPoints - 1;

positions.resize(numberOfPoints);
velocities.resize(numberOfPoints);
forces.resize(numberOfPoints);
edges.resize(numberOfEdges);

for (size_t i = 0; i < numberOfPoints; ++i)
{
positions[i].x = -static_cast(i);
}

for (size_t i = 0; i < numberOfEdges; ++i)
{
edges[i] = Edge{ i, i + 1 };
}
}

void ExportStates(Array1& x, Array1& y) const
{
x.Resize(positions.size());
y.Resize(positions.size());

for (size_t i = 0; i < positions.size(); ++i)
{
x[i] = positions[i].x;
y[i] = positions[i].y;
}
}

protected:
void OnAdvanceTimeStep(double timeIntervalInSeconds) override
{
size_t numberOfPoints = positions.size();
size_t numberOfEdges = edges.size();

// Compute forces
for (size_t i = 0; i < numberOfPoints; ++i)
{
// Gravity force
forces[i] = mass * gravity;

// Air drag force
Vector3D relativeVel = velocities[i];
if (wind != nullptr)
{
relativeVel -= wind->Sample(positions[i]);
}
forces[i] += -dragCoefficient * relativeVel;
}

for (size_t i = 0; i < numberOfEdges; ++i)
{
size_t pointIndex0 = edges[i].first;
size_t pointIndex1 = edges[i].second;

// Compute spring force
Vector3D pos0 = positions[pointIndex0];
Vector3D pos1 = positions[pointIndex1];
Vector3D r = pos0 - pos1;
double distance = r.Length();
if (distance > 0.0)
{
Vector3D force = -stiffness * (distance - restLength) * r.Normalized();
forces[pointIndex0] += force;
forces[pointIndex1] -= force;
}

// Add damping force
Vector3D vel0 = velocities[pointIndex0];
Vector3D vel1 = velocities[pointIndex1];
Vector3D relativeVel0 = vel0 - vel1;
Vector3D damping = -dampingCoefficient * relativeVel0;

forces[pointIndex0] += damping;
forces[pointIndex1] -= damping;
}

// Update states
for (size_t i = 0; i < numberOfPoints; ++i)
{
// Compute new states
Vector3D newAcceleration = forces[i] / mass;
Vector3D newVelocity = velocities[i] + timeIntervalInSeconds * newAcceleration;
Vector3D newPosition = positions[i] + timeIntervalInSeconds * newVelocity;

// Collision
if (newPosition.y < floorPositionY)
{
newPosition.y = floorPositionY;

if (newVelocity.y < 0.0)
{
newVelocity.y *= -restitutionCoefficient;
newPosition.y += timeIntervalInSeconds * newVelocity.y;
}
}

// Update states
velocities[i] = newVelocity;
positions[i] = newPosition;
}

// Apply constraints
for (size_t i = 0; i < constraints.size(); ++i)
{
size_t pointIndex = constraints[i].pointIndex;
positions[pointIndex] = constraints[i].fixedPosition;
velocities[pointIndex] = constraints[i].fixedVelocity;
}
}
};

예시에서 구현한 사항은 뉴턴 운동 제2 법칙, 중력, 훅의 법칙, 시간 통합 그리고 장애물이다. 유체 솔버 역시 예시와 마찬가지 방법으로 진행한다. 우선 고려해야 할 것은 물리적 상태를 어떻게 나타낼 것인가다, 이후 힘을 계산하고 움직임을 계산하고 제약 조건과 장애물을 적용하면 유체 솔버를 얻을 수 있다.

---

직면한 문제와 미래

---

지금까지 구현한 방식은 SPH, PCISPH solver를 비롯하여 Upwind/ENO/FMM Level set solver, Level set-based smoke solver, PIC/FLIP/APIC solver 등이 있으며 MPS/MPM solver, FAB solver, IVOCK fluid solver, Multiphase fluid solver 등은 아직 구현하지 못했다.

옥찬호가 유체역학 엔진을 만진 이유는 게임에 적용하기 위해서였다. 그러나 아직은 계산 속도 부분에서 문제를 안고 있다. 너무 느리기 때문이다. 비록 CPU 병렬을 통해 비약적인 성능 향상을 끌어냈지만, 게임에서 사용하기엔 아직 갈 길이 너무 멀다. 최근 머신러닝 기반 실시간 유체 시뮬레이션으로 계산 속도 문제를 해결하기 위한 연구가 진행되고 있다.

또한, 렌더링의 문제도 존재한다. 현재 CubbyFlow는 시뮬레이션만 수행하는 라이브러리이기 때문에 렌더러와 결합하는 작업이 필요하다. 그래서 DX11, DX12, OpenGL, Vulcan과 결합하려는 시도가 이어지고 있다. 파이선을 제외하고 지원하지 않는 다른 언어 API도 단점 중 하나다.

이처럼 유체역학 엔진을 만들기는 쉽지 않다. CPU/GPU 병렬프로그래밍을 통해 속도를 비약적으로 향상할 수 있었으나 아직 게임에서 실시간으로 유체역학을 사용하기는 무리다. 그러나 머신러닝을 활용해 사용 가능성을 연구하는 중이다.