《传热学/流体力学》中几个简单演示程序-Voronoi

周星星9527

发布于 2019-09-03 10:59:53

1.3K0

发布于 2019-09-03 10:59:53

文章被收录于专栏：javascript趣味编程

BGM:

继之前通过测量一些简单算法在微信小程序上运行时间衡量手机综合性能：不服跑个分！-解Laplace偏微分方程测手机/PC性能。此次，又做了大量《传热学》或《流体力学》相关算例补充（热辐射和对流的例子尚没有完成）。

先从吃饭说起，如何就近“觅食”？Voronoi算法为“懒癌”晚期患者就近“觅食”提供了最优化方法。关于Voronoi图或者Delaunay图，之前提过一次，原文在这里。看动画体验下Delaunay三角化：

可以到小程序中体验该交互程序：Delaunay三角化。可见这样得到了一套三角网格单元，类似有限元三角单元网格。事实上Delaunay三角化是网格剖分的一类非常常见的方法，一种二维有限元三角网格剖分思路如下图：

有了网格才能基于该网格离散各类PDE。终于引到本文主题：现将几个《传热学》相关的小程序总结如下，可在微信中点击体验：

《传热学》相关小程序演示动画如下（其中下图1D非稳态导热计算发散，调小时间步长后重新计算，结果收敛！）：

《(计算)流体力学》中的几个小程序，可在微信中点击体验：

Blasius偏微分方程求解速度边界层（理论这里）
理想流体在管道中的有势流动（源码戳这）
涡量-流函数法求解顶驱方腔流动（源码戳这）
SIMPLE算法求解顶驱方腔流动（源码戳这）
Lattice Boltzmann Method计算绕流演示（可能无法在手机上演示参考源码见文尾）

关于《(计算)流体力学》相关的几个小程序演示动画如下：

顺便，《(热工过程)自动控制》中关于PID控制器的仿真可点击此处体验：PID控制演示小程序，（PID控制相关视频见：基础/整定/重要补充）。动画如下：

最后，Lattice Boltzmann Method是研究流体流动热门方向，简单理论入门视频如下：

例如计算绕流的流动如下图：

将其移植到微信小程序，将障碍物右下方10×10=100个节点的速度可视化，在微信小程序开发工具演示效果如下：

可见随着程序迭代，速度场不再杂乱无章。运行一段时间后，速度场如下：

然而在手机微信小程序上查看，该程序并不能正常运行和显示，调试过程头很疼经历如下图，不知道为什么？

点击体验：LBM小程序（很可能无法运行）。这里有一份源于网络的Matlab脚本实现的LBM圆柱绕流源代码：

% =========================================================================
% Channel flow past a cylinderical obstacle, using a LB method
% =========================================================================
% Lattice Boltzmann sample in Matlab
% Original implementaion of Zou/He boundary condition
% =========================================================================
% =========================================================================
clear all 
close all 
clc 
% =========================================================================
 
% GENERAL FLOW CONSTANTS --------------------------------------------------
lx = 400;                                                                  % number of cells in x-direction
ly = 100;                                                                   % number of cells in y-direction
obst_x = lx/5+1;                                                           % position of the cylinder; (exact y-symmetry is avoided)
obst_y = ly/2+3;                         
obst_r = ly/10+1;                                                          % radius of the cylinder
uMax = 0.1;                                                                % maximum velocity of Poiseuille inflow
Re = 100;                                                                  % Reynolds number
nu = uMax * 2.*obst_r / Re;                                                % kinematic viscosity
omega = 1. / (3*nu+1./2.);                                                 % relaxation parameter
maxT = 400000;                                                             % total number of iterations
tPlot = 50;                                                                % cycles
 
% =========================================================================
% D2Q9 LATTICE CONSTANTS --------------------------------------------------
t     = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; 
cx    = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; 
cy    = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; 
opp   = [ 1, 4, 5, 2, 3, 8, 9, 6, 7]; 
col   = [2:(ly-1)]; 
in    = 1;                                                                 % position of inlet
out   = lx;                                                                % position of outlet
[y,x] = meshgrid(1:ly,1:lx);                                               % get coordinate of matrix indices
obst  = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;                        % Location of cylinder
 
obst(:,[1,ly]) = 1;                                                        % Location of top/bottom boundary
bbRegion = find(obst);                                                     % Boolean mask for bounce-back cells
 
 
% =========================================================================
% INITIAL CONDITION: Poiseuille profile at equilibrium --------------------
L      = ly-2;  
y_phys = y-1.5; 
 
ux = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys); 
uy = zeros(lx,ly); 
rho = 1; 
 
for i=1:9 
    cu = 3*(cx(i)*ux+cy(i)*uy); 
 
    fIn(i,:,:) = rho .* t(i) .*( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); 
end 
 
% MAIN LOOP (TIME CYCLES)--------------------------------------------------
for cycle = 1:maxT 
 
    % MACROSCOPIC VARIABLES
    rho = sum(fIn); 
    ux = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; 
    uy = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; 
 
    % MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS -------------------------
 
    % Inlet: Poiseuille profile -------------------------------------------
    y_phys = col-1.5; 
    ux(:,in,col) = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys); 
    uy(:,in,col) = 0; 
    rho(:,in,col) = 1 ./ (1-ux(:,in,col)) .* (sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,7,8],in,col))); 
 
    % Outlet: Constant pressure -------------------------------------------
    rho(:,out,col) = 1; 
    ux(:,out,col) = -1 + 1 ./ (rho(:,out,col)) .* (sum(fIn([1,3,5],out,col)) + 2*sum(fIn([2,6,9],out,col))); 
    uy(:,out,col) = 0; 
    % MICROSCOPIC BOUNDARY CONDITIONS: INLET (Zou/He BC)
    fIn(2,in,col) = fIn(4,in,col) + 2/3*rho(:,in,col).*ux(:,in,col); 
    fIn(6,in,col) = fIn(8,in,col) + 1/2*(fIn(5,in,col)-fIn(3,in,col)) ...
                                  + 1/2*rho(:,in,col).*uy(:,in,col) ...
                                  + 1/6*rho(:,in,col).*ux(:,in,col); 
                               
    fIn(9,in,col) = fIn(7,in,col) + 1/2*(fIn(3,in,col)-fIn(5,in,col)) ...
                                  - 1/2*rho(:,in,col).*uy(:,in,col) ...
                                  + 1/6*rho(:,in,col).*ux(:,in,col); 
      
    % MICROSCOPIC BOUNDARY CONDITIONS: OUTLET (Zou/He BC)
    fIn(4,out,col) = fIn(2,out,col) - 2/3*rho(:,out,col).*ux(:,out,col); 
 
    fIn(8,out,col) = fIn(6,out,col) + 1/2*(fIn(3,out,col)-fIn(5,out,col)) ...
                                    - 1/2*rho(:,out,col).*uy(:,out,col) ...
                                    - 1/6*rho(:,out,col).*ux(:,out,col); 
 
    fIn(7,out,col) = fIn(9,out,col) + 1/2*(fIn(5,out,col)-fIn(3,out,col)) ...
                                    + 1/2*rho(:,out,col).*uy(:,out,col) ...
                                    - 1/6*rho(:,out,col).*ux(:,out,col); 
    % COLLISION STEP ------------------------------------------------------
for i=1:9 
    cu = 3*(cx(i)*ux+cy(i)*uy); 
    fEq(i,:,:) = rho .* t(i) .* ( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); 
    fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:)); 
end 
 
% =========================================================================
% OBSTACLE (BOUNCE-BACK)
for i=1:9 
    fOut(i,bbRegion) = fIn(opp(i),bbRegion); 
end 
 
% =========================================================================
% STREAMING STEP
for i=1:9 
    fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]); 
end 
 
% =========================================================================
% VISUALIZATION
if (mod(cycle,tPlot)==1) 
    u = reshape(sqrt(ux.^2+uy.^2),lx,ly); 
    u(bbRegion) = nan; 
    imagesc(u');
    axis equal off; drawnow
end
end