通信人家园

 找回密码
 注册

只需一步,快速开始

短信验证,便捷登录

搜索

军衔等级:

  新兵

注册:2011-3-28
跳转到指定楼层
1#
发表于 2011-3-28 21:56:33 |显示全部楼层 |倒序浏览
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A LATTICE BOLTZMANN FORMULATION FOR THE ANALYSIS OF
% RADIATIVE HEAT TRANSFER PROBLEMS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in matlab
% Copyright (C) 2010-12-10 14:53:49 Zhao Lei .
% Address:
% E-mail: haorizi2008.cool@163.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clf
% RADIATION CONSTANS
ly     = 51;
lx     = 501;                        
delta_x= 1./(lx-1);
delta_y= 1./(ly-1);
delta_t= 1.0e-004;
sigma = 5.67e-008;          % constant s is the geometric distance
Thot  = 1;
Tcold = 0;
Ib    = sigma*Thot.^4/pi ;  % blackbody intensity
G     = 4*pi.*Ib ;          % volumetric absorption
IEq   = G/(4*pi);           % equilibrinm PDFs
beta  = 1.0;                % extinction coefficient
n   =8;
alpha = pi/4;      %
tPlot = 100;       % iterations between successive graphical outputs
tStatistics = 10;  % iterations between successive file accesses
Kn = delta_x*beta <= 0.05;  % Kn number is smaller than a threshold value that ensures stability and accuracy
% D2Q8 LATTICE CONSTANTS
cx  = [  1,   0,   -1,   0,   1,   -1,  -1,   1  ];
cy  = [  0,   1,    0,  -1,   1,   1,   -1,  -1  ];
t   = ones(1,8)*alpha/(2*pi);       % weight
az  = [0:2*pi/n:2*pi];      
[y,x] = meshgrid(1:ly,1:lx);

ux  = delta_x/delta_t;   % velocity
uy  = delta_y/delta_t;
omega = sqrt((cx*ux)*(cx*ux)'+(cx*uy)*(cx*uy)')*beta;    % relaxation frequency

% INITIAL CONDITION FOR RADIATION:  
IIn  = ones(8,lx,ly);   %8 is the direction number,2 is r
IOut = ones(8,lx,ly);
ts = 0;
% MAIN LOOP
while (ts<80000 && Kn <= 0.05 && abs((abs(IOut)-abs(IIn))./abs(IIn)))<1e-6 || ts<100
      
    %  BOUNDARY CONDITIONS FOR RADIATION:
    for i=1:8
        IIn (i , : ,1) = sigma*Thot^4/pi;  % sigma is Stefan-Boltzmann constant
        IOut(i , 1 ,:) = 0;
        IOut(i ,lx ,:) = 0;
        IOut(i ,: ,ly) = 0;
    end
      % COLLISION STEP
    for i=1:8
       % IEq   = zeros(1,lx,ly);  %让平衡分布函数IEq初始为同维的零矩阵
        IEq (i,:,:)  = IEq(i,:,:)+ reshape(t*reshape(IIn,8,lx*ly),1,lx,ly);       % equilibrinm PDFs
       % IOut(i,:,:) = zeros(1,lx,ly);
        IOut(i,:,:) = IIn(i,:,:) - omega .* (IIn(i,:,:)-IEq(i,:,:));
    end
   
      % STREAMING STEP RADIATION
    for i=1:8
       IIn(i,:,:) = circshift(IOut(i,:,:), [0,cx(i),cy(i)]);
       tx  = pi*sin(alpha/2)*cos(az(i));     %weight x direction
       ty  = pi*sin(alpha/2)*sin(az(i));     %weight y direction
    end
    ts = ts+1;
end
    % radiative heat fluxes
    Qx = sum(tx.*IOut);
    Qy = sum(ty.*IOut);
  
    % VISUALIZATION
    if (mod(ts,tPlot)==0)
       Q = reshape(sqrt(Qx.^2+Qy.^2),lx,ly);
       imagesc(Q');colorbar
       title('heat flux');
       axis equal off; drawnow
    end


黑体部分我想加上,老是出错,collision step 中 IEq (i,:,:)  = IEq(i,:,:)+ reshape(t*reshape(IIn,8,lx*ly),1,lx,ly);   出错,不会改啊,忘各位伸出援手

举报本楼

您需要登录后才可以回帖 登录 | 注册 |

版规|手机版|C114 ( 沪ICP备12002291号-1 )|联系我们 |网站地图  

GMT+8, 2025-7-21 13:15 , Processed in 0.148879 second(s), 18 queries , Gzip On.

Copyright © 1999-2025 C114 All Rights Reserved

Discuz Licensed

回顶部