6.2 粘性流体涡量-流函数法求解流场

求解N-S方程最大难处就是压力和速度场是耦合的。如何处理压力是核心问题,涡量-流函数方法很好的解决了该难题。将求解速度场和压力场转变为求解涡量和流函数。

而速度场与流函数关系如下:

Matlab file exchange上一个涡量-流函数方法计算流动的例子,使用Matlab计算流体流动,代码如下:

clear;%参数设置Re=10;                        L=1;                          n=100;dh=L/n;psi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0;   for k=1:100000    err=0;    for i=2:n    xi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;    xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;    end         for j=2:n    xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2;    xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2;    end        for i=2:n    for j=2:n        u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);        v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));        err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4-psi(i,j);        psi(i,j)=psi(i,j)+rho*err1;        err2=(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1))/4 ...            -Re*dh*(u(i,j)*(xi(i+1,j)-xi(i-1,j))+v(i,j)*(xi(i,j+1)-xi(i,j-1)))/8-xi(i,j);        xi(i,j)=xi(i,j)+rho*err2;        temp=max(abs(err1),abs(err2));        if err<temp        err=temp;        end    end    end    if err<1e-6    break;endend  kerrrho%psicontour(psi,100);

算法详细请查看张德良《计算流体力学教程》,尤其是要关注到边界及角部的涡量计算设定。

我用javascript写了类似的程序,并用js做了后处理,流场结果如下,Contour就是流函数,而流函数的等值线就是流线:

涡量:

X方向速度:

Y方向速度:

感兴趣的读者可以自行实现程序开发(正文完)

原文发布于微信公众号 - 传输过程数值模拟学习笔记(SongSimStudio)

原文发表时间:2019-04-09

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券