【游戏流体力学基础及Unity代码（七）】车流量问题，非线性水波以及burgers方程
封面是水波的干涉效果。建议阅读本章内容前先去看看一个视频《如何解决交通堵塞》，油管或是B站翻译版车流量问题
你想想，你带着朋友，坐着汽车，吃着火锅唱着歌，突然就堵车了。为何堵车如此难以解决呢？
首先，把公路简化成一条单行道，并且划分成许多区域，假设车辆会从左边的区域驶往右边的区域。
假设每个网格中汽车的速度的u_i范围是0~1，汽车的密度rho_i范围是0~128。那么u_i*rho_i就是在下个时刻从网格i进入网格i+1的数量。并且假设速度恒定，写成离散公式就是如下，其实就是一维对流方程：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Crho%5E%7Bt%2B1%7D_i+-+%5Crho%5E%7Bt%7D_i%7D%7B%5CDelta+t%7D+%3D+c%5Cfrac%7B%5Crho%5Et_%7Bi-1%7D-%5Crho%5E%7Bt%7D_%7Bi%7D%7D%7B%5CDelta+x%7D+%5C%5C 写成python代码如下：
import numpy as np
length = 8
iterations = 10
h = np.zeros((iterations,length))
h = h = 64
for k in range(0,iterations-1):
for i in range(1,length):
h = h + (h - h) * 0.5用Spyder这个超级漂亮的IDE查看车辆密度变换如下：
可以看到车流在缓缓向右移动，并且单个网格的车辆密度没原来那么高了。
但是路要一步步走，汽车要是开得太快了，会撞到别人。周围车辆少的时候，可以开快一点，周围车辆多的时候，就要慢一点。这时候速度不再是恒定的，而是与密度有关系，比如下式：
https://www.zhihu.com/equation?tex=u+%3D+u_%7Bmax%7D%281+-+%5Cfrac%7B%5Crho%7D%7B%5Crho_%7Bmax%7D%7D%29+%5C%5C umax就是1，rho_max在这里就是128，于是原来的方程变成了下面这样非线性的，因为速度是密度的函数：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%5Crho%7D%7B%5Cpartial+t%7D+%2B+u%28%5Crho%29%5Cfrac%7B%5Cpartial+%5Crho%7D%7B%5Cpartial+x%7D+%3D+0%5Ctag%7B1%7D%5C%5C
写成代码就是
import numpy as np
length = 8
iterations = 10
h = np.zeros((iterations,length))
u = np.zeros((iterations,length))
h = h = 64
umax = 0.5
hmax = 128.0
for k in range(0,iterations-1):
for i in range(0,length-1):
u = (1 - h / hmax)* umax
for i in range(1,length):
h = h + h * u - h * u最终结果如下：
把1式变成一般形式，就是经典的无粘性Burgers方程，这是一个非线性方程。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+t%7D+%2B+u%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%3D+0+%5C%5C 加上一个粘性项，就是粘性Burger方程。
https://www.zhihu.com/equation?tex=%7B%5Cpartial+u+%5Cover+%5Cpartial+t%7D+%2B+u+%7B%5Cpartial+u+%5Cover+%5Cpartial+x%7D+%3D+%5Cnu+%7B%5Cpartial%5E2+u+%5Cover+%5Cpartial+x%5E2%7D+%5C%5C 不仅是流体，隔壁量子力学，光学要研究的物质许多都有非线性特征。比如薛定谔方程(Schrdinger equation)：
https://www.zhihu.com/equation?tex=-%5Cfrac%7B%5Chbar%5E2%7D%7B2m%7D+%5Cfrac%7Bd%5E2+%5Cpsi%7D%7Bdx%5E2%7D+%2B+V%5Cpsi+%3D+E%5Cpsi%5C%5C 2016年就有篇用薛定谔方程模拟烟雾的siggraph《Schrodinger's Smoke》。不过现在可以把这三个公式忘了，我们继续研究如何开车。
红灯与绿灯
当前方十字路口出现红灯时，前方车辆就会停下来，后面的车辆的也会跟着停下来。
这时候，堵的地方，也就说密度最大的地方，和不堵的地方，也就说密度还没有达到最大的地方，两者之间的界限处密度会急剧变化。出现这种物理量急剧变化的地方，就叫它激波(Shock Wave)，如下图。超音速飞机会产生音爆也是因为这种激波哦。
如上图第零秒。车辆的速度向右，激波的右边车辆密度达到最大，而激波的左边密度没有达到最大。研究激波形成后车辆如何变化的问题就是黎曼问题(Riemann Problem)。关于这个问题可在各种有关Burgers方程的介绍中看到，看看《Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction》 by Eleuterio F. Toro 也是个不错的选择。
这里我们只讨论最简单的情况，激波的速度应该是多少？比如激波右边车辆密度最大，为3，激波左边为1。那么结果应该如上图。令f1代表激波左的速度，f2是激波右边的速度。rho_1是左边密度，rho_2是右边密度。当激波左边速度为2时，结果就是下式，对应上图中间的情况
https://www.zhihu.com/equation?tex=s+%3D+%5Cfrac%7Bf%28%5Crho%29_1+-+f%28%5Crho%29_2%7D%7B%5Crho_1+-+%5Crho_2%7D+%3D+%5Cfrac%7B2+-+0%7D%7B1+-+3%7D+%3D+-1+%5C%5C 上面这个公式也可以像下面这样推导。假设一段连续的区间，激波的位置是x0
这个区间上总的速度流量，等于激波左边的速度总量，加上激波右边的速度总量
https://www.zhihu.com/equation?tex=%5Cint_%7Bx1%7D%5E%7Bx2%7Du%28x%2Ct%29dx+%3D+%28x_0+-+x_1%29u_L+%2B+%28x_2+-+x_0%29u_R+%5C%5C 对时间微分，就得到了激波位置的变化量，也就是激波的速度s
https://www.zhihu.com/equation?tex=%5Cfrac%7Bd%7D%7Bdt%7D%5Cint_%7Bx1%7D%5E%7Bx2%7Du%28x%2Ct%29dx+%3D+s%28u_L+-+u_R%29%5C%5C 把上面的式子写成激波速度的函数，这就是Rankine-Hugoriot condition。f在这里就是速度流量。
https://www.zhihu.com/equation?tex=s+%3D+%5Cfrac%7Bf%28u%29_1+-+f%28u%29_2%7D%7Bu_1+-+u_2%7D+%5C%5C 绿灯时，前方车辆依次继续前进，密度从左到右逐渐降低，并且没有那么陡峭，形成了膨胀波(expansion wave)。
关于激波和膨胀波同样推荐看王洪伟老师的讲解https://www.bilibili.com/video/BV1pW411X7xE?p=14，通俗易懂。这个问题在流体力学中，就是大坝破坏(dam break)问题。左侧的水和右边的土地原本被大坝隔开，现在大坝突然消失，水会从左边进入右边，这时候流体的运动要比这里的车流量复杂得多，下篇将会讨论。
如何设计道路与规则，才能让通行效率最高，这就是车流量问题(Traffic Flow)。
你给我翻译翻译什么叫非线性
线性的挺简单。比如一个水波是sin(x)，另一个水波也是sin(x)，那么它们如果遇上了，线性叠加后就变成了2sin(x)。然后真实情况往往是这样的：
动图截取自https://www.youtube.com/watch?v=G-xHbcMm_Gs。两个非常平缓的水波相遇后却得到一个非常陡峭的水波。又比如下图，截取自https://www.youtube.com/watch?v=Iuv6hY6zsd0&t=346s。用球不断振动水面，就可以创造出两个水波。
水波单独看起来没什么神奇的，当个水波相交的时候却像是下面这样，相交处的水波并不是原来两个波简单的线性叠加，由于destructive interference
最后变成这样的好看到可以做壁纸的东西：
巧了，隔壁光学中的杨氏双缝干涉(Young`s double slit interference)的结果看起来和这个居然一模一样，图我就不放了。看来这种非线性模式很值得研究。其实海岸边形成的破碎的浪也是因为这种线性。在深水中，水波并不高。而进入浅水区后，海岸的阻力开始影响水波的底部，这样一来水波的底部移动得慢，而水波上部移动得快，水波的前部移动得慢，而水波后边移动得快，于是水波逐渐变得又窄又高，最后水波上部前倾，最终坍塌，水浪就这样被打在沙滩上。
色散
除了非线性外，水波还有一个重要的特征。例如向水中丢一块石头，水面会产生一系列波，如下。
波长的波传播的快，比如2处。波长短的波传播得慢，比如1处。虽然一开始长短波都聚集在一起，但随着时间的推移，它们会逐渐分开。比如下图，截取自https://www.youtube.com/watch?v=EIqKG5TiSYs，最终的结果是上方的红色波，红色波在某个时刻形成了很高的波，但最后开始分散成波长不同的波。这就是色散(dispersion)。
巧了，隔壁光学也是这样玩的。阳光近似白色，实际上是由各种不同波长的光组成的。下雨之后出现阳光时，还在空气中的水滴会折射阳光。阳光进入雨滴后，色散导致不同波长的光传播的方向不一样，于是就出现了彩虹。
孤立波
很多介绍非线性波的书上都会记载英国科学家John Scott Russell在1884年观察的奇特现象，他在行驶的船的前面看到了这样的水波，水波以很快的速度离开船头，而且形状几乎没什么改变。Russell想着“别急，让水波继续前进一会儿”，于是骑马沿着河道跟踪了两公里后，这个水波才最终消失。这个水波看起来像是这样的：
这种波又叫孤立波(solitary wave)，或者叫做孤立子(soliton)。这种现象广泛存在于水波中，这是由于水波的非线性和色散共同造成的。孤立波的第一个特点是水波可以保持前进而形状基本不变。第二个特点是速度取决于水波的形状，以及水面的深度。第三个特点是两个孤立波在相遇后仍然能保持原来的形状，如下，截取自https://github.com/W-J-Trenberth/Finite-difference-for-KdV
之后1895年，荷兰数学家D.Korteweg和他的学生G. de Veries求得了单向运动的孤立波运动方程，即Kdv方程。一般形式如下，也就是非线性和色散的结合：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+t%7D+%2B+6u%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%2B+%5Cfrac%7B%5Cpartial%5E3+u%7D%7B%5Cpartial+x%5E3%7D%3D+0+%5C%5C 使用Lax-Wendroff Scheme写成离散形式如下。很多计算流体力学的书上都有Lax-Wendroff以及其它方法的解释，本系列不会再说明了。
https://www.zhihu.com/equation?tex=%5Cbegin%7Bequation%7D+%5Cbegin%7Bsplit%7D+U_j%5E%7Bn%2B1%7D+%3D+U_j%5E%7Bn%7D+%26-+3%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D+U_j%5En%5CBig%5B+%5Cleft%28U_%7Bj%2B1%7D%5En+%2B+U_%7Bj%7D%5En%5Cright%29+-+3%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D+%5Cbig%28+%5Cleft%28+U_%7Bj%2B1%7D%5En+%5Cright%29%5E2+%2B+%5Cleft%28+U_%7Bj%7D%5En+%5Cright%29%5E2++%5Cbig%29%5C%5C+%26+-+%5Cleft%28U_%7Bj%7D%5En+%2B+U_%7Bj-1%7D%5En%5Cright%29+%2B+3%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D++%5Cbig%28+%5Cleft%28+U_%7Bj%7D%5En+%5Cright%29%5E2+-+%5Cleft%28+U_%7Bj-1%7D%5En+%5Cright%29%5E2++%5Cbig%29%5CBig%5D+%5C%5C+%26+-+%5Cfrac%7B%5CDelta+t%7D%7B+2%5CDelta+x%5E3%7D+%5Cleft%28U_%7Bj%2B2%7D%5En+-+2U_%7Bj%2B1%7D%5En+%2B+2U_%7Bj-1%7D%5En+-+U_%7Bj-2%7D%5En+%5Cright%29+%5Cend%7Bsplit%7D+%5Cend%7Bequation%7D
写成python代码如下：
import numpy as np
L = 20 # 区域长度
grid_points = 101#网格数量
delta_x = L / (grid_points-1)
delta_t = 0.000004
num_tsteps = int((1 - 0) / delta_t)
u = np.zeros((num_tsteps,grid_points))#初始化水波
x = np.linspace(-L/4, 3*L/4, grid_points)
for i in range(0,grid_points):
u = 0.5 * 3.0 *(1 / np.cosh((np.sqrt(3.0)/2)*x))**2
for k in range(0,num_tsteps-1):#Lax-Wendroff Scheme
Ui = u
Ui_p1 = np.roll(u, -1)
Ui_p2 = np.roll(u, -2)
Ui_m1 = np.roll(u, 1)
Ui_m2 = np.roll(u, 2)
u = Ui - 3 * (delta_t / delta_x) * Ui * ((Ui_p1 - Ui_m1) - \
3 * (delta_t / delta_x) * (Ui_p1**2 + Ui_m1**2))\
- delta_t / (2 * delta_x**3) * (Ui_p2 - Ui_m2 - 2.0 * (Ui_p1 - Ui_m1))
#这个变量仅仅是为了方便查看，每隔1000个dt采样一次，否则原来的变化太慢了
us = np.zeros((250,grid_points))
for i in range(0,250):
us = u用spyder查看us变量，横轴就是x轴，纵轴是时间。一开始最水波高峰出现第25个网格为1.5，
但是17 * 1000次迭代后，水波最高峰出现在了第26个网格，大约是1.5
此后也是最高峰将右移，但水波形状基本上没有改变。也就是说非线性现象和色散现象达到了一种平衡，没有让水波聚集或扩散开来，而是一直维持着原来的形状。如下图
这份代码修改自https://github.com/TiagoJCor/Numerical-Methods-KDV，也可以用nbviewer查看notebookhttps://nbviewer.jupyter.org/github/TiagoJCor/Numerical-Methods-KDV/blob/master/numerical_KDV.ipynb。
但船在行驶的时候，不仅前方会生成孤立波，两边也会以及后方也会产生各种水波。这就是船行波(bow wave)。
更多有关船行波的介绍可看https://www.youtube.com/watch?v=D14QuUL8x60以及《3D Modeling of Nonlinear Wave Phenomena on Shallow Water Surfaces》的第2.6节https://b-ok.global/book/3492383/24b31d。
但是想要模拟出真实的浅水效果，kdv方程还远远不够。不过弄清楚本篇的概念后，下一篇就可以用浅水波方程(shallow water equation)模拟真正的雨落池塘，用unity复现evan在2011年在webgl上的水面模拟http://madebyevan.com/webgl-water/，以及洪水冲破堤坝后的流动。敬请期待！
上一篇：光影帽子：【游戏流体力学基础及Unity代码（六）】用NavierStokes方程模拟粘性染料流动
推荐阅读
https://github.com/clawpack/riemann_book，notebook很不错，网页版在这里http://www.clawpack.org/riemann_book/html/Index.html
《Free-surface flow shallow-water dynamics》by Nikolaos D. Katopodes。图很多，对各种概念解释得很清楚，非常推荐
《3D Modeling of Nonlinear Wave Phenomena on Shallow Water Surfaces》by Iftikhar B. Abbasov 同样对公式解释得很不错，非常推荐
《Nonlinear Dispersive Waves Asymptotic Analysis and Solitons》by MARK J. ABLOWITZ.
《Solitary Waves in Dispersive Complex Media Theory, Simulation, Applications》 by Vasily Y. Belashov, Sergey V. Vladimirov
《Partial Differential Equations and Solitary Waves Theory》by Abdul-Majid Wazwaz
页:
[1]