generate random size sphere which are randomly located inside a cylinder.

i want to generate the random size shpere which are randomly located inside a cylinder. distribution follow the gaussian distribution for both location and radius of sphere. Any one please help me on this problem?

Answer 1

[Walter Roberson]

Code is at https://www.mathworks.com/matlabcentral/answers/461635-randomly-generated-spheres-in-a-volume-in-matlab#answer_374712

You might need to adjust the details of the generation to be gaussian.

[M bd]

thank you for you answer.

but above code work on cube, i need to change it inside a cylinder. how to change it?

[Walter Roberson]

Expand the comments. I showed how to do it in https://www.mathworks.com/matlabcentral/answers/461635-randomly-generated-spheres-in-a-volume-in-matlab#comment_1610158

[M bd]

can you please put exect code. i tried but showing some error.

[Walter Roberson]

wanted = 250;
CylRadius = 40; CylHeight = 100;
MaxR = 10;
XYZR = GetRandomSpheresCylinder(wanted, [CylRadius, CylHeight], MaxR);
figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:wanted
    thisR = XYZR(k,4);
    surf(X .* thisR + XYZR(k, 1), Y .* thisR + XYZR(k, 2), Z .* thisR + XYZR(k, 3));
end
xlim auto; ylim auto; zlim auto; axis equal

function XYZR = GetRandomSpheresCylinder(nWant, Size, MaxRadius)
% INPUT:
%   nWant:  Number of spheres
%   Size:  Dimension as [1 x 2] double vector, [r h]
%   MaxRadius: maximum Radius of spheres, peak at 1/2 this
% OUTPUT:
%   P:      [nWant x 3] matrix, centers
a = 4; b = 4;  %beta parameters

CR = Size(1);   %cylinder radius
CH = Size(2);   %cylinder height
sph_r2 = sqrt(rand(nWant,1)) * CR;
sph_angle = rand(nWant,1) .* 2*pi;
sph_size = betarnd(a,b, nWant,1) .* MaxRadius;
vertical = rand(nWant,1) .* (CH - 2*sph_size) + sph_size;
[x, y, z] = pol2cart(sph_angle, sph_r2, vertical);
XYZR = [x, y, z, sph_size];
end

[M bd]

Thank you so much @Walter Roberson

i wnat to ask another question.

let if we want to add volume fraction of sphere with respect to cylinder, then how can i proceed?

[Walter Roberson]

First you would have to decide whether you want overlaps to be permitted or not. The current code generates overlaps (a lot of them!) and would have to be redesigned a bit to remove overlaps.

If you do permit overlaps, then calculating the volume fraction accurately can be... difficult. The task becomes much easier if you are willing to quantize the volume.

Calculating volumes when overlaps is not permitted is easier.

Asking to generate without overlaps until a certain volume fraction is reached can take a long time.

If you do want to generate without overlaps, then you typically occupy the positions for the largest spheres quickly: you stop being able to fit large spheres fairly soon. You can continue generating on your size distribution, but the larger sizes will get rejected as not fitting. The longer you go, the more likely that you can only fit spheres on the smaller scale. And that means that your distribution goals are unlikely to be met... unless your distribution goals were "power law".

Generation of random objects (spheres) filling other objects (cylindar) tends to follow one of several algorithms:

• Each object is locked in place as soon as it is generated. Tentative new positions and sizes are generated for new objects, and if they do not fit there, then that tentative assignment fails and you go on to generate new tentative positions and sizes;
• Each object is locked in place as soon as it is generated. Tentative new positions and sizes are generated for new objects, and if they do not fit there, the program tries to find a different location where an object that size would fit. This tends to get closer to size distribution goals
• Objects are not locked in place. Tentative new positions and sizes are generated throughout the volume for new objects, and if they do not fit there, then the program tries moving existing objects around to make the new one fit, such as by exerting a "force" that tries to move existing objects (and that can push against other objects moving them... but eventually you hit the walls and that pushes back...
• Objects are not locked in place. Tenative new sizes are generated, but the object is "dropped" from the top and allowed to fall, pushing on existing objects, they move, etc., walls get hit, and so on. The advantage of this algorithm is that if there is not enough room to squeeze a new sphere into a lower level, that it just sits on top
• Each object is locked in place as soon as it is generated. Tenative new sizes are generated, but the object is "dropped" from the top and allowed to fall. When it touches something on fewer than 3 contact points beheath it, it can slide or tumble to fill an existing hole, but it never pushes anything out of the way. This algorithm can be easier than continually recalculating how hundreds of existing objects push against each other as room is to be made for a new object
• Objects are not locked in place. A new object is placed on top of an existing object, and then the containing object is "shaken" to see if there is room for it (or other earlier objects) to fall

Volume fraction and object size distribution goals are harder to reach if objects cannot move.

If objects can move, then smaller objects tend to migrate downwards and larger objects tend to migrate upwards: https://en.wikipedia.org/wiki/Granular_convection

[M bd]

i wnat to generate it without overlap.

[M bd]

@Walter Roberson can you please solve this issue to add volume fraction on above code without overlapping of sphere and also sphere should not touch the cylinder boundary?

[Walter Roberson]

and also sphere should not touch the cylinder boundary

Let the variable Margin be the minimum distance between any sphere and the cylinder boundary. Then

sph_r2 = sqrt(rand(nWant,1)) * (CR - Margin);
sph_angle = rand(nWant,1) .* 2*pi;
sph_size = betarnd(a,b, nWant,1) .* MaxRadius;
vertical = rand(nWant,1) .* (CH - 2*sph_size - 2*Margin) + sph_size + Margin;

[Walter Roberson]

can you please solve this issue to add volume fraction on above code without overlapping of sphere

Using which of the 6 algorithms I described is required?

If it is the either of the first two (objects locked in place) then you need to describe how you want to deal with the problem I discussed that 'that means that your distribution goals are unlikely to be met... unless your distribution goals were "power law".'

If it is any of the last 4 algorithms, involving moving objects around after they are generated, then my response would be "No, I am not willing to do that, it is too much work."

Question: is there a minimum gap between spheres? You said "without overlapping" but what about if they are touching?

[M bd]

@Walter Roberson it will be fine if any of the first two algorithm and better if they would not touching.

[M bd]

@Walter Roberson can you please solve this issue to add volume fraction on above code without overlapping of sphere with any of first two alogirithm?

[Walter Roberson]

@M bd

I had some rush dental treatment yesterday, and I spent the evening recovering. Now I am behind in responding to people, and it may take me a while to get around to writing the code. I would recommend that you get a start by taking the code I posted and start modifying it.

Given any existing set of sphere centers and their radiuses, if you have a tentative new point and radius for it, you can tell if the new point would overlap with the existing points if pdist2() of the distances between the centers is less than the sum of the proposed radius and the radius of the existing point.

overlaps = pdist2([newx, newy, newz], [existing_x(:), existing_y(:), existing_z(:)]) - new_radius - existing_radius(:)
if all(overlaps > 0); this proposal fits; end

[M bd]

@Walter Roberson i hope that you will recover soon.

On for your above suggestion, how to define existing coordinate and new coordinate in the code?

[Image Analyst]

@M bd you keep track as you place the spheres, perhaps something like

numberWanted = 30
numberPlaced = 0;
existing_x = zeros(1, numberWanted);
existing_y = zeros(1, numberWanted);
existing_z = zeros(1, numberWanted);
for k = 1 : 100000 % Try at least 100 thousand times before quitting.
    % Get the proposed location and radius.
    thisR = XYZR(k,4);
    newx = XYZR(k, 1);
    newy = XYZR(k, 2);
    newz = XYZR(k, 3);
    % See if it overlaps any prior ones
    if numberPlaced > 1
        overlaps = pdist2([newx, newy, newz], [existing_x(:), existing_y(:), existing_z(:)]) - new_radius - existing_radius(:)
    else
        overlaps = 1; % Definitely want to place the very first one.
    end
    if all(overlaps > 0)
        % This proposed location works because there are no overlapping spheres.
        numberPlaced = numberPlaced + 1; % Increment the number that we've placed.
        % Draw it.
        surf(X .* thisR + XYZR(k, 1), Y .* thisR + XYZR(k, 2), Z .* thisR + XYZR(k, 3));
        % Log this location as an existing one.
        existing_x(numberPlaced) = XYZR(k, 1);
        existing_y(numberPlaced) = XYZR(k, 2);
        existing_z(numberPlaced) = XYZR(k, 3);
    end    
end
% Crop to as many as we could place, in case we couldn't place as many as we wanted.
existing_x = existing_x(1 : numberPlaced);
existing_y = existing_y(1 : numberPlaced);
existing_z = existing_z(1 : numberPlaced);

[M bd]

@Image Analyst and @Walter Roberson how could i define existing radius and new radius as i didn't define it before.

[Walter Roberson]

Image Analyst's suggested code here acts to filter the XYZR coordinates produced by GetRandomSpheresCylinder() -- assuming that the nWant for the call to that function was at least 100000.

[Image Analyst]

The number wanted is independent of the number of tries to place it, however you need to try at least as many times as particles you have. But if you have say a rectangular box, and you managed to fit one sphere in it, it will keep trying to fit more in. However you don't want it to keep trying forever is there is just no way that any other sphere could possibly fit - that would just be an infinite loop. So it says to try 100000 times before giving up and assuming that no more could fit anywhere, and you'd be left with just the one in there. So in the example I said to place 30 spheres. But if you could only place 25, and after 100000 other trial locations, it was not able to place a single additional sphere, it would exit and you'd have just the 25 that did fit. Of course we'll have to try at least 30 times because if all spheres are going to be able to fit and there is never any overlap, we need to try that many times to make sure they're all placed. But no sense beating a dead horse -- if after 100 thousand tries you cannot find a location to fit in another one, it will quit.

Does that explain it?

[M bd]

@Image Analyst but how can i get volume fraction of sphere.

i want input as volume fration of sphere.

can you please put here whole code? i am getting error after placing your code on walters code.

[Walter Roberson]

Where Image Analyst coded

existing_z = zeros(1, numberWanted);

after that add

existing_R = zeros(1, numberWanted);

Then where he had

existing_z(numberPlaced) = XYZR(k, 3);

add after that

existing_R(numberPlaced) = XYZR(k, 4);

and where he coded

existing_z = existing_z(1 : numberPlaced);

add after that

existing_R = existing_R(1 : numberPlaced);

Now to calculate volume fraction, take

total_occupied_volume = sum(4/3*pi*existing_R.^3);

and referring back to the variables that were already defined that looked like

CylRadius = 40; CylHeight = 100;

then

total_cyl_volume = pi .* CylRadius.^2 .* CylHeight;

and now you can calculate the volume fraction.

If you need to calculate the volume fraction as you go, then you can keep a running total of the cube of existing_R() that has been assigned so far, after which the occupied volume is 4/3*pi times that sum of cubes.

[Walter Roberson]

can you please put here whole code?

He is not likely to do that. It is your assignment, and he already gave you links to places where people have posted sphere packing code, and he already gave you very solid outline of ensuring there are no overlaps and keeping the entries that are ok. He has done more than his share of your assignment.

If you post your code, with error, then one of the volunteers might explain to you why you got the error, but we are waiting for you to put in visible effort, so we are probably not going to fix the code for you, just tell you the cause of the error.

Answer 2

[Image Analyst]

Click on the tag links on the right for sphere packing and circle packing.