在计算机图形学中，模拟真实世界的流体动力学是一个复杂但非常有趣的挑战。本文将探讨如何通过数学模型和编程技术实现交互式水面效果。将从流体动力学的基本原理出发，逐步深入到具体的编程实现。
首先，需要理解流体的一些基本假设：
基于这些假设，可以使用简化的流体力学方程。在二维情况下，有以下方程：
du(x,t)/dt + u(x,t).du(x,t)/dx + g.dh(x,t)/dx = 0
dp(x,t)/dt + d(u(x,t)p(x,t))/dx = d(h(x,t)-b(x))/dt + p(x,t).du(x,t)/dx + u(x,t).dp(x,t)/dx = 0
其中，g 是重力加速度，h(x, t) 是液体表面的高，b(x) 是容器底部的高度，p(x, t) = h(x, t) - b(x) 是液体的深度，u(x, t) 是垂直柱体的液体水平速度。第一个方程表示牛顿定律（F=ma），描述了液体的运动；第二个方程表示体积守恒。
这两个方程是非线性的，但如果液体速度较小且深度变化缓慢，可以进行以下近似：
du(x,t)/dt + g.dh(x,t)/dx = 0
dh(x,t)/dt + p(x).du(x,t)/dx = 0
通过对第一个方程关于x求导，第二个方程关于t求导，可以得到：
d^2u(x,t)/dxdt + g.d^2h(x,t)/dx^2 = 0
d^2h(x,t)/dt^2 + p(x).d^2u(x,t)/dt^2 = 0
由于u是一个合适的函数（其二阶导数是连续的），其二阶偏导数相等，可以得到一个不包含u的方程：
d^2h(x,t)/dt^2 = gp(x)*d^2h(x,t)/dx^2
这个方程是一个波动方程，表示高度随时间和x坐标的变化。在三维情况下，方程的形式为：
d^2h(x,y,t)/dt^2 = gp(x,y) * (d^2h(x,y,t)/dx^2 + d^2h(x,y,t)/dy^2)
为了能够使用高度网格，需要一个离散公式。此外，这个方程是非线性的，因为p(x, y)的存在。为了简化，假设p是常数，即无论深度如何，波速都是恒定的。得到离散化的方程：
d^2h(x,y,t)/dt^2 = gp(d^2h(x,y,t)/dx^2 + d^2h(x,y,t)/dy^2)
可以使用中心差分法来离散化这个方程的项：
d^2h(x,y,t)/dt^2 => (Dh(x,y,t+1) - Dh(x,y,t))/Dt^2
d^2h(x,y,t)/dx^2 => (Dh(x+1,y,t) - Dh(x,y,t))/Dx^2
d^2h(x,y,t)/dy^2 => (Dh(x,y+1,t) - Dh(x,y,t))/Dy^2
通过设置Dr = Dx = Dy = 网格分辨率，得到方程的离散形式：
(h(x,y,t+1) + h(x,y,t-1) - 2h(x,y,t))/Dt^2 = gp/Dr^2 . (h(x+1,y,t)+h(x-1,y,t)+h(x,y+1,t)+h(x,y-1,t)-4h(x,y,t))
可以使用更紧凑的符号表示这个递推关系：
[0 1 0] * h(x,y,t+1) = gpDt^2/Dr^2 * [1 -4 1] * h(t) - h(x,y,t-1) + 2h(x,y,t)
然后可以得到：
h(x,y,t+1) = gpDt^2/Dr^2 * [1 0 1] * h(t) - h(x,y,t-1) + (2 - 4gpDt^2/Dr^2) * h(x,y,t)
通过设置gpDt^2/Dr^2 = 1/2，可以简化这个关系，得到：
h(x,y,t+1) = 1/2 * (h(x+1,y,t) + h(x-1,y,t) + h(x,y+1,t) + h(x,y-1,t) - h(x,y,t-1))
这个递推关系有一个2的步长：t+1时刻的高度与t和t-1时刻的高度有关。可以使用两个高度矩阵H1和H2来实现这一点：
H2[x, y] = 1/2(H1[x+1, y] + H1[x-1, y] + H1[x, y+1] + H1[x, y-1]) - H2[x, y]
然后可以在每一步交换这两个矩阵。计算一个点的新高度只需要4次加法，1次减法和1次右移（用于除以2）。
从得到的结果中，减去1/2n的结果（即这个结果右移n位），以实现一定的阻尼（n=4给出相当美观的结果，n < 4给出更多的粘性，n > 4给出更多的流动性）。
需要注意的是，这些点的高度是有符号的，0是静止的中性水平。从高度矩阵中，可以很容易地构建一个多边形表示，通过考虑每个网格盒作为一个四边形（或2个三角形），其4个顶点的高度由相邻4个网格盒的高度给出。
在示例中，使用GL_TRIANGLE_STRIP进行网格化，并使用一些多通道效果来实现现实感。
首先，根据顶点法线按比例扰动第一组纹理坐标（标志的纹理坐标）。这看起来像是折射。
/*
计算归一化
*/
for (x = 0; x < FLOTSIZE*2; x++) {
for (y = 0; y < FLOTSIZE*2; y++) {
sqroot = sqrt(_snormal[x][y][0]*_snormal[x][y][0] +
_snormal[x][y][1]*_snormal[x][y][1] + 0.0016f);
_snormaln[x][y][0] = _snormal[x][y][0]/sqroot;
_snormaln[x][y][1] = _snormal[x][y][1]/sqroot;
_snormaln[x][y][2] = 0.04f/sqroot;
//
perturbate coordinates of background
//
mapping with the components X,Y of normals...
//
simulate refraction
_newuvmap[x][y][0] = _uvmap[x][y][0] + 0.05f*_snormaln[x][y][0];
_newuvmap[x][y][1] = _uvmap[x][y][1] + 0.05f*_snormaln[x][y][1];
}
}
然后使用一个假的环境映射公式计算第二组纹理坐标（天空的纹理坐标）。
//
really simple version of a fake envmap generator
for (x = 0; x < FLOTSIZE; x++) {
for (y = 0; y < FLOTSIZE; y++) {
//
trick : xy projection of normals ->
//
assume reflection in direction of the normals
//
looks ok for non-plane surfaces
_envmap[x][y][0] = 0.5f + _snormaln[x][y][0]*0.45f;
_envmap[x][y][1] = 0.5f + _snormaln[x][y][1]*0.45f;
}
}
然后使用多纹理和混合将纹理混合在一起。
glEnable(GL_BLEND);
//
use texture alpha-channel for blending
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glActiveTextureARB(GL_TEXTURE0_ARB);
glBindTexture(GL_TEXTURE_2D, 2);
//
2nd texture -> background ..
glActiveTextureARB(GL_TEXTURE1_ARB);
glBindTexture(GL_TEXTURE_2D, 1);
//
2nd texture -> envmap
glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
/*
enable texture mapping and specify ourself texcoords */
glDisable(GL_TEXTURE_GEN_S);
glDisable(GL_TEXTURE_GEN_T);
for (x = 0; x < strip_width; x++) {
glBegin(GL_TRIANGLE_STRIP);
//
WARNING: glTexCoord2fv BEFORE glVertex !!!
glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x+1][1]);
glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x+1][1]);
glVertex3fv(_sommet[x+1][1]);
//
otherwise everything is scrolled !!!
for (y = 1; y < strip_width; y++) {
glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x][y]);
glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x][y]);
glVertex3fv(_sommet[x][y]);
glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x+1][y+1]);
glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x+1][y+1]);
glVertex3fv(_sommet[x+1][y+1]);
}
glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x][y]);
glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x][y]);
glVertex3fv(_sommet[x][y]);
glEnd();
}