종합 실험실 · 한 항씩 쌓는 유체역학

0장의 실험실로 돌아왔다. 처음 그 상자에서 물을 저었을 때는 실험 창 왼쪽 위의 칩 다섯 개가 무엇인지 몰랐다. 이제는 하나하나가 무슨 일을 하는지 안다. 0장은 마지막 장에 이르면 쌓아 온 항들이 그대로 하나의 유체 시뮬레이터가 된다고 했고, 이 장이 그 마지막 장이다. 먼저 아직 따로 다루지 않은 항, 바깥에서 주는 힘을 마저 본다. 그다음 실험실의 항을 하나씩 끄면서 항마다 하는 일을 다시 확인하고, 모든 항이 함께 일하는 흐름 하나, 날개를 지나는 흐름을 읽는다. 끝으로 실험실을 움직이는 엔진의 코드를 단계마다 펼쳐 본다. 방정식 막대에서는 모든 항이 강조되어 있다.

바깥에서 주는 힘

방정식 오른쪽 끝의 $\frc{\mathbf{f}}$ 는 바깥에서 물에 주는 힘을 단위 질량당 값으로 적은 것이다. 0장에서 물을 저은 손이 바로 이 항이라고 했다. 엔진은 손이 지나간 자리 둘레의 물에 손이 움직인 방향으로 속도를 더한다. 손에서 멀어질수록 빨리 작아지는 종 모양으로 더하므로, 손 바로 밑의 물이 가장 세게 밀린다. 바깥 힘에는 물 전체에 한꺼번에 작용하는 것도 있다. 중력이 그렇고, 데워진 물이 둘레의 찬물보다 가벼워 떠오르게 하는 부력도 바깥 힘이다. 이 장의 실험실에는 젓는 손의 힘만 있다.

물 전체를 고르게 미는 힘이 늘 흐름을 만드는 것은 아니다. 컵에 담긴 물은 중력을 받으면서도 가만히 있다. 아래로 갈수록 압력이 높아져서, 그 압력의 차이가 물을 위로 밀어 중력과 꼭 맞서기 때문이다. 3장의 압력은 물이 몰리거나 비지 않도록 어디서나 맞춰지는 장이었다. 아래로 몰리려는 물을 떠받치는 것도 그 압력이다. 엔진에서도 그렇다. 0장의 닫힌 상자에 물 전체를 아래로 미는 힘을 주면, 압력이 아래로 갈수록 높아지는 층을 이루고 물은 움직이지 않는다.

곁주석

0장의 상자(480×270 격자, 원기둥 하나, 스텝 1/120초)에 아래로 1의 체적력을 주고 엔진을 10초 돌려 보았다. 압력 $p/\rho$ 는 높이 1을 내려갈 때마다 1씩 높아졌고(0.001% 안에서), 물은 거의 가만히 있었다. 가장 빠른 물도 10초 내내 0.008을 넘지 않았다. 한 스텝 동안 힘이 보태는 속도(힘 1에 스텝 1/120초를 곱한 값)에 가까운 크기다. 스텝을 1/60초로 하면 0.017, 1/240초로 하면 0.004로, 스텝에 비례해 바뀌었다. 같은 상자에서 압력 투영을 끄면 1.5초 만에 가장 빠른 물이 1.4의 빠르기로 바닥을 향해 몰렸고, 그 뒤로 더 빨라지지는 않았다. 부력은 물의 온도나 소금기에 따라 정해지는 힘이다. 이 엔진의 잉크는 흐름을 바꾸지 않는 표시일 뿐이라서(0장) 부력을 넣지 않았다.

항을 하나씩 끄기

아래 실험은 0장의 실험실을 다시 꺼낸 것이다. 두 가지가 다르다. 하나는 물이 왼쪽에서 들어와 오른쪽으로 나간다는 점이다. 그래서 무대에 놓은 장애물마다 뒤로 꼬리가 생긴다. 다른 하나는 항마다 스위치가 있다는 점이다. 외력, 이류, 점성, 압력의 네 스위치로 항을 하나씩 끄고 켤 수 있고, 실험 창 왼쪽 위의 칩은 언제나 지금 켜져 있는 항만 보여 준다. 대본은 원기둥 앞을 같은 손길로 몇 번 저으면서 항을 하나씩 끈다. 그다음은 직접 해 볼 차례다. 마우스로 끌면(휴대폰에서는 손가락으로 쓸면) 물을 젓고, ‘손’을 ‘그리기’로 바꾸면 끌고 간 길을 따라 장애물이 그려진다. ‘장애물’로 원기둥, 판, 콩 모양을 고를 수도 있고, ‘젓기’로 장애물을 끌어 옮길 수도 있다. 장애물을 바꾸어도 그린 선은 남는다.

곁주석

값은 모두 무차원이다. 물은 0.3의 빠르기로 들어오고 원기둥의 지름은 0.16이다. 점성 $\nu$ 의 처음 값 0.0005에서 원기둥의 레이놀즈 수는 100쯤이고, ‘점성 ν’로 $\nu$ 를 0.0001에서 0.002까지(Re로 480에서 24까지) 바꿀 수 있다. 압력 스위치 하나가 $\prs{-\nabla p/\rho}$ 와 $\cont{\nabla\cdot\mathbf{u} = 0}$ 을 함께 끄고 켜는 것은, 물이 몰리거나 비지 않게 지키는 것이 바로 압력이 하는 일이기 때문이다(3장). 잉크는 어느 항을 끄든 흐름에 실려 다닌다. 엔진이 이류 항과 따로 잉크를 스텝마다 옮기기 때문이다. WebGPU가 없는 브라우저에서는 녹화한 영상이 나오고, 영상에서는 스위치를 바꾸거나 저을 수 없다.

처음에는 모든 항이 켜져 있다. 대본이 원기둥 앞을 아래에서 위로 저으면 바깥 힘이 물을 밀어 소용돌이를 만들고, 그 소용돌이는 흐름에 실려 원기둥을 돌아 떠내려간다.

외력을 끄면 같은 손길로 저어도 아무 일도 일어나지 않는다. 물은 밀리지 않고 잉크도 생기지 않는다. 엔진에서는 손의 잉크도 이 단계가 넣기 때문이다. 0장에서 물을 저은 손이 바로 이 항이었다.

대본은 이류를 끄면서 와도 층을 켠다. 주황과 하늘색은 물이 도는 방향이고, 같은 색의 칩과는 상관이 없다. 저어서 생긴 소용돌이가 흐름에 실려 가지 않고 저은 자리에 선다. 3초 동안 그 중심이 옮겨 간 거리는 0.013뿐이었다. 이류를 켠 채 같은 손길로 저으면 소용돌이는 2초 만에 0.51을 떠내려간다. 그래도 잉크와 입자는 멈추지 않고 서 있는 소용돌이 사이로 흘러 지나간다. 5장의 소용돌이 쌍이 이류 항을 끄자 그 자리에 멈춘 것과 같다. 흐름이 다른 것을 실어 나르는 일은 그대로이고, 흐름이 자기 속도를 실어 나르는 일만 꺼진 것이다.

이류를 다시 켜고 점성을 0.0005에서 0.002로 네 배 올리면, 원기둥의 Re가 100쯤에서 24로 내려간다. 같은 손길로 저어 만든 소용돌이가 금방 잦아든다. 저은 지 2초 뒤에 가장 센 곳의 와도는 처음의 10%였다. 처음 값 0.0005의 점성에서는 26%가 남는다. 이웃끼리 속도를 고르게 나누는 점성이 소용돌이의 도는 움직임을 둘레로 나누어 퍼뜨리기 때문이다(4장).

점성을 끄면 소용돌이가 더 오래 남는다. 2초 뒤에도 37%가 남는다. 원기둥에서 떨어져 나가는 층도 더 잘게 말려 소용돌이를 흘린다. 그래도 점성을 끈 물에서조차 소용돌이는 조금씩 약해진다. 점성이 없는 평면 흐름이라면 물 덩어리마다 제 와도를 그대로 지니고 다니므로, 가장 센 곳도 약해지지 않아야 한다. 이것은 5장에서 본 엔진의 수치 오차, 한 스텝을 단계로 나누어 계산하는 데서 오는 오차다(곁주석).

곁주석

대본의 손길로 점성을 끈 물을 저으면, 2초 뒤 가장 센 곳의 와도가 1/60초 스텝(이 실험실)에서 37%, 1/120초 스텝(0장의 실험실)에서 50%, 1/240초 스텝에서 62% 남았다. 격자를 가로세로 절반으로 줄이면 40%, 두 배로 늘리면 33%였다. 스텝의 길이가 남는 몫을 바꾸는 폭이 격자가 바꾸는 폭보다 훨씬 크고, 격자를 촘촘히 해도 남는 몫이 늘지 않는다. 그래서 5장에서 본 대로, 한 스텝을 단계로 나누는 데서 오는 오차로 읽는다.

압력을 끄면 원기둥 뒤가 비고, 빈 자리의 가장자리에 입자가 몰려 밝은 줄이 생긴다. 2장에서 투영을 껐을 때처럼 물이 한쪽에서는 몰리고 다른 쪽에서는 빠져나간다. 압력이 있을 때는 원기둥 앞에서 압력이 높아져 물을 미리 옆으로 비켜 보내고, 뒤에서는 낮아져 둘레의 물을 끌어들인다(3장). 압력이 없으니 물은 원기둥에 막히기만 할 뿐 돌아가지 못하고, 원기둥 뒤를 채우는 물이 없다. 압력을 다시 켜면 발산이 곧 다시 0에 가까워지고, 몇 초 뒤에는 처음처럼 원기둥 뒤에 꼬리가 다시 생긴다.

곁주석

압력을 끈 지 3.5초 뒤, 칸의 크기 $h$ 를 곱하고 들어오는 빠르기 $U$ 로 나눈 발산은 −0.68에서 0.09 사이였고, 다시 켠 지 1.5초 뒤에는 0.000005 안이었다. 스위치의 조합 16가지를 저마다 30초씩 저으며 돌려도 흐름은 터지지 않았다. 가장 빠른 것은 외력만 켠 조합으로 들어오는 빠르기의 26배였다. 엔진은 50배를 넘으면 터진 것으로 본다.

이제 직접 해 보자. 항을 켜고 끄며 저어 보고, ‘손’을 ‘그리기’로 바꾸어 흐름을 가로지르는 선을 그어 보자. 선이 지나간 칸이 고체가 되고, 물은 그 둘레로 돌아간다. 원기둥을 없애고 흐름을 세로로 가로지르는 선을 그었더니, 3초 뒤 선 바로 뒤에서는 물이 거꾸로 흐르고(들어오는 빠르기의 0.4배쯤), 선 끝을 돌아 나가는 물은 들어오는 물의 1.8배로 빨라졌다.

곁주석

그린 선은 굵기가 0.03이고, 한 번 그을 때 점을 48개까지 받는다. 선은 여섯 개까지 남고, 일곱째를 그으면 가장 먼저 그은 선이 사라진다. ‘처음부터’를 누르면 모두 지워진다. 키보드만 쓸 때는 ‘장애물’로 모양을 고른다. 녹화한 영상에서는 그릴 수 없다.

날개는 어떻게 뜨는가

비행기의 날개는 공기를 가르며 앞으로 나아갈 뿐인데 무거운 비행기를 떠받친다. 그 힘이 양력이다. 아래 실험에서는 물이 왼쪽에서 날개를 향해 흘러온다. 날개가 공기 속을 나아가는 것과 같은 일을, 날개를 제자리에 두고 물을 흘려서 본다. 날개의 단면은 1930년대에 NACA(미국 항공자문위원회)가 식으로 정해 둔 날개꼴 가운데 하나인 NACA 0006이다. 위아래가 대칭이고, 가장 두꺼운 곳이 앞전에서 시위(앞전과 뒷전을 잇는 길이)의 30% 되는 곳에 있으며 두께는 시위의 6%다. 날개가 다가오는 흐름에 대해 기운 각도를 받음각이라 하고, ‘받음각’으로 0도, 4도, 8도, 16도를 고를 수 있다. 처음은 8도다. 3초마다 흰 잉크 띠 하나가 날개 앞을 가로질러 놓인다. 같은 순간에 출발한 물을 표시하는 띠다.

곁주석

값은 모두 무차원이다. 물은 0.5의 빠르기 $U$ 로 들어오고, 날개의 시위 $c$ 는 0.4, 점성 $\nu$ 는 0.0002다. 시위로 잰 레이놀즈 수 $Uc/\nu$ 는 1000이다. 실제 비행기 날개는 이것의 백 배에서 만 배인 $10^5$ 에서 $10^7$ 쯤의 Re로 난다. Re가 낮아서 날개 표면에 붙어 느려진 물의 층이 두껍다. 흐름은 2차원이라 끝없이 긴 날개의 한 단면이고, 날개 끝에서 말려 나가는 소용돌이는 없다. 데스크톱의 480×270 격자에서 시위는 108칸, 두께는 6.5칸이고, 휴대폰의 240×135 격자에서는 그 절반이다. 날개는 칸을 고체로 표시해 그리므로 표면이 계단 모양이다. 위아래의 미끄러운 벽은 날개의 회전축(시위의 4분의 1 되는 곳)에서 0.5씩 떨어져 있어서, 끝없이 넓은 공기 속의 날개와 꼭 같지는 않다. 흰 띠는 높이 0.3에서 0.7까지, 날개 앞 0.38쯤에 한꺼번에 놓는 잉크 점 열여섯 개이고, 물을 밀지 않는다. 날개를 끌어 옮기거나 물을 저을 수는 없다. WebGPU가 없는 브라우저에서는 녹화한 영상이 나오고, 영상에서는 받음각을 바꿀 수 없다.

물이 흐르기 시작하면 뒷전에서 소용돌이 하나가 말려 떨어져 나가 흘러간다. 와도 층에서 주황, 곧 시계 반대 방향으로 도는 소용돌이다. 이 소용돌이가 떠나면서 날개 둘레에는 반대 방향, 곧 시계 방향으로 도는 흐름이 남는다. 물이 흐르기 시작한 지 0.5초 뒤, 날개를 두른 상자의 둘레에서 잰 순환은 시계 방향으로 0.038이었고, 날개 뒤를 두른 상자의 둘레에서 잰 순환은 시계 반대 방향으로 0.038이었다. 크기가 같고 방향이 반대다. 날개 둘레의 이 시계 방향 흐름은 날개 위의 물에는 흐름 방향으로, 아래의 물에는 흐름을 거스르는 방향으로 보태져서, 위의 흐름을 빠르게 하고 아래의 흐름을 느리게 한다. 처음 소용돌이를 출발 소용돌이라 한다.

곁주석

순환은 5장에서 소용돌이의 세기를 잰 값이다. 고리를 따라 물이 흐르는 빠르기의 평균에 고리의 길이를 곱한다. 점성이 없는 유체에서는 물과 함께 움직이는 고리 위의 순환이 시간이 지나도 바뀌지 않는다. 켈빈의 순환 정리다. 처음의 고른 흐름에서는 어느 고리의 순환도 0이었으니, 출발한 뒤에도 날개 둘레의 순환과 출발 소용돌이의 순환을 더하면 0이어야 한다. 이 실험은 점성이 있지만, 물이 흐르기 시작한 지 0.5초 뒤에 잰 두 순환은 0.1% 안에서 크기가 같았다.

이제 흰 띠를 따라가 보자. 띠는 앞전에서 날개의 위아래로 갈라진다. 날개 위로 지나간 물이 먼저 뒷전에 닿는다. 그때 날개 아래로 지나간 물은 아직 시위의 83% 되는 곳에 있다. 띠를 놓은 지 1.42초 만에 위쪽 물이, 1.55초 만에 아래쪽 물이 뒷전에 닿았다. 0.13초 차이다. 아래쪽 물은 날개를 지난 뒤에도 끝내 따라잡지 못해서, 뒷전에서 0.3 지난 곳에서도 여전히 0.13초 늦다. 같이 출발한 물이 같이 도착하지 않는다.

이 차이는 날개 위의 흐름이 아래보다 빠르기 때문이다. 시위 가운데에서 표면으로부터 0.1 떨어진 곳의 빠르기는 위가 $U$ 의 1.16배, 아래가 1.03배였다. 압력 층으로 바꾸어 보면 날개 위는 압력이 낮고(하늘색) 아래는 높다(주황). 같은 곳의 압력 $p/\rho$ 는 위가 −0.035, 아래가 0.007이었다. 날개 둘레에서 압력이 가장 높은 곳은 앞전 바로 아래, 흐름이 날개에 부딪혀 멈추는 곳이고(0.12), 가장 낮은 곳은 앞전 바로 위다(−0.14). 날개의 표면 전체에 걸쳐 압력이 미는 힘을 더하면 위를 향한다. 이것이 양력이다. 양력은 다가오는 흐름의 $\tfrac12\rho U^2$ 과 시위로 나눈 양력 계수로 견준다.

C_L = \frac{L}{\tfrac{1}{2}\rho U^2 c}

$L$ 은 날개의 폭 1당 양력이다. 이 실험의 8도에서 $C_L$ 은 0.49였다.

곁주석

얇은 날개 이론은 높은 Re에서 얇은 날개의 양력 계수를 받음각(라디안)의 $2\pi$ 배로 어림하고, 8도라면 0.88이다. 이 실험의 0.49는 그보다 작다. Re 1000에서 표면의 느린 층이 두꺼운 것이 한 까닭으로 보이지만, 가까운 벽과 계단 표면도 얽혀 있어서 이 실험이 그 몫을 따로 가려내지는 않았다. 휴대폰의 성긴 격자에서는 0.42였고, 그래도 위쪽 물이 0.1초 먼저 뒷전에 닿았다.

빠른 흐름과 낮은 압력이 함께 나타나는 것은 3장의 베르누이 식 그대로다. 베르누이 식은 네 조건, 곧 정상, 비점성, 비압축, 같은 유선 위라는 조건에서 맞는다. 물이 흐르기 시작한 지 몇 초가 지나 흐름이 자리를 잡은 뒤이고(정상), 표면에 붙은 느린 층 밖이며(점성이 거의 일하지 않는다), 물의 밀도는 일정하다(비압축). 넷째 조건도 맞는다. 날개 위와 아래를 지나는 유선은 모두 날개 앞의 고른 흐름에서 왔고, 거기서는 모든 유선의 압력과 빠르기가 같다. 그래서 상수 $C$ 가 모든 유선에서 같고, 위아래를 서로 견줄 수 있다. 실제로 위, 아래, 날개 앞에서 잰 $p/\rho + \tfrac12 |\mathbf{u}|^2$ 은 날개 앞의 값과 $\tfrac12 U^2$ 의 2.1%(위), 2.4%(아래)만큼 달랐다(위와 아래 사이는 4.5%). 위아래의 압력 자체는 $\tfrac12 U^2$ 의 34%만큼 달랐으니, 이 차이들은 그에 견주면 작다.

날개 뒤의 흐름은 아래로 꺾인다. 뒷전에서 0.2 뒤, 뒷전 높이에서 0.3 위까지의 물은 평균 $U$ 의 0.035배로 아래를 향해 흘렀다. 날개가 물을 아래로 꺾었다는 뜻이고, 뉴턴의 제3법칙(작용과 반작용)에 따라 물은 날개를 위로 민다. 압력과 꺾인 흐름은 양력의 서로 다른 두 설명이 아니라 한 흐름의 두 모습이다. 날개 위아래의 압력 차는 날개를 위로 미는 동시에 둘레의 물을 아래로 꺾는다. 물이 날개를 미는 힘이 그 압력이고, 날개가 물을 미는 힘이 그 반작용이다.

받음각을 0도로 바꾸면 날개가 바뀐 흐름에 맞추느라 소용돌이를 하나 흘리고, 몇 초 뒤에 흐름이 다시 자리를 잡는다. 기울지 않은 대칭 날개는 위와 아래가 거울에 비친 것처럼 같다. 띠는 위아래로 함께 뒷전에 닿고(한 프레임, 1/60초 안에서), 양력 계수는 0.001보다 작다. 양력이 없다. 날개를 기울여야 양력이 생긴다.

16도로 기울이면 흐름이 위 표면을 끝까지 따라가지 못하고 떨어져 나가고, 날개 뒤로 큰 소용돌이를 흘린다. 날개 위 표면의 뒤쪽 절반 바로 위(표면에서 0.06 안)에서 가장 느린 곳의 흐름 방향 속도를 재 보면, 4도에서는 아직 앞으로 흘렀지만 $U$ 의 0.0035배로 거의 멎어 있었고, 8도에서는 $U$ 의 0.035배로 거꾸로 흘렀으며, 16도에서는 $U$ 의 0.42배로 거꾸로 흘렀다(16도의 값은 기울인 지 3.4초 뒤의 것이고, 소용돌이를 흘리는 동안 시간에 따라 바뀐다). 실제 날개는 받음각을 너무 키우면 흐름이 위 표면에서 떨어져 나가 양력이 크게 줄어든다. 이를 실속이라 한다.

오해 하나

양력을 두고 널리 퍼진 설명이 있다. 날개 위쪽이 더 볼록해서 길이가 길고, 앞전에서 갈라진 물이 뒷전에서 다시 만나야 하므로 위쪽 물이 더 빨리 가야 한다는 설명이다. 이 설명은 틀렸다. 무엇보다 이 실험에서 같이 출발한 물은 같이 도착하지 않았다. 위로 간 물이 먼저 닿았고, 아래로 간 물은 끝내 따라잡지 못했다. 둘이 뒷전에서 다시 만나야 할 까닭은 식 어디에도 없다. 설명이 기대는 길이의 차이도 맞지 않는다. 이 실험의 NACA 0006은 위아래 표면의 길이가 똑같은 대칭 날개라서, 그 설명대로라면 8도에서도 양력이 없어야 한다. 그런데 양력 계수가 0.49였다. 종이비행기의 평평한 날개도 날고, 비행기는 뒤집혀서도 난다. 날개 위의 흐름이 아래보다 빠르고 압력이 낮다는 것은 맞다. 그 빠르기는 실제로 이 설명이 내놓는 것보다 훨씬 커서, 같은 시간에 도착한다고 셈한 빠르기로는 실제보다 훨씬 작은 양력이 나온다. 빠르기를 정하는 것은 날개 둘레의 흐름 전체, 곧 날개가 흐름을 아래로 꺾는 방식이다.

곁주석

NASA 글렌 연구소의 항공 입문서 ‘Equal Transit Theory Interactive’ 쪽이 이 설명의 문제를 같은 근거로 정리해 두었다. 대칭 날개와 평평한 날개도 양력을 내고, 뒤집혀 나는 비행기가 있으며, 앞전에서 가까이 있던 공기 분자가 뒷전에서 함께 있지 않고, 이 설명으로 셈한 양력은 실제보다 훨씬 작다는 것이다.

엔진의 한 스텝

0장에서 말한 대로, 쌓아 온 항들이 그대로 이 실험실의 엔진이 된다. 엔진은 시간을 짧게 끊어 한 스텝씩 나아가고(이 장의 실험에서는 한 스텝이 1/60초다), 스텝마다 단계 다섯을 정해진 차례로 돈다. ①에서 ④까지는 식의 항을 맡고(④는 두 항을 함께 맡는다), ⑤는 무대의 가장자리와 장애물에서 속도를 정한다. 아래 그림이 그 한 스텝이다. 단계를 고르면 위 줄에서 그 단계가 읽기만 하는 값은 점선 테두리로, 쓰는 값은 채워서 표시된다. 읽고 나서 그 자리에 다시 쓰는 값도 채워진다. 각 단계의 실제 코드는 그림 아래에 같은 번호로 실었다. 코드는 그래픽 장치에서 돌아가는 WGSL(WebGPU의 셰이더 언어)로 쓰여 있고, 여기 실린 것은 엔진이 쓰는 파일 그대로다.

곁주석

엔진은 속도를 칸의 가장자리(면)에, 압력과 잉크를 칸의 가운데에 둔다. 가로 속도는 칸의 왼쪽과 오른쪽 면에, 세로 속도는 아래와 위 면에 있다. 이렇게 엇갈려 놓으면 칸 하나로 들어오고 나가는 양을 면의 속도로 바로 셀 수 있다(2장의 상자). 그림은 한 스텝을 5초에 한 바퀴로 늘려 보여 준다. 경계 단계는 다른 단계 뒤마다 한 번씩 돌지만 그림에는 한 번만 그렸다. 꺼진 항의 단계는 건너뛰고(②, ③, ④는 그 뒤의 경계 단계까지 함께), 잉크를 옮기는 일은 언제나 돈다.

① 외력

바깥 힘 $\frc{\mathbf{f}}$ 를 더한다. 물 전체에 미는 힘(체적력)에 한 스텝의 시간을 곱해 모든 면의 속도에 더하고, 젓는 손이 미는 힘과 잉크는 손이 지난 자리를 가운데로 한 종 모양으로 더한다. 장애물이나 무대의 가장자리에 닿은 면은 건드리지 않는다. 그 면의 속도는 ⑤ 경계 단계가 정한다.

\frc{\mathbf{f}}

forces.wgsl

src/lab/engine/shaders/forces.wgsl

// Stage 1 · external force f: body force plus Gaussian "splats" (stirring) of velocity and dye.
struct Splat {
  pos: vec2f,
  vel: vec2f,
  dye: vec4f,
  radius: f32,
  pad0: f32,
  pad1: f32,
  pad2: f32,
};

@group(0) @binding(1) var<storage, read_write> uF: array<f32>;
@group(0) @binding(2) var<storage, read_write> vF: array<f32>;
@group(0) @binding(3) var<storage, read_write> dyeF: array<vec4f>;
@group(0) @binding(4) var<storage, read> splats: array<Splat>;
@group(0) @binding(5) var<storage, read> solid: array<u32>;

fn splatWeight(p: vec2f, s: Splat) -> f32 {
  var d = p - s.pos;
  let L = domainSize();
  if (periodicX()) { d.x = d.x - L.x * round(d.x / L.x); }
  if (periodicY()) { d.y = d.y - L.y * round(d.y / L.y); }
  return exp(-dot(d, d) / (s.radius * s.radius));
}

@compute @workgroup_size(8, 8)
fn forces(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let p = vec2f(f32(i), f32(j) + 0.5) * P.h;
    var du = P.gx * P.dt;
    for (var s = 0u; s < P.splatCount; s++) { du += splats[s].vel.x * splatWeight(p, splats[s]); }
    uF[iU(i, j)] += du;
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let p = vec2f(f32(i) + 0.5, f32(j)) * P.h;
    var dv = P.gy * P.dt;
    for (var s = 0u; s < P.splatCount; s++) { dv += splats[s].vel.y * splatWeight(p, splats[s]); }
    vF[iV(i, j)] += dv;
  }
  if (i < P.nx && j < P.ny && solid[iC(i, j)] == 0u) {
    let p = (vec2f(f32(i), f32(j)) + 0.5) * P.h;
    var dd = vec4f(0.0);
    for (var s = 0u; s < P.splatCount; s++) { dd += splats[s].dye * splatWeight(p, splats[s]); }
    dyeF[iC(i, j)] = clamp(dyeF[iC(i, j)] + dd, vec4f(0.0), vec4f(1.0));
  }
}

② 이류

흐름이 제 속도를 실어 나르는 $\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}}$ 를 푼다. 면마다, 이 스텝 동안 그 자리에 온 물이 어디에서 왔는지 흐름을 거슬러 짚고, 그곳의 속도를 옮겨 온다(advect_vel.wgsl). 짚을 때는 반 스텝 거슬러 간 곳의 속도로 다시 한 번 짚는다. 같은 셰이더를 거꾸로 한 번 더 돌려 옮겨 온 값을 흐름을 따라 되돌려 보면, 처음 값과 조금 어긋난다. 둘째 셰이더 correct_vel.wgsl이 그 어긋난 만큼의 절반을 고친다. 고친 값이 옮겨 올 때 쓴 둘레 값의 범위를 벗어나면 고치기 전의 값을 쓴다. 잉크는 같은 방법으로, 제 몫의 두 셰이더로 옮긴다(advect_dye.wgsl, correct_dye.wgsl). 다만 고친 값이 그 범위를 벗어나면 고치기 전의 값 대신 범위의 끝 값으로 자른다. 잉크는 스텝의 맨 끝, 압력 투영까지 마친 흐름으로 옮기고, 이류 항을 꺼도 언제나 옮긴다. 이 실험실에서 이류를 꺼도 잉크가 흐른 까닭이다.

\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}}

advect_vel.wgsl · correct_vel.wgsl · advect_dye.wgsl · correct_dye.wgsl

src/lab/engine/shaders/advect_vel.wgsl

// Stage 2 · inertia (u.grad)u, part 1: semi-Lagrangian transport of the velocity field.
// DT_SIGN = +1: source u^n traced backward -> hat.  DT_SIGN = -1: source hat traced forward -> back.
// Both passes trace through u^n (uIn / vIn), as MacCormack requires.
override DT_SIGN: f32 = 1.0;

@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> srcU: array<f32>;
@group(0) @binding(4) var<storage, read> srcV: array<f32>;
@group(0) @binding(5) var<storage, read_write> dstU: array<f32>;
@group(0) @binding(6) var<storage, read_write> dstV: array<f32>;
@group(0) @binding(7) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn advect_vel(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  let dts = DT_SIGN * P.dt;
  if (i <= P.nx && j < P.ny) {
    let k = iU(i, j);
    if (uFixed(i, j)) {
      dstU[k] = srcU[k];
    } else {
      dstU[k] = sampleU_src(trace(vec2f(f32(i), f32(j) + 0.5) * P.h, dts));
    }
  }
  if (i < P.nx && j <= P.ny) {
    let k = iV(i, j);
    if (vFixed(i, j)) {
      dstV[k] = srcV[k];
    } else {
      dstV[k] = sampleV_src(trace(vec2f(f32(i) + 0.5, f32(j)) * P.h, dts));
    }
  }
}

src/lab/engine/shaders/correct_vel.wgsl

// Stage 2 · inertia, part 2: MacCormack correction  u_new = hat + (u^n - back) / 2,
// falling back to hat where the result would leave the range of the values it was interpolated from.
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read_write> hatU: array<f32>;
@group(0) @binding(4) var<storage, read_write> hatV: array<f32>;
@group(0) @binding(5) var<storage, read> backU: array<f32>;
@group(0) @binding(6) var<storage, read> backV: array<f32>;
@group(0) @binding(7) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn correct_vel(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let k = iU(i, j);
    let r = rangeU(trace(vec2f(f32(i), f32(j) + 0.5) * P.h, P.dt));
    let c = hatU[k] + 0.5 * (uIn[k] - backU[k]);
    hatU[k] = select(c, hatU[k], c < r.x || c > r.y);
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let k = iV(i, j);
    let r = rangeV(trace(vec2f(f32(i) + 0.5, f32(j)) * P.h, P.dt));
    let c = hatV[k] + 0.5 * (vIn[k] - backV[k]);
    hatV[k] = select(c, hatV[k], c < r.x || c > r.y);
  }
}

src/lab/engine/shaders/advect_dye.wgsl

// Dye transport (visualisation): the same semi-Lagrangian trace as the velocity, cell-centred RGBA.
override DT_SIGN: f32 = 1.0;

@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> dyeSrc: array<vec4f>;
@group(0) @binding(4) var<storage, read_write> dyeDst: array<vec4f>;

@compute @workgroup_size(8, 8)
fn advect_dye(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let p = (vec2f(f32(id.x), f32(id.y)) + 0.5) * P.h;
  dyeDst[iC(id.x, id.y)] = sampleDye(trace(p, DT_SIGN * P.dt));
}

src/lab/engine/shaders/correct_dye.wgsl

// Dye MacCormack correction, clamped per channel to the interpolation range (no new extrema).
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> dyeSrc: array<vec4f>;
@group(0) @binding(4) var<storage, read_write> dyeHat: array<vec4f>;
@group(0) @binding(5) var<storage, read> dyeBack: array<vec4f>;

@compute @workgroup_size(8, 8)
fn correct_dye(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let k = iC(id.x, id.y);
  let p = (vec2f(f32(id.x), f32(id.y)) + 0.5) * P.h;
  let r = rangeDye(trace(p, P.dt));
  dyeHat[k] = clamp(dyeHat[k] + 0.5 * (dyeSrc[k] - dyeBack[k]), r.lo, r.hi);
}

③ 점성

이웃과 속도를 고르게 나누는 $\vis{\nu\nabla^2\mathbf{u}}$ 를 푼다. 엔진은 스텝이 끝난 뒤의 새 속도로 이웃과의 차이를 셈하는 식을 푼다. 새 속도가 식의 양쪽에 다 있어서 한 번에 풀리지 않으므로, 칸을 바둑판처럼 두 색으로 나누어 한 색씩 번갈아 고치기를 되풀이한다. 되풀이할 횟수는 $\nu\,\Delta t/h^2$ 이 클수록 많이 필요하다. 이 장의 실험실은 점성에 따라 4번에서 25번(처음 값에서는 7번)을 되풀이한다.

\vis{\nu\nabla^{2}\mathbf{u}}

diffuse.wgsl

src/lab/engine/shaders/diffuse.wgsl

// Stage 3 · viscosity nu*lap(u): one red-black Gauss-Seidel sweep of the implicit step
//   (1 + 4a) u - a * sum(u_nb) = u0,   a = nu * dt / h^2
// Neighbours beyond walls come from the ghost rules in fetchU_d / fetchV_d.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> uD: array<f32>;
@group(0) @binding(2) var<storage, read_write> vD: array<f32>;
@group(0) @binding(3) var<storage, read> u0: array<f32>;
@group(0) @binding(4) var<storage, read> v0: array<f32>;
@group(0) @binding(5) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn diffuse(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (((i + j) & 1u) != COLOR) { return; }
  let a = P.nu * P.dt / (P.h * P.h);
  let ii = i32(i);
  let jj = i32(j);
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let s = fetchU_d(ii - 1, jj) + fetchU_d(ii + 1, jj) + fetchU_d(ii, jj - 1) + fetchU_d(ii, jj + 1);
    uD[iU(i, j)] = (u0[iU(i, j)] + a * s) / (1.0 + 4.0 * a);
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let s = fetchV_d(ii - 1, jj) + fetchV_d(ii + 1, jj) + fetchV_d(ii, jj - 1) + fetchV_d(ii, jj + 1);
    vD[iV(i, j)] = (v0[iV(i, j)] + a * s) / (1.0 + 4.0 * a);
  }
}

④ 압력 투영

$\prs{-\nabla p/\rho}$ 와 $\cont{\nabla\cdot\mathbf{u} = 0}$ 을 함께 푼다. 먼저 칸마다 면의 속도로 들어오고 나가는 양을 세어 발산을 잰다(divergence.wgsl). 그다음 그 발산을 지울 압력을 푼다. 3장에서 압력이 어디서나 한꺼번에 맞춰진다고 한 것은, 칸마다 이웃 칸의 압력과 이어진 식 하나씩을 모든 칸에서 함께 풀어야 한다는 뜻이다. 엔진은 이 식을 격자를 절반씩 성기게 줄여 가며 푸는 것이 기본이다. 촘촘한 격자에서 잔물결 같은 오차를 줄이고, 남은 오차를 성긴 격자로 옮겨 큰 무늬의 오차를 싸게 줄인 뒤, 그 고침을 다시 촘촘한 격자로 가져온다(다중 격자법, mg_로 시작하는 다섯 파일). 설정이 한 격자에서 되풀이해 푸는 방식일 때나, 격자를 절반씩 줄여도 가장 성긴 격자가 충분히 작아지지 않을 때(칸 4096개 이하가 되지 않을 때)는 pressure.wgsl이 한 격자에서 되풀이해 푼다. 마지막으로 gradient.wgsl이 압력의 기울기만큼 속도를 빼서, 물이 몰리거나 비지 않는 흐름을 만든다.

\mathbf{u} \leftarrow \mathbf{u} - \Delta t\,\frac{\nabla p}{\rho}

\prs{-\frac{\nabla p}{\rho}},\quad \cont{\nabla\cdot\mathbf{u} = 0}

divergence.wgsl · mg_smooth.wgsl · mg_residual.wgsl · mg_restrict.wgsl · mg_coarsest.wgsl · mg_prolong.wgsl · pressure.wgsl · gradient.wgsl

src/lab/engine/shaders/divergence.wgsl

// Stage 4a · mass conservation check: net outflow of each fluid cell, (u_E - u_W + v_N - v_S) / h.
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;
@group(0) @binding(4) var<storage, read_write> divOut: array<f32>;

@compute @workgroup_size(8, 8)
fn divergence(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let i = id.x;
  let j = id.y;
  let k = iC(i, j);
  if (solid[k] != 0u) {
    divOut[k] = 0.0;
    return;
  }
  divOut[k] = (uIn[iU(i + 1u, j)] - uIn[iU(i, j)] + vIn[iV(i, j + 1u)] - vIn[iV(i, j)]) / P.h;
}

src/lab/engine/shaders/mg_smooth.wgsl

// Multigrid smoother: one red-black Gauss-Seidel half-sweep of  sum(p_nb - p) / h^2 = f  on one level.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;

@compute @workgroup_size(8, 8)
fn mg_smooth(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny || ((id.x + id.y) & 1u) != COLOR) { return; }
  let k = id.y * L.nx + id.x;
  if (mask[k] != 0u) { return; }
  let s = levelStencil(i32(id.x), i32(id.y));
  if (s.y == 0.0) { return; } // sealed off on all sides: nothing to balance
  p[k] = (s.x - L.h * L.h * L.fScale * f[k]) / s.y;
}

src/lab/engine/shaders/mg_residual.wgsl

// Multigrid residual  r = f - sum(p_nb - p) / h^2  on one level (0 in solid and sealed cells).
@group(0) @binding(1) var<storage, read> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;
@group(0) @binding(4) var<storage, read_write> r: array<f32>;

@compute @workgroup_size(8, 8)
fn mg_residual(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny) { return; }
  let k = id.y * L.nx + id.x;
  let s = levelStencil(i32(id.x), i32(id.y));
  if (mask[k] != 0u || s.y == 0.0) {
    r[k] = 0.0;
    return;
  }
  r[k] = L.fScale * f[k] - (s.x - s.y * p[k]) / (L.h * L.h);
}

src/lab/engine/shaders/mg_restrict.wgsl

// Multigrid restriction: a coarse cell's right-hand side is the mean residual of its fine children inside
// the grid (one or two of them on the last row or column of an odd side); its correction starts at 0.
@group(0) @binding(1) var<storage, read> r: array<f32>;
@group(0) @binding(2) var<storage, read_write> fc: array<f32>;
@group(0) @binding(3) var<storage, read_write> pc: array<f32>;

@compute @workgroup_size(8, 8)
fn mg_restrict(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.cnx || id.y >= L.cny) { return; }
  var sum = 0.0;
  var count = 0.0;
  for (var b = 0u; b < 2u; b++) {
    for (var a = 0u; a < 2u; a++) {
      let i = 2u * id.x + a;
      let j = 2u * id.y + b;
      if (i < L.nx && j < L.ny) {
        sum += r[j * L.nx + i];
        count += 1.0;
      }
    }
  }
  let k = id.y * L.cnx + id.x;
  fc[k] = sum / count;
  pc[k] = 0.0;
}

src/lab/engine/shaders/mg_coarsest.wgsl

// Multigrid coarsest level: L.iterations red-black SOR sweeps by one workgroup of 256 invocations, with a
// storage barrier between colours so every cell sees its neighbours' new values.
@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;

@compute @workgroup_size(256)
fn mg_coarsest(@builtin(local_invocation_index) t: u32) {
  let n = L.nx * L.ny;
  for (var it = 0u; it < L.iterations; it++) {
    for (var color = 0u; color < 2u; color++) {
      for (var k = t; k < n; k += 256u) {
        let i = k % L.nx;
        let j = k / L.nx;
        if (((i + j) & 1u) == color && mask[k] == 0u) {
          let s = levelStencil(i32(i), i32(j));
          if (s.y > 0.0) {
            p[k] = mix(p[k], (s.x - L.h * L.h * L.fScale * f[k]) / s.y, L.omega);
          }
        }
      }
      storageBarrier();
    }
  }
}

src/lab/engine/shaders/mg_prolong.wgsl

// Multigrid prolongation: each fluid fine cell adds the bilinear interpolation of the coarse correction,
// its weights renormalised over the fluid coarse cells, so nothing leaks in from a solid one. The parent
// cell always carries a weight of at least 9/16, so the weights never sum to 0.
@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> pc: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;
@group(0) @binding(4) var<storage, read> maskC: array<u32>;

// Coarse cell index, or -1 beyond a non-periodic edge or inside a solid.
fn coarseCell(I: i32, J: i32) -> i32 {
  let cnx = i32(L.cnx);
  let cny = i32(L.cny);
  var ii = I;
  var jj = J;
  if (I < 0 || I >= cnx) {
    if (L.edgeL != LEVEL_PERIODIC) { return -1; }
    ii = ((I % cnx) + cnx) % cnx;
  }
  if (J < 0 || J >= cny) {
    if (L.edgeB != LEVEL_PERIODIC) { return -1; }
    jj = ((J % cny) + cny) % cny;
  }
  let k = jj * cnx + ii;
  if (maskC[u32(k)] != 0u) { return -1; }
  return k;
}

@compute @workgroup_size(8, 8)
fn mg_prolong(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny) { return; }
  let k = id.y * L.nx + id.x;
  if (mask[k] != 0u) { return; }
  let g = vec2f(f32(id.x), f32(id.y)) * 0.5 - 0.25; // fine cell centre in coarse index space
  let c0 = vec2i(floor(g));
  let t = g - vec2f(c0);
  var sum = 0.0;
  var weight = 0.0;
  for (var b = 0; b < 2; b++) {
    for (var a = 0; a < 2; a++) {
      let K = coarseCell(c0.x + a, c0.y + b);
      if (K >= 0) {
        let w = select(1.0 - t.x, t.x, a == 1) * select(1.0 - t.y, t.y, b == 1);
        sum += w * pc[u32(K)];
        weight += w;
      }
    }
  }
  if (weight > 0.0) { p[k] += sum / weight; }
}

src/lab/engine/shaders/pressure.wgsl

// Stage 4b · pressure: one red-black SOR sweep of  sum(p_nb - p) / h^2 = div(u*) / dt.
// Walls, inflow and obstacles are Neumann (no term); outflow edges are Dirichlet with ghost p = 0.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> pR: array<f32>;
@group(0) @binding(2) var<storage, read> divIn: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;

// Neighbour contribution: x = pressure, y = 1 when the neighbour enters the stencil.
fn neighbour(i: i32, j: i32, side: u32) -> vec2f {
  let outsideX = i < 0 || i >= i32(P.nx);
  let outsideY = j < 0 || j >= i32(P.ny);
  if ((outsideX && !periodicX()) || (outsideY && !periodicY())) {
    return vec2f(0.0, select(0.0, 1.0, side == EDGE_OUTFLOW));
  }
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  let k = iC(u32(ii), u32(jj));
  if (solid[k] != 0u) { return vec2f(0.0, 0.0); }
  return vec2f(pR[k], 1.0);
}

@compute @workgroup_size(8, 8)
fn pressure(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny || ((id.x + id.y) & 1u) != COLOR) { return; }
  let k = iC(id.x, id.y);
  if (solid[k] != 0u) { return; }
  let i = i32(id.x);
  let j = i32(id.y);
  let s = neighbour(i - 1, j, P.edgeL) + neighbour(i + 1, j, P.edgeR)
        + neighbour(i, j - 1, P.edgeB) + neighbour(i, j + 1, P.edgeT);
  if (s.y == 0.0) { return; } // a cell sealed off on all sides has nothing to balance
  let gaussSeidel = (s.x - (P.h * P.h / P.dt) * divIn[k]) / s.y;
  pR[k] = mix(pR[k], gaussSeidel, P.omega);
}

src/lab/engine/shaders/gradient.wgsl

// Stage 4c · pressure push: u -= (dt/h) * grad p on every free face, making the field divergence-free.
@group(0) @binding(1) var<storage, read_write> uG: array<f32>;
@group(0) @binding(2) var<storage, read_write> vG: array<f32>;
@group(0) @binding(3) var<storage, read> pIn: array<f32>;
@group(0) @binding(4) var<storage, read> solid: array<u32>;

// Pressure of cell (i, j); beyond a non-periodic edge the ghost pressure is 0 (used by outflow faces only).
fn pAt(i: i32, j: i32) -> f32 {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); } else if (i < 0 || i >= i32(P.nx)) { return 0.0; }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); } else if (j < 0 || j >= i32(P.ny)) { return 0.0; }
  return pIn[iC(u32(ii), u32(jj))];
}

fn uProjected(i: u32, j: u32) -> bool {
  if (!periodicX() && i == 0u) { return P.edgeL == EDGE_OUTFLOW && !isSolid(0, i32(j)); }
  if (!periodicX() && i == P.nx) { return P.edgeR == EDGE_OUTFLOW && !isSolid(i32(P.nx) - 1, i32(j)); }
  return !uFixed(i, j);
}

fn vProjected(i: u32, j: u32) -> bool {
  if (!periodicY() && j == 0u) { return P.edgeB == EDGE_OUTFLOW && !isSolid(i32(i), 0); }
  if (!periodicY() && j == P.ny) { return P.edgeT == EDGE_OUTFLOW && !isSolid(i32(i), i32(P.ny) - 1); }
  return !vFixed(i, j);
}

@compute @workgroup_size(8, 8)
fn gradient(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  let s = P.dt / P.h;
  if (i <= P.nx && j < P.ny && uProjected(i, j)) {
    uG[iU(i, j)] -= s * (pAt(i32(i), i32(j)) - pAt(i32(i) - 1, i32(j)));
  }
  if (i < P.nx && j <= P.ny && vProjected(i, j)) {
    vG[iV(i, j)] -= s * (pAt(i32(i), i32(j)) - pAt(i32(i), i32(j) - 1));
  }
}

⑤ 경계

무대의 가장자리와 장애물에 닿은 면의 속도를 정한다. 벽을 가로지르는 면은 0이라서 물이 벽을 뚫고 흐르지 못한다. 입구의 면은 들어오는 빠르기, 출구의 면은 바로 안쪽 면의 값을 따르고(다만 ④ 압력 투영 뒤에는 투영이 정한 값을 그대로 둔다), 장애물에 닿은 면은 장애물이 움직이는 속도를 따른다. 벽을 따라 미끄러지는 속도는 이 단계가 아니라, 다른 단계가 벽 너머의 값을 셈할 때 벽의 종류에 따라 정해진다. 이 단계는 다른 단계가 돌 때마다 그 뒤에 돈다. 다른 단계가 경계에서 무엇을 셈하든, 다음 단계는 언제나 경계가 맞춰진 속도에서 시작한다.

boundary.wgsl

src/lab/engine/shaders/boundary.wgsl

// Stage 5 · boundaries: edge faces (walls, inflow, outflow, periodic copies) and faces touching obstacles.
// Tangential wall conditions are applied when sampling (ghost values in wgsl.ts), not stored here.
override SET_OUTFLOW: bool = true; // false after projection: outflow faces keep their projected values

@group(0) @binding(1) var<storage, read_write> uB: array<f32>;
@group(0) @binding(2) var<storage, read_write> vB: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;
@group(0) @binding(4) var<storage, read> solidVel: array<vec2f>;

fn solidVelAt(i: i32, j: i32) -> vec2f {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  return solidVel[iC(u32(ii), u32(jj))];
}

fn uTouches(i: i32, j: i32) -> bool { return isSolid(i - 1, j) || isSolid(i, j); }
fn vTouches(i: i32, j: i32) -> bool { return isSolid(i, j - 1) || isSolid(i, j); }

// Velocity of a face touched by an obstacle: mean of the adjacent solid cells.
fn uSolid(i: i32, j: i32) -> f32 {
  let a = isSolid(i - 1, j);
  let b = isSolid(i, j);
  if (a && b) { return 0.5 * (solidVelAt(i - 1, j).x + solidVelAt(i, j).x); }
  if (a) { return solidVelAt(i - 1, j).x; }
  return solidVelAt(i, j).x;
}

fn vSolid(i: i32, j: i32) -> f32 {
  let a = isSolid(i, j - 1);
  let b = isSolid(i, j);
  if (a && b) { return 0.5 * (solidVelAt(i, j - 1).y + solidVelAt(i, j).y); }
  if (a) { return solidVelAt(i, j - 1).y; }
  return solidVelAt(i, j).y;
}

// Final value of an interior face after this pass, computed instead of read when this pass rewrites it.
fn uFinal(i: i32, j: i32) -> f32 {
  if (uTouches(i, j)) { return uSolid(i, j); }
  return uB[iU(u32(i), u32(j))];
}

fn vFinal(i: i32, j: i32) -> f32 {
  if (vTouches(i, j)) { return vSolid(i, j); }
  return vB[iV(u32(i), u32(j))];
}

fn edgeValue(kind: u32, inflowValue: f32, inner: f32, current: f32) -> f32 {
  if (kind == EDGE_INFLOW) { return inflowValue; }
  if (kind == EDGE_OUTFLOW) { return select(current, inner, SET_OUTFLOW); }
  return 0.0; // no-slip and slip walls: no flow through the wall
}

@compute @workgroup_size(8, 8)
fn boundary(@builtin(global_invocation_id) id: vec3u) {
  let i = i32(id.x);
  let j = i32(id.y);
  let nx = i32(P.nx);
  let ny = i32(P.ny);
  if (i <= nx && j < ny) {
    let k = iU(id.x, id.y);
    if (periodicX() && i == nx) {
      uB[k] = uFinal(0, j);
    } else if (uTouches(i, j)) {
      uB[k] = uSolid(i, j); // obstacles win over edge rules, also on edge faces
    } else if (!periodicX() && i == 0) {
      uB[k] = edgeValue(P.edgeL, P.inflowL, uFinal(1, j), uB[k]);
    } else if (!periodicX() && i == nx) {
      uB[k] = edgeValue(P.edgeR, P.inflowR, uFinal(nx - 1, j), uB[k]);
    }
  }
  if (i < nx && j <= ny) {
    let k = iV(id.x, id.y);
    if (periodicY() && j == ny) {
      vB[k] = vFinal(i, 0);
    } else if (vTouches(i, j)) {
      vB[k] = vSolid(i, j);
    } else if (!periodicY() && j == 0) {
      vB[k] = edgeValue(P.edgeB, P.inflowB, vFinal(i, 1), vB[k]);
    } else if (!periodicY() && j == ny) {
      vB[k] = edgeValue(P.edgeT, P.inflowT, vFinal(i, ny - 1), vB[k]);
    }
  }
}

셰이더 앞에 붙는 공통 코드

엔진은 셰이더 파일마다 앞에 공통 코드를 붙여서 돌린다. 대부분의 셰이더는 COMMON(격자의 크기, 한 스텝의 시간, 점성, 가장자리의 종류, 체적력 같은 값을 담은 구조체와, 면과 칸의 번호를 셈하는 함수)으로 시작한다. 고체 칸을 알아보는 SOLID와, ② 이류가 흐름을 거슬러 짚는 TRACE가 필요한 셰이더에 붙는다. 다중 격자법의 셰이더는 COMMON 대신 한 층의 격자 크기를 담은 LEVEL_WGSL로 시작한다. 그 밖에 아무 자리의 속도를 둘레 면의 값으로 어림하는 함수, 잉크를 어림하는 함수, 한 층의 압력 식은 셰이더마다 쓰는 값의 이름에 맞춰 wgsl.ts의 함수가 써 넣는다.

COMMON · SOLID · TRACE · LEVEL_WGSL

src/lab/engine/wgsl.ts · COMMON

struct Params {
  nx: u32, ny: u32, h: f32, dt: f32,
  nu: f32, omega: f32, edgeL: u32, edgeR: u32,
  edgeB: u32, edgeT: u32, inflowL: f32, inflowR: f32,
  inflowB: f32, inflowT: f32, gx: f32, gy: f32,
  splatCount: u32, pad0: u32, pad1: u32, pad2: u32,
};
@group(0) @binding(0) var<uniform> P: Params;

const EDGE_NOSLIP: u32 = 0u;
const EDGE_SLIP: u32 = 1u;
const EDGE_INFLOW: u32 = 2u;
const EDGE_OUTFLOW: u32 = 3u;
const EDGE_PERIODIC: u32 = 4u;

fn periodicX() -> bool { return P.edgeL == EDGE_PERIODIC; }
fn periodicY() -> bool { return P.edgeB == EDGE_PERIODIC; }
fn wrapi(i: i32, n: i32) -> i32 { return ((i % n) + n) % n; }

// u face (i, j): x = i*h, y = (j+0.5)*h | v face (i, j): x = (i+0.5)*h, y = j*h | cell (i, j): centre
fn iU(i: u32, j: u32) -> u32 { return j * (P.nx + 1u) + i; }
fn iV(i: u32, j: u32) -> u32 { return j * P.nx + i; }
fn iC(i: u32, j: u32) -> u32 { return j * P.nx + i; }

fn domainSize() -> vec2f { return vec2f(f32(P.nx), f32(P.ny)) * P.h; }

// Keeps a traced point inside the domain (wraps along periodic axes).
fn confine(p: vec2f) -> vec2f {
  let L = domainSize();
  var q = p;
  if (periodicX()) { q.x = q.x - floor(q.x / L.x) * L.x; } else { q.x = clamp(q.x, 0.0, L.x); }
  if (periodicY()) { q.y = q.y - floor(q.y / L.y) * L.y; } else { q.y = clamp(q.y, 0.0, L.y); }
  return q;
}

src/lab/engine/wgsl.ts · SOLID

fn isSolid(i: i32, j: i32) -> bool {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  if (ii < 0 || jj < 0 || ii >= i32(P.nx) || jj >= i32(P.ny)) { return false; }
  return solid[iC(u32(ii), u32(jj))] != 0u;
}

// A u face (i, j) lies between cells (i-1, j) and (i, j). Fixed faces are owned by the boundary kernel:
// edge faces, the duplicate face nx on a periodic axis, and every face touching an obstacle.
fn uFixed(i: u32, j: u32) -> bool {
  if (periodicX()) {
    if (i == P.nx) { return true; }
  } else if (i == 0u || i == P.nx) {
    return true;
  }
  return isSolid(i32(i) - 1, i32(j)) || isSolid(i32(i), i32(j));
}

fn vFixed(i: u32, j: u32) -> bool {
  if (periodicY()) {
    if (j == P.ny) { return true; }
  } else if (j == 0u || j == P.ny) {
    return true;
  }
  return isSolid(i32(i), i32(j) - 1) || isSolid(i32(i), i32(j));
}

src/lab/engine/wgsl.ts · TRACE

// dtSigned > 0 traces back in time (semi-Lagrangian departure point); dtSigned < 0 traces forward.
fn trace(p: vec2f, dtSigned: f32) -> vec2f {
  let mid = confine(p - 0.5 * dtSigned * velAt(p));
  return confine(p - dtSigned * velAt(mid));
}

src/lab/engine/wgsl.ts · LEVEL_WGSL

struct Level {
  nx: u32, ny: u32, cnx: u32, cny: u32,
  h: f32, fScale: f32, edgeL: u32, edgeR: u32,
  edgeB: u32, edgeT: u32, omega: f32, iterations: u32,
};
@group(0) @binding(0) var<uniform> L: Level;

const LEVEL_OUTFLOW: u32 = 3u;
const LEVEL_PERIODIC: u32 = 4u;

다시, 식 전체

이 교재가 한 장에 하나씩 쌓아 온 식이다.

\mat{\frac{D\mathbf{u}}{Dt}} = \frac{\partial \mathbf{u}}{\partial t} + \adv{(\mathbf{u}\cdot\nabla)\mathbf{u}} = \prs{-\frac{\nabla p}{\rho}} + \vis{\nu\nabla^2\mathbf{u}} + \frc{\mathbf{f}}, \qquad \cont{\nabla\cdot\mathbf{u} = 0}

$\vel{\mathbf{u}}$ 는 흐름의 속도다. 물 덩어리를 따라가며 볼 수도, 한 점에서 지켜볼 수도 있다(1장).
$\cont{\nabla\cdot\mathbf{u} = 0}$ 은 들어온 만큼 나가야 한다는 질량 보존이다. 압력을 끄면 물이 몰리거나 빈다(2장).
$\prs{-\nabla p/\rho}$ 는 물이 몰리거나 비지 않도록 어디서나 한꺼번에 맞춰져 물을 미는 압력이다(3장). 날개를 떠받치는 것도 이 힘이다.
$\vis{\nu\nabla^2\mathbf{u}}$ 는 이웃과 속도를 고르게 나누는 점성이다. 저어 만든 소용돌이를 잦아들게 한다(4장).
$\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}}$ 는 흐름이 제 속도를 실어 나르는 일이다. 끄면 소용돌이가 그 자리에 선다(5장).
$\frc{\mathbf{f}}$ 는 바깥에서 주는 힘이다. 끄면 저어도 물이 움직이지 않는다(이 장).

한 문장으로

나비에-스토크스 방정식의 항은 저마다 한 가지 일을 맡고, 엔진은 스텝마다 그 일을 단계로 나누어 차례로 돌며, 날개의 양력처럼 모든 항이 함께 만드는 흐름도 그 항들의 일로 읽을 수 있다.

생각해 볼 질문

종합 실험실에서 압력 스위치만 끄면 원기둥 뒤가 비고 물이 한쪽에 몰린다. 압력 항이 켜져 있을 때는 무엇이 이 일을 막을까?

답 보기

압력이다. 압력은 물이 몰리거나 비지 않도록 어디서나 한꺼번에 맞춰진다(3장). 원기둥 앞에서는 압력이 높아져 다가오는 물을 미리 옆으로 비켜 보내고, 원기둥 뒤에서는 압력이 낮아져 둘레의 물을 끌어들인다. 엔진에서는 ④ 압력 투영 단계가 발산을 재고, 그 발산을 지울 압력을 풀어, 그 기울기만큼 속도를 뺀다. 이 장의 실험실에서 압력을 다시 켜자 발산이 0.000005 안으로 돌아왔다.

생각해 볼 질문

이 장의 날개는 위아래가 대칭이다. 받음각을 −8도로 하여 앞전을 아래로 숙이면, 양력은 어느 쪽을 향하고, 앞에서 같이 출발한 띠는 어느 쪽이 먼저 뒷전에 닿을까?

답 보기

대칭 날개를 −8도로 숙인 흐름은 8도의 흐름을 위아래로 뒤집은 것이다. 그래서 날개 아래의 흐름이 빠르고 압력이 낮아져, 양력은 아래를 향한다. 양력 계수는 8도의 0.49와 크기가 같고 부호가 반대다. 띠는 아래로 지나간 물이 먼저, 0.13초쯤 앞서 뒷전에 닿고, 날개 뒤의 흐름은 위로 꺾인다. 날개의 모양은 그대로인데 양력의 방향이 바뀌었다. 양력을 정하는 것은 표면의 길이가 아니라 날개가 흐름을 꺾는 방향이다.