guanghua.m

用matlab编程模拟光滑表面的气体流动效应

function RenxiaoTxtr2
clear,clc;close all;
CurrentWP = pwd;
%第一步生成一个高斯分布白噪声二维随机序列，均值为0，方差为1，单位微米，长宽都为252微米。计算其复立叶变换
rand('seed',2008);
eta = normrnd(0,0.1,252,252);
%  zxg=xcorr2(eta);
A = fft2(eta);
% 第二步根据指定的自相关函数，得到滤波器输出信号的功率谱密度，同时确定输入序列的功率谱密度...
% βx,βy表示沿x轴,y轴方向上的表面轮廓的自相关长度,单位微米。各向同性的粗糙表面,βx=βy=β
% 表面高度标准差0.05微米。模拟表面的长宽都是252微米。
for betay =30;
    R = zeros(252,252);
    betax =3;        
    for taox = 1:252
        for taoy = 1:252  
            R(taox,taoy) = (0.1)^2*exp(-2.3*sqrt(((taox-126)/betax)^2+((taoy-126)/betay)^2));
        end
    end
    Rtxty = R;
    Gz = fft2(Rtxty);
    C = 0.01; % C等于(0.5)^2 ，等于输入序列eta的方差
    %第三步计算滤波器的传递函数
    Hwxwy = sqrt(Gz/C);
    % 矩阵的各元素求平方根, sqrtm是矩阵求平方根
    %第四步计算输入序列经过二维滤波器后的输出序列的富丽叶变换
    Zwxwy = Hwxwy .* A;
    %第五步计算滤波器的传递函数对其进行富丽叶逆变换得到表面的高度分布函数
    zwxwy = ifft2(Zwxwy);
    Z=real(zwxwy);
%     xg=xcorr2(zwxwy);
    H=0.8*ones(252,252); % 定平均膜厚是1微米                        
    z=H+Z;     % z代表局部膜厚
%      for i=1:252;
%       for j=1:252;
%         if Z(i,j)>=-0.1   
%             Z1(i,j)=z(i,j);
%         else  Z1(i,j)=-0.09;  
%         end
%     h=Z1;
%     L0(i,j)=(h(i,j))^3;  
%       end
%      end  
%     l(i,j)=(h(i,j))^3;
    for i=1:252;
      for j=1:252;
        if z(i,j)>0
            z1(i,j)=z(i,j);
        else  z1(i,j)=0.01;  
        end
      end
    end  
end 
hh=z1;
% hh=H
% HH=reshape(hh,252*252,1)
% save HH.dat HH -ascii; 
%上述生成粗糙表面数据
lnx=252;                     
lny=251 ;                    
dx=1;     %划分的网格间距为1微米
dy=1;
%定义压强的的初值
pin=1;%入口压力，单位是帕斯卡
pout=0.5; %出口压力 
% maxl=100;
% I=250;
% J=250;
% hx=1/I;
% hy=1/J;

% for j=1:250
%  for i=1:250
%   X0(i,j)=pin-(i-1)*(pin-pout)/(250-1);
%   p(i,j)=X0(i,j);
%  end
% end

for j=1:lny+1
 for i=1:lnx
  ppo(i,j)=pin-(i-1)*(pin-pout)/(lnx-1);
  pp(i,j)=ppo(i,j);
 end
end
 for i=1:lnx
  hh(i,1)= hh(i,3);
  hh(i,lny+1)=hh(i,lny-1);
 end
 
 kiter=1;
err=1;  
if err>0.0000000000001
   for j=2:lny
      for i=2:lnx-1
%            b=((Pi0_j+P0i_j)/hx^2+(Pi_j0+Pi_0j)/hy^2+3/hi_j*((hi0_j-h0i_j)*(Pi0_j-P0i_j)/(2*hx)^2+...
%                 (hi_j0-hi_0j)*(Pi_j0-Pi_0j)/(2*hy)^2));
          
          b1=(ppo(i+1,j)+ppo(i-1,j))/dx^2;
          b2=(ppo(i,j+1)+ppo(i,j-1))/dy^2;
          a=-2./dx^2-2./dy^2;
          c=(ppo(i+1,j)-ppo(i-1,j))^2/(4*dx^2)+(ppo(i,j+1)-ppo(i,j-1))^2/4/dy/dy;  
          bx=((hh(i+1,j)-hh(i-1,j))*(ppo(i+1,j)-ppo(i-1,j)))/(4*dx^2); 
          by=((hh(i,j+1)-hh(i,j-1))*(ppo(i,j+1)-ppo(i,j-1)))/(4*dy^2); 
          b3=3*(bx+by)/hh(i,j);    
          b=b1+b2+b3;
          pp(i,j)=(-b-sqrt(b*b-4.*a*c))/2/a;  
%           if pp(i,j)>1
%             pp(i,j)=pp(i+1,j);
%         else  pp(i,j)= pp(i,j);  
%         end
      end
         pp(1,j)=pin;
         pp(lnx,j)=pout;
   end
         for i=1:lnx   
          pp(i,1)=pp(i,3);
          pp(i,lny+1)=pp(i,lny-1);
         end  
%        err=0.001;
         for j=1:lny+1
           for i=1:lnx
%            err=max(abs((ppo(i,j)-pp(i,j))/pp(i,j)));
           err=abs((ppo(i,j)-pp(i,j))/pp(i,j));
           ppo(i,j)=pp(i,j);
          %  d(i,j)=(ppo(i,j))^3;
        L(i,j)=(hh(i,j))^3;
            end
         end 
%      for i=1:252;
%       for j=2:252;
% %          ppo(i,j)>pin
%         if ppo(i,j)>pin 
%         ppp(i,j)=(pin-pout)/2; 
%         elseif ppo(i,j)<pout
%          ppp(i,j)=(pin-pout)/2;
%         else  ppp(i,j)=ppo(i,j); 
%         end        
%       end 
%      end       
   for j=1:lny
      for i=2:lnx-1
     M(i,j)=(ppo(i+1,j)-ppo(i-1,j))/2/dx;  

      end
   end    
   for j=1:lny
   M(lnx,j)=M(lnx-1,j);
   end
   for i=1:lnx     
   M(i,lny+1)=M(i,lny);
   end 
%    for j=1:lny+1
%    for i=1:lnx
%      if abs(M(i,j))>=0.002;  
%      M(i,j)=-0.002; 
%      else
%       M(i,j)= M(i,j);   
%      end   
%    end
%    end
   
 for j=1:lny+1
      for i=1:lnx
  N(i,j)=L(i,j)*M(i,j);
    N(i,j)=L(i,j)*ppo(i,j)*M(i,j);
      end
 end
%  S=sum(N(:)); 
S=sum(N); 
 Pm=mean(ppo(:));%平均压力。
 fai=S/(Pm*0.8^3*(pout-pin));
%   fai=S/0.2^3/(pout-pin);   %Φ是流量因子
  
% if(abs(kiter,100)=0)then  
% write(*,*)kiter,err 
%  do i=1,lnx
%  write(5,*)ppo(i,25)
% %  enddo
% % !call Plot(pp,lnx,lny+2,kiter)
%  endif   
  kiter=kiter+1   
end   
% M=(ppo(i+1,j)-ppo(i-1,j))^2/(4*dx^2)