Matlab生成不同粒径和形状的二维离散颗粒随机分布

在岩土力学、泥沙运行、多相流颗粒物输送问题，常采用基于无网格法的光滑粒子法（SPH）或离散元法，而这些方法常需要对所研究问题建立初始的颗粒随机分布，在此以一个简单的案例加以说明，仅作抛砖引玉：

案例基本情况：10*10区域，随机均布50*50=2500个粒子，坐标按随机产生，颗粒为4种类型，大小分为5种，然后用matlab编制小程序进行图形绘制。

主要程序如下：
clc;clear;close all;

Nx=50;
Ny=50;
N=Nx*Ny;  %点的数量
Lx=10;  %x方向长度
Ly=10;  %y方向长度
NGeo=4; %颗粒的形状种类数量
Nrad=5;   %颗粒的粒径种类数量
Nrank=3;  %颗粒的级配种类数量
rbanjing=50;  %颗粒基准半径大小

a=rand(N,1); %在0-1之间获取随机数，共N个数  
AAx0 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标x
AAy0 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标y

AAx1 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标x
AAy1 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标y

AAx2 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标x
AAy2 = rand(Nx,Ny);  %生成Nx*Ny的均布散点坐标y

AAx =AAx0;  %生成Nx*Ny的均布散点坐标x
AAy =AAy0;  %生成Nx*Ny的均布散点坐标y

aLx=reshape(AAx,[N,1])*Lx; %获取x方向的粒子随机坐标
aLy=reshape(AAy,[N,1])*Ly; %获取y方向的粒子随机坐标

b1 = randi([1,NGeo],N,1); %获取颗粒形状种类分布
b2 = rbanjing*randi([1,Nrad],N,1); %获取颗粒粒径种类分布
b3 = randi([1,Nrank],N,1); %获取颗粒级配种类分布

b4=num2str(zeros(N,1));  %生成初始字符数组
Bb=num2str(zeros(NGeo,1)); 
kk1=0;
kk2=0;
kk3=0;
kk4=0;

for i=1:N
    
  if b1(i) == 1
      b4(i)='o';       %圆形
      kk1=kk1+1;
      aLx1(kk1)=aLx(i);
      aLy1(kk1)=aLy(i);
      b21(kk1)=b2(i);
      b41(kk1)=b4(i);
      Bb(1)='o';
  elseif b1(i) == 2
      b4(i)='d';       %菱形
      kk2=kk2+1;
      aLx2(kk2)=aLx(i);
      aLy2(kk2)=aLy(i);
      b22(kk2)=b2(i);
      b42(kk2)=b4(i);  
      Bb(2)='d';
  elseif b1(i) == 3
      b4(i)='s';       %方形
      kk3=kk3+1;
      aLx3(kk3)=aLx(i);
      aLy3(kk3)=aLy(i);
      b23(kk3)=b2(i);
      b43(kk3)=b4(i);      
      Bb(3)='s';
  elseif b1(i) == 4    
      b4(i)='h';       % 六角形
      kk4=kk4+1;
      aLx4(kk4)=aLx(i);
      aLy4(kk4)=aLy(i);
      b24(kk4)=b2(i);
      b44(kk4)=b4(i); 
      Bb(4)='h';
  end    
end
figure;
set(gcf,'outerposition',get(0,'screensize'));

for ii=1:NGeo  %根据粒子种类分别绘制散点图

    if ii==1
      scatter(aLx1',aLy1',b21',Bb(1)); hold on;
    elseif ii==2
      scatter(aLx2',aLy2',b22',Bb(2)); hold on;
    elseif ii==3
      scatter(aLx3',aLy3',b23',Bb(3)); hold on;
    elseif ii==4
      scatter(aLx4',aLy4',b24',Bb(4));
    end

end