%基于微粒间相互作用建模的仿真

%各变量的定义：mu--黏度μ,phi--浓度Φ，phi_w(t)--浓度Φw(t)，
%db--微粒直径d,dc--通道直径D,q--流量q，abw，abb--相互作用系数，
%do--有效直径Dopen(t),area_c--通道横截面积Ac，L--通道长度L,N--通道个数N
%C--系统阻尼C， rb(1)--覆盖在壁上的粒子数，p(1)--压降p,Vt--固相总量
%area_b--阻塞面积，delta_t--Δt， n_it--迭代次数， Nbb，Nbw--覆在微粒顶部和壁上的微粒数
%R(n)--流体阻抗Z(t)，gamma(n)--剪切率γ，
clear all;
%-设置图表格式
%set(0,'DefaultAxesBox','On');       %-设置坐标轴
%set(0,'DefaultAxesFontSize',16);        %-字体大小
%set(0,'DefaultAxesFontName','Arial');    %-字体
%set(0,'DefaultAxesLineWidth',2);         %-线宽 
%set(0,'DefaultLineLineWidth',2);
%set(0,'DefaultAxesGridLineStyle',':');      %-网格
%set(0,'DefaultTextFontSize', 16);          %-设置文本格式
%set(0,'DefaultTextFontName','Ariaf');
 %悬浮液参数
mu=1.006e-3;    %viscosity of water 黏度[kg m-1 s-1] 
db=1.0*1e-6;     %微粒直径 [m]
q=1.0*1e-9/60;     % 流量 [ul/min]
phi=1e-2;       %solid concentration颗粒的浓度 [-]
abw=1e10;      %壁与微粒之间相互作用系数 [m"-2]
abb=1e10;       %微粒间相互作用系数 [m^-2]
%通道参数
N=1;         % 通道个数
dc=5e-6;       %通道入口直径 [m]
area_c= pi*dc^2/4;    %通道横截面积[m^2]
L=5.5e-4;     %length of a channel [m]
C=5e-13;   %系统阻尼[m^4 s^2/kg]
%初始条件
rb(1)=0;    %覆盖在壁上的微粒初始值[-]
area_b(1)=rb(1)*pi*2/3*dc*db;  %resulting blocked area 阻塞面积[m^2]
p(1)=0;     %-初始压强值 [Pa]
Vt(1)=0;    %-初始固相体积
%iteration  循环 迭代
delta_t=1;  %时间步长值 time step[s]
n_it=1000;      %迭代次数 number of iterations

for n=1:n_it;
  if area_b(n)>area_c   % 如果通道已经阻塞 
    do(n)=do(n-1);      % 有效直径已经为0，Dopen(t)---->0
    R(n)= R(n-1);   %通道的流体阻抗为无穷大，Z(t)----->∞
  else
    do(n)=2*((area_c-area_b(n))/pi)^(1/2);      %有效直径Dopen[m]
    if do(n)>0
    R(n)= mu*128/pi*L/dc(n)^4;      % 通道的流体阻抗
    else
      R(n)=R(n-1);
    end
  end
  qc(n)=N*(q-p(n)/R(n));    %流向阻尼的流量[m^3/s]
      p(n+1)= p(n)+qc(n)*delta_t/C;      % 压降增量 [Pa]
  Vt(n+1)=Vt(n)+N*phi*p(n)/R(n)*delta_t;    % 通过通道的颗粒的总量
  gamma(n)=32/pi*(p(n)/R(n))/do(n)^3;     % 剪切率 measure for the shear rate[1/s]
  if do(n)<=db   %如果阻塞，  
  phi_w(n)=phi;    %近壁区的浓度-阻塞=bulk
  else
    phi_w(n)=do(n)/db*phi/(2*do(n)/db-1)^(1/2); % 近壁区浓度 
  end
  Nbw(n)=abw*db^2*gamma(n)*(1-rb(n))*pi*dc/db*phi_w(n)*delta_t; % 覆在壁上的粒子数
  Nbb(n)=abb*db^2*gamma(n)*rb(n)*pi*do(n)/db*phi_w(n)*delta_t;% 覆在原来粒子上顶的粒子数
rb(n+1)=rb(n)+db/(pi*dc)*Nbw(n); % 新的覆盖率
   if rb(n+1)>1    %  ？？？覆盖率>=1时，阻塞已经形成？
  rb(n+1)=1;
   end
  area_b(n+1)=area_b(n)+pi/4*db^2*(Nbw(n)+Nbb(n));%新的阻塞区域面积
end
%Vt(n+1);
%画出结果
time_n1=linspace(0, n_it*delta_t, n_it+1);
figure%(1):
clf;
hold on;
plot(time_n1,p/1000,'b-');
axis([0 1000 0 1 ]);
grid on;
xlabel('t [sec]');
ylabel('p [kPa]')