以梦为马_ck
2021-07-23 21:51
采纳率: 100%
浏览 151

MATLAB实现长方体容器内球形颗粒堆积

右手坐标系;容器的中心点为原点;
容器尺寸;球的体积分数;球的半径可以调整;
有过程的动态演示;
最后可以输出球的位置信息并保存。
PS:希望能写注释比较详细!

  • 点赞
  • 收藏

3条回答 默认 最新

  • soar3033 2021-07-24 11:07
    已采纳

    动态绘图,并存储。请采纳

    img

    img

    img

    xl=200;%长
    yl=200;%宽
    zl=100;%高
    xmin=-xl/2;%计算x轴最负值,因为0是原点。下同
    ymin=-yl/2;
    zmin=-zl/2;
    r=10;%球半径
    n=40;%球数量
    for i = 1:n%循环生成初始位置
        if i==1%如果是第一个,直接随机生成
            x(i)=rand()*(xl-2*r)+xmin+r;
            y(i)=rand()*(yl-2*r)+ymin+r;
            z(i)=rand()*(zl-2*r)+zmin+r;
        else%如果不是第一个
            while 1%进入死循环,并生成位置
                x(i)=rand()*(xl-2*r)+xmin+r;
                y(i)=rand()*(yl-2*r)+ymin+r;
                z(i)=rand()*(zl-2*r)+zmin+r;
                st=0;%标记
                for j=1:i-1%遍历所有已生成的球心
                    if ((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2)<(2*r)^2 %如果干涉
                        st=1;%更改标记
                        break    %跳出for循环
                    end    
                end
                if st==0%如果for循环完毕后,st依然为0,表明没有干涉
                    break%则推出死循环
                end
            end
        end
    end
    while 1%循环,用于绘制动态球下落
        st2=0;%标记,用于识别是否所有球都下落完
        for i= 1:n%遍历球
            st=0;%标记
            for j=1:n%遍历球
                if i~=j%如果ij不是一个球,则进行判断,如果两球干涉,且i球高于j球,则停止下落
                    if ((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2)<(2*r)^2 && z(i)>z(j)
                        st=1;%修改下落标记
                        break    
                    end
                end
                if z(i)-zmin<r %如果i球已经到底面
                    st=1;%修改下落标记
                    break
                end    
            end
            if st==0%如果标记位0,则下落
                z(i)=z(i)-1;%下落计算,步长为1,可更改,越小越慢
                st2=1;%更改st2标记,表明有球下落,循环继续
            end
        end
        if st2==0%如果无球下落,则跳出循环
            break
        end
        for i =1:n%遍历所有球绘图
            [xx,yy,zz]=ellipsoid(x(i),y(i),z(i),r,r,r);%设定绘制的球心和半径
            surf(xx,yy,zz)%绘制球
            if i==1%当绘制第一个球的时候
                hold on%打开hold on,否则只显示一个球
                axis([xmin,xmin+xl,ymin,ymin+yl,zmin,zmin+zl])%设置坐标轴
            end
        end
        hold off%每次绘制完一轮球,关闭hold,以清屏
        pause(0.01)%延时
    end
    aa=[x;y;z];%把球的坐标输入到aa
    t=table(aa);%把aa变成表t
    tmp={'x','y','z'};%设置表头
    t.Properties.RowNames=tmp;%设置表头
    writetable(t,'1.xlsx','WriteRowNames',true);%输出到excel保存
    
    点赞 1 打赏 评论
  • 王大师王文峰 2021-07-24 01:25

    采纳一下,谢谢

    function spheres
    clc,clear;
    n=50;  %球的个数
    x1=0;x2=100;  %x轴范围
    y1=0;y2=100;  %y轴范围
    z1=0;z2=20;   %z轴范围
    r=2;   %球的半径
    A=[x2-x1;y2-y1;z2-z1];
    A=A-2*r;
    centre=zeros(3,n-1);%保存圆心坐标
    for i=1:n
        c=rand(3,1);
        c=c.*A; %圆心坐标
        c=c+[x1;y1;z1]+r;
        if i==1
            [x,y,z]  = ellipsoid(c(1),c(2),c(3),r,r,r);
            surf(x,y,z,ones(size(x))) %画出来球
            hold on
        else
            s=0;
            while(s==0)
                c=rand(3,1);
                c=c.*A;    %圆心坐标
                c=c+[x1;y1;z1]+2;
                s=judge(i,c,centre,r);   %判断圆心是否可以  即新画出的球不会与之前的重叠
            end
            [x,y,z]  = ellipsoid(c(1),c(2),c(3),r,r,r);
            surf(x,y,z,ones(size(x))) %画出来球
            hold on
        end
         
    end
     
    axis([0 100 0 100 0 20]);
    function s=judge(i,c,centre,r)
    for j=1:i-1
        temp=centre(:,j);
        f=c-temp;
        d=norm(f,2);
        if d<2*r
            s=0;  %表示新圆心不可以在c表示的点
           return
        end
    end
    s=1;%表示新圆心可以
    
    
    
    点赞 打赏 评论

相关推荐 更多相似问题