如果流体没有粘性,也就是理想流体,我们就可以借助流函数求解速度场,再通过流函数的导数运算计算速度场。流函数满足拉普拉斯方程,谢天谢地,幸好是拉普拉斯方程,求解简单的很,完全与前面的温度场扩散方程一样了:二维流动的计算域网格中,如网格均匀,则中间节点的流函数为四个直接相邻节点流函数的平均值,这与温度场完全一致。如果网格不均匀,稍微要复杂一点点。
其中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写了一个理想流体流动的例子,流动和残差结果如下:
请读者自行完成程序开发(正文完)
本文分享自 传输过程数值模拟学习笔记 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!