6.1 理想流体简单有势流动

如果流体没有粘性,也就是理想流体,我们就可以借助流函数求解速度场,再通过流函数的导数运算计算速度场。流函数满足拉普拉斯方程,谢天谢地,幸好是拉普拉斯方程,求解简单的很,完全与前面的温度场扩散方程一样了:二维流动的计算域网格中,如网格均匀,则中间节点的流函数为四个直接相邻节点流函数的平均值,这与温度场完全一致。如果网格不均匀,稍微要复杂一点点。

其中f(x,y)=0时为拉普拉斯方程,否则为泊松方程,流函数与速度关系如下:

Matlab file exchange上一个计算理想流体流动的例子,使用Matlab计算流体流动,代码如下:

% Copyright 2013 The MathWorks, Inc.
%% Initialization% store input parameters with shorter namesD1 = 0.2; % mD2 = 0.1; % mL1 = 0.5; % mL2 = 0.5; % mU1 = 0.2; % m/snx = 51;  % number of x nodesny = 11;  % number of y nodesnmax = 10000; % maxsimum number of iterations% compute grid lengthdx = (L1+L2)/(nx-1);dy = max(D1,D2)/(ny-1);% nondimensionalizationdX = dx/D1;dY = dy/D1;PSI = ones(ny,nx);%% Boundary conditionsPSI(1,:) = 0; % bottomPSI(floor(D1/dy+1),1:floor(L1/dx+1)) = 1; % left part of topPSI(floor(D2/dy+1),floor(L1/dx+1):nx) = 1; % right part of topfor j = 1:floor(D1/dy+1)    PSI(j,1) = (j-1)*dY; % inletendfor j = 1:floor(D2/dy+1)    PSI(j,nx) = D1/D2*(j-1)*dY; % outletend%% Iteration to solve for stream function PSIerr = 1e-6; % convergence criterian = 0;while n < nmax    n = n + 1;    tempPSI = PSI;    % left part of the channel    for i = 2:floor(L1/dx+1)        for j = 2:(floor(D1/dy+1)-1)            PSI(j,i) = 1/(2/dX^2+2/dY^2)*((tempPSI(j,i+1)+PSI(j,i-1))/dX^2+(tempPSI(j+1,i)+PSI(j-1,i))/dY^2);        end    end    % right part of the channel    for i = floor(L1/dx+1):nx-1        for j = 2:(floor(D2/dy+1)-1)            PSI(j,i) = 1/(2/dX^2+2/dY^2)*((tempPSI(j,i+1)+PSI(j,i-1))/dX^2+(tempPSI(j+1,i)+PSI(j-1,i))/dY^2);        end    end    % checking for convergence    if max(max(abs(PSI-tempPSI))) <= err        break    endend%% Print convergence result to screeenif n < nmax    fprintf('Converged. Total number of iterations: %g\n',n);else    fprintf('NOT Converged! Change n_x, n_y, n_max or convergence criteria.\n')end%% Solve for velocity u and v from stream functionu = zeros(ny,nx);v = zeros(ny,nx);psi = PSI*U1*D1; % dimentionalizationfor j = 2:ny-1    u(j,:) = (psi(j+1,:)-psi(j-1,:))/2/dy;endfor i = 2:nx-1    v(:,i) = -(psi(:,i+1)-psi(:,i-1))/2/dx;end% velocity at boundaries (inviscid, so slip occurs)u(1,:) = u(2,:);u(ny,:) = u(ny-1,:);v(:,1) = v(:,2);v(:,nx) = v(:,nx-1);%% Plot results[X, Y] = meshgrid(0:dx:(L1+L2),0:dy:max(D1,D2));% contour plot of stream functionax1 = subplot(2,1,1);contourf(ax1,X,Y,PSI),colormap(ax1,'jet')% vector plot of velocity u and vax2 = subplot(2,1,2);quiver(ax2,X,Y,u,v),axis([ax1 ax2],[0 L1+L2 0 max(D1,D2)]),% plot the shape of the channel based on input parametersif D1 < D2    cornerposition = [0 D1 L1 D2-D1];elseif D1 > D2    cornerposition = [L1 D2 L2 D1-D2];else    cornerposition = [0 D1 1 1];endrectangle('Parent',ax1,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])rectangle('Parent',ax2,'Position',cornerposition,'FaceColor',[0.94 0.94 0.94])set(get(ax1,'XLabel'),'String','x (m)')set(get(ax1,'YLabel'),'String','y (m)')set(get(ax1,'Title'),'String','Normalized Stream Function \Psi')set(get(ax2,'XLabel'),'String','x (m)')set(get(ax2,'YLabel'),'String','y (m)')set(get(ax2,'Title'),'String','Velocity Profile')

我用javascript写了一个理想流体流动的例子,流动和残差结果如下:

请读者自行完成程序开发(正文完)

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

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

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

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券