% The MIT License (MIT) % Copyright (c) 2016 Andrea Chiarelli, Andrew Dawson, Alvaro Garcia % The University of Nottingham % Permission is hereby granted, free of charge, to any person obtaining a % copy of this software and associated documentation files (the "Software"), % to deal in the Software without restriction, including without limitation % the rights to use, copy, modify, merge, publish, distribute, sublicense, % and/or sell copies of the Software, and to permit persons to whom the % Software is furnished to do so, subject to the following conditions: % The above copyright notice and this permission notice shall be included % in all copies or substantial portions of the Software. % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS % OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, % FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE % AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER % LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING % FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER % DEALINGS IN THE SOFTWARE. clc clear all close all %%%%% INPUT VARIABLES START %%%%% % Parameter to choose based on your needs target_planar_void_ratio=40; % Planar void content as a % % Initialisation of the domain: CIRCLE % All sizes in metres (m) origin_of_cartesian_axes=0; domain_diameter=0.102; % Number of packed domains you wish to generate number_of_domains=1; for domain=1:1:number_of_domains % Number of particles in the first generation seeded number_of_particles=500; % Particles added in each new generation new_particles_per_generation=500; % Maximum size of a particle maximum_radius=0.0035; % Initial size of the particles in metres (m) minimum_radius_seeded=0.1e-4; %%%%% INPUT VARIABLES END %%%%% % Generation of a vector of radii, one row for each particle radius=minimum_radius_seeded*ones(number_of_particles,1); % Generation of the centres with random positions, one row for each particle centre(:,1)=minimum_radius_seeded+(domain_diameter-2*minimum_radius_seeded)*rand(number_of_particles,1); centre(:,2)=minimum_radius_seeded+(domain_diameter-2*minimum_radius_seeded)*rand(number_of_particles,1); % Generation of the angle for the equation for parametric circumferences angles_for_parametric_circle=0:0.01:2*pi; % Initialisation of the first generation of particles x_coordinates_circles=zeros(length(angles_for_parametric_circle), number_of_particles); y_coordinates_circles=x_coordinates_circles; % Counter for later use counter=number_of_particles; % Beginning of the packing process tic % For timekeeping iteration=2; planar_void_ratio(2)=100; % This needs to start at 2 because in the loops below it compares the current void ratio to its previous value while planar_void_ratio(iteration)~=planar_void_ratio(iteration-1) && planar_void_ratio(iteration)>target_planar_void_ratio iteration=iteration+1; for k=1:number_of_particles % Condition 1: Is the distance between particle k and all other % particles smaller than the sum of their radii plus the minimum % radius growth? This must return 1, because in order to grow this % condition must be 0 when particle k is compared to all other % particles but when it will be compared to itself the result will % be 1. condition(1)=sum(sqrt((centre(:,1)-centre(k,1)).^2+(centre(:,2)-centre(k,2)).^2)<(radius(k)+minimum_radius_seeded+radius)); % Conditions 2: If the radius grows by 'minimum_radius_seeded', % will it still be inside the circle defined above? condition(2)=sqrt((centre(k,1)-domain_diameter/2)^2+(centre(k,2)-domain_diameter/2)^2)target_planar_void_ratio && quit_the_calculation==0 i_new_generation=counter; if planar_void_ratio(iteration)<35 && increase_particles_per_generation==0 new_particles_per_generation=new_particles_per_generation+1000; increase_particles_per_generation=1; end if planar_void_ratio(iteration)<25 && increase_particles_per_generation==1 new_particles_per_generation=new_particles_per_generation+5000; increase_particles_per_generation=2; end if planar_void_ratio(iteration)<15 && increase_particles_per_generation==2 new_particles_per_generation=new_particles_per_generation+10000; increase_particles_per_generation=3; end % The next condition makes sure that no new generation of particles is empty while counter==i_new_generation for q=i_new_generation+1:i_new_generation+new_particles_per_generation % A new particle is initialised centre(q,1)=minimum_radius_seeded+(domain_diameter-2*minimum_radius_seeded)*rand; centre(q,2)=minimum_radius_seeded+(domain_diameter-2*minimum_radius_seeded)*rand; radius(q)=minimum_radius_seeded; % Is the particle acceptable? If yes, it is added to the vectors of % centres and radii if sum(sqrt((centre(:,1)-centre(q,1)).^2+(centre(:,2)-centre(q,2)).^2)<(radius(q)+minimum_radius_seeded+radius))==1 counter=counter+1; centre(counter,:)=centre(q,:); radius(counter)=radius(q); end end end % The vectors of centres and radii are cleaned from the unused values centre(counter+1:end,:)=[]; radius(counter+1:end)=[]; iteration=iteration+1; % This is re-initialised due to the counters used below planar_void_ratio(iteration)=100; while planar_void_ratio(iteration)~=planar_void_ratio(iteration-1) && quit_the_calculation==0 iteration=iteration+1; for k=i_new_generation+1:counter % For the meaning of these conditions see above. condition(1)=sum(sqrt((centre(:,1)-centre(k,1)).^2+(centre(:,2)-centre(k,2)).^2)<(radius(k)+minimum_radius_seeded+radius)); condition(2)=sqrt((centre(k,1)-domain_diameter/2)^2+(centre(k,2)-domain_diameter/2)^2)