matlab生成50个三维随机骨料模型代码，骨料为八面体，并以重力形式自然堆叠在一立方米的模型内，并展示最终成果

以下是一个可能的MATLAB代码实现，用于生成50个三维随机骨料模型，并以重力形式自然堆叠在一个立方体模型内：

% 设置模型参数
num_models = 50;   % 模型数量
num_particles = 500;   % 每个模型中的粒子数量
particle_size = 0.1;    % 粒子直径
gravity = [0, 0, -9.8]; % 重力加速度
box_size = 1;    % 立方体模型尺寸

% 预分配存储空间
models = cell(num_models, 1);

% 循环生成模型
for i = 1:num_models
    % 生成随机粒子位置
    positions = (rand(num_particles, 3) - 0.5) * box_size;
    
    % 生成八面体形状的粒子
    [x, y, z] = sphere(6);
    x = x(:) * particle_size / 2;
    y = y(:) * particle_size / 2;
    z = z(:) * particle_size / 2;
    vertices = [x, y, z];
    faces = delaunay(x, y, z);
    
    % 随机旋转八面体粒子
    rotations = rand(num_particles, 3) * 360;
    quaternions = eul2quat(rotations, 'XYZ');
    rotated_vertices = quatrotate(quaternions, vertices);
    
    % 将旋转后的粒子加入模型中
    model = struct('vertices', [], 'faces', []);
    for j = 1:num_particles
        model.vertices = [model.vertices; positions(j, :) + rotated_vertices];
        model.faces = [model.faces; faces + (j - 1) * size(vertices, 1)];
    end
    
    % 自然堆叠粒子
    model.vertices = natural_stacking(model.vertices, model.faces, gravity);
    
    % 存储模型
    models{i} = model;
end

% 自然堆叠函数
function stacked_vertices = natural_stacking(vertices, faces, gravity)
    % 计算每个面的法向量
    normals = cross(vertices(faces(:, 2), :) - vertices(faces(:, 1), :), ...
                    vertices(faces(:, 3), :) - vertices(faces(:, 1), :));
    face_centers = (vertices(faces(:, 1), :) + vertices(faces(:, 2), :) + vertices(faces(:, 3), :)) / 3;
    
    % 初始化每个粒子的速度和加速度
    velocities = zeros(size(vertices));
    accelerations = repmat(gravity, size(vertices, 1), 1);
    
    % 模拟时间步长
    dt = 0.01;
    for i = 1:1000
        % 计算每个粒子的受力
        forces = zeros(size(vertices));
        for j = 1:size(faces, 1)
            force = normals(j, :) * dot(gravity, normals(j, :));
            forces(faces(j, :), :) = forces(faces(j, :), :) + repmat(force / 3, 3, 1);
        end
        
        % 更新速度和位置
        velocities = velocities + accelerations * dt;
        vertices = vertices + velocities * dt;
        
        % 计算每个粒子的碰撞检测
        for j = 1:size(vertices, 1)
            if vertices(j, 3) < 0
                vertices(j, 3) = 0;
                velocities(j, :) = -velocities(j, :) * 0.5;
            end
        end
        
        % 更新加速度
        accelerations = forces / repmat(sum(forces, 2), 1, 3);
    end
    
    stacked_vertices = vertices;
end

% 显示模型
figure;
axis equal;
hold on;
for i = 1:num_models
    model = models{i};
    patch('Vertices', model.vertices, 'Faces', model.faces, 'FaceColor', 'r', 'EdgeColor', 'none');
end
xlim([-box_size/2, box_size/2]);
ylim([-box_size/2, box_size/2]);
zlim([0, box_size]);
xlabel('X');
ylabel('Y');
zlabel('Z');
view(3);

这个代码将生成50个随机的三维骨料模型，其中每个模型包含500个八面体形状的粒子。粒子的位置和旋转是随机生成的，并且每个模型都以重力自然堆叠在一个立方体模型内。最后，代码将显示所有模型的最终成果。

运行代码后，您应该会看到一个立方体模型，其中包含了50个自然堆叠的骨料模型。每个模型都是以红色呈现的，立方体模型的边界是蓝色的。