好奇韦神求解的纳维斯托克斯方程长什么样？

用微信扫码二维码

分享至好友和朋友圈

北大韦东奕用数学一天搞定六个博士四个月都没搞定的物理问题，一个科技公司用PS5做了一个cluster，用来模拟某些产品的物理性能。但是面临一个问题，“模型越复杂，模拟失真越高”。他们用了六个博士，搞了四个月也没有完全解决。总是出现模拟和实际测试有很大区别的问题。他们总结主要是他们对纳维斯托克斯方程的处理有问题，但是又搞不清楚哪里有问题。

五一劳动节前，科研团队中有一个人说他同学在北大，可以联系一下韦东奕，韦神，让他帮忙看看。结果一天时间韦神就把全部方程发过来了。科研团队从下午调试到凌晨。早上六点多的时候，他们发现实验成功了。和过往的真实数据匹配99.8%。韦东奕太神了！科技公司想给韦东奕钱，韦神不要。说太简单了，没必要要钱。最后说了半天，给韦神的公交卡充了钱。

好奇韦神求解的纳维斯托克斯方程长什么样？

纳维（Navier）-斯托克斯（Stokes）方程（简称N-S方程）

看起来，是不是感觉超级复杂？接下来，请跟随着模小数一起使用MATLAB软件编程，进行数学建模与仿真模拟，数值求解纳维（Navier）-斯托克斯（Stokes）方程！

对于纳维（Navier）-斯托克斯（Stokes）方程，我们首先将区间[0，1]，划分成101个格点，前后再多加2个格点，为了实现outlet边界条件。首先是设置变量，MATLAB源程序代码如下：

clear

clc

%Set the parameters

xmin = 0;

xmax = 1;

N = 100;

dx = (xmax-xmin)/N;

tmin = 0;

tmax = 0.5;

dt = 0.001;

v = 1;

%Set the initial condition

x = linspace(xmin-dx, xmax+dx, N+3);

u0 = [ones(1,5), zeros(1,103-5)];

utemp = u0;

uplus = u0;

nstep = tmax/dt;

t = 0;

然后编写主循环的思路是：对于时间步数循环，计算边界点的值，根据差分格式，更新每个点的值，时间和空间+1。MATLAB源程序代码如下：

forn = 1:nstep

%calculate the boundary condition

utemp(1) = utemp(3);

utemp(N+3) = utemp(N+1);

%calculate the new value

fori = 2:N+2

uplus(i)= utemp(i) - v*dt*(utemp(i+1) - utemp(i-1))/(2*dx);

%YOU ONLY NEED TO CHANGE THIS LINE!!!

end

%update time step and u

t = t + dt;

utemp = uplus;

%visualization

plot(x,utemp,'bo-','markerfacecolor','b');

shg;

pause(dt);

end

运行结果：

（可以插入附件中的动画或者视频）

所以，早些年，大量研究在寻找好的格式，和格式稳定性理论，其中著名的格式：Lax-Friedrichs(LF) method

MATLAB源程序代码如下：

uplus(i)= (utemp(i+1) + utemp(i-1))/2 - v*dt*(utemp(i+1) -utemp(i-1))/(2*dx);

运行结果

还有一种格式Upwindmethod，格式：

该方法稳定且扩散性较小，但仍然只有一阶准确。可以使用高阶空间插值方案进行改进。这是大名鼎鼎的Godnouv的前身。

MATLAB源程序代码如下：

%This is a demo of First Order Upwind for Advection Equation

clear

clc

%Set the parameters

xmin = 0;

xmax = 1;

N = 100;

dx = (xmax-xmin)/N;

tmin = 0;

tmax = 0.5;

dt = 0.009;

v = 1;

%Set the initial condition

x = linspace(xmin-dx, xmax+dx, N+3);

u0 = exp( -200*(x-0.25).^2 );

%u0 = [ones(1,5), zeros(1,103-5)];

utemp = u0;

uplus = u0;

nstep = tmax/dt;

t = 0;

forn = 1:nstep

%calculate the boundary condition

utemp(1) = utemp(3);

utemp(N+3) = utemp(N+1);

%calculate the new value

fori = 2:N+2

uplus(i)= utemp(i) - v*dt*(utemp(i) - utemp(i-1))/dx;

end

%update time step and u

t = t + dt;

utemp = uplus;

%calculate exact solution

exact = exp(-200*(x-0.25-v*t).^2);

%visualization

plot(x,exact,'r-')

holdon

plot(x,utemp,'bo-','markerfacecolor','b');

holdoff

shg;

pause(dt);

end

Lax-Wendroffmethod：

这个格式的动机，对于时间我们希望更高的精度：

然后呢，因为时间的差分可以近似时间的偏导:

总结：韦神，YYDS！

纳维-斯托克斯方程是流体流动建模的核心。在特定的边界条件（如入口、出口和壁）下求解这些方程，可以预测给定几何体中的流体速度和压力。由于这些方程本身的复杂性，我们只能得到非常有限的解析解。例如，对于两个平行板之间的流动或圆管内的流动，方程的求解会相对容易一些；但对于更为复杂的几何结构，求解方程会非常困难。

戳我戳我

温馨提示：微信公众号信息流改版，每个用户可以设置常读订阅号，这些订阅号将以大卡片的形式展示。因此，如果不想错过“校苑数模”的文章，第一时间获得国赛资讯、答案查询、你一定要进行以下操作：进入“校苑数模”公众号 → 点击右上角的 ··· 菜单 → 选择「设为星标」

特别声明：以上内容(如有图片或视频亦包括在内)为自媒体平台“网易号”用户上传并发布，本平台仅提供信息存储服务。

Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.