Measure Logic
Problem Statement
Note
The project file for this example may be viewed/run in PFC.[1] The data files used are shown at the end of this example.
The measure logic (see the measure create
command) is exercised with PFC to compute porosity. Two regular packings composed of monosized balls organized in crystallike patterns and a polydisperse packing are generated. The porosity is computed using measurement circles (or spheres in 3D) with different radii. The accuracy of the numerical results is evaluated through comparison to analytical values.
Closed Form solution
Two crystallike packings are considered. Both packings consist of equalradius balls in contact with no overlap. The first packing is a simple cubic arrangement (see Fig. 1 for a 2D view) in which each ball is in contact with four neighbors in 2D (six in 3D). The second packing is a hexagonal arrangement (see Fig. 2) in which each ball is in contact with six neighbors in 2D (twelve in 3D). The hexagonal arrangement is the closest of all possible regular twodimensional or threedimensional packings. For both packings, the porosity is independent of the size of the balls, provided that the packing is extended over a very large region so that the effect of boundaries is negligible (Refer to [Deresiewicz1958a] for a thorough discussion of packings in both two and three dimensions).
The density, \(D\), of a packing is defined as the ratio of the area of space occupied by solid matter, \(A_s\), to the total area, \(A\). If we assume that the total area comprises both solid and void areas and that there is no overlap, then the density is related to the porosity, \(n\), by the relation
For the simple cubic packing, we have, in two dimensions:
and, in three dimensions:
For the hexagonal packing, in two dimensions, we have:
and, in three dimensions:
In the case of a polydisperse assembly, a specific procedure is implemented to compute the porosity (see the get_poros function in “PolydisperseAssembly.dat (2D)” and “PolydisperseAssembly.dat (3D)”) using the ratio between the sum of the surface (volume in 3D) of each ball in the model over the total surface (volume in 3D) of the sample.
Numerical Model
Several parameterized FISH functions are used to generate the measurement circles (spheres in 3D) that will be used to compute the porosity of the packings.
These functions are defined in {”MeasureLogic.dat (2D)” for 2D; “MeasureLogic.dat (3D)” for 3D}.
The function define_meas_parameters takes two arguments: nmeas, whereby the user can define the number of measurement regions to be created within the packing; and pcntg, through which the tolerance value can be defined in the measure create
command. This last parameter determines the error (< pcntg) tolerance accepted during porosity computation in the measurement regions. See the reference page for details.
The measurement regions have random size and random position in the model. This is defined in the create_meas_regions function.
To verify the accuracy of the porosity computation, a check_error function with three arguments is defined: val, the reference value of the porosity, which is computed through Equations (1), (2), (3), (4) in the case of the regular assemblies, or directly using Equation (5) in the case of the polydisperse assembly, to be compared to the numeric results; t, a string variable which allows definition of the name of the table that will contain the values of the errors committed in the porosity measurement; and average_radius, the average radius of the balls comprising the assembly which will be used to normalize the diameters of the measurement regions while postprocessing the results of the simulation.
In the case of the regular assemblies, a parameterized FISH function is used to construct each packing using the ball generate
command in {”RegularAssembly.dat (2D)” for 2D; “RegularAssembly.dat (3D)” for 3D}. This function has three arguments: size, resolution, and packing type (size, nres, itype). Two models are built for each packing. These models contain 9 rows and 9 columns of balls (9x9x9 in 3D) and 30 rows and columns (30x30x30 in 3D), respectively.
A screenshot of the cubic and hexagonal arrangements (9x9 packings, low resolution, 2D view) is shown in Fig. 1 and Fig. 2.
To build a polydisperse assembly, a procedure is implemented in a FISH function whose sole argument is the size of the packing. The ball generate gauss
option is used, thus ensuring a particle size distribution to follow a Gaussian distribution. A cloud of nonoverlapping balls is generated. The functions are all defined in {”PolydisperseAssembly.dat (2D)” for 2D; “PolydisperseAssembly.dat (3D)” for 3D}.
Through the define_meas_parameters
function in “MeasureLogic.dat (2D)”, the measurement parameters are defined. Shown in Fig. 4 and Fig. 5 are screenshots of the measurement regions (nmeas = 5000) as they have been installed to measure the porosity of the cubic regular packing and the polydisperse packing.
Results and Discussion
The numerical results are compared to the reference values. The average errors committed using nmeas = 5000 measurement circles and pcntg = 0.0 are reported (Table 1 and Table 2).
In the case of regular assemblies, for the 2D model, we have:
Cubic Arrangement 
Hexagonal Arrangement 


Low resolution 
High Resolution 
Low resolution 
High Resolution 

9x9 
30x30 
9x9 
30x30 

Exact Porosity 
0.2146 
0.2146 
0.0931 
0.0931 
Average Error [%] 
3.87% 
0.88% 
3.64% 
0.74% 
Similar results are obtained in the 3D model:
Cubic Arrangement 
Hexagonal Arrangement 


Low resolution 
High Resolution 
Low resolution 
High Resolution 

9x9 
30x30 
9x9 
30x30 

Exact Porosity 
0.4764 
0.4764 
0.2595 
0.2595 
Average Error [%] 
2.24% 
0.26% 
1.71% 
0.20% 
The higher the resolution, the lower the average error committed in the porosity measurement. Raw data are stored in four tables (see the table
command reference page), one for each examined packing.
Fig. 6 shows the average error, plotted in a logarithmic scale as function of the diameter of the circle employed to measure the porosity, normalized by the average radius of balls. Again, it can be seen that the lowest errors (< 1%, negative values on the logarithmic axis) correspond to the highest resolutions and to the largest measurement regions.
We can make the following observations:
the measure of the porosity is locally always correct, even more so considering that a tolerance equal to 0 is employed;
the accuracy of the local measure compared to the global porosity, computed analytically, depends on the representativeness of the area intercepted by the measurement region with respect to the global area, in terms of proportion void/solid;
it is then clear how a larger measurement region has more chances to provide an accurate measure; and
it is even clearer how a higher resolution ensures the measurement regions to be more representative of the overall structure.
We do the same comparison in the case of the polydisperse assembly. The error is plotted in Fig. 7.
The relatively higher errors, with respect to the regular samples, are due to a nonhomogeneous value of porosity throughout the sample. Then, in this case, the employment of larger measurement regions ensures the measures to be less dispersed, even if it seems harder to get rid of a small error (\(\simeq 2\) %).
References
Deresiewicz, H. “Mechanics of Granular Matter”, in Advances in Applied Mechanics, Vol. 5, pp. 233306. H. L. Dryden and Th. Von Karman, eds. New York: Academic Press, 1958.
Data Files
MeasureLogic.dat (3D)
; fname: MeasureLogic.dat (3D)
;
; Exercise the Mechanical Ball integration to the Measure Logic with PFC.
; Porosity is computed using measurement spheres with random radii.
; Results for different model resolutions are compared against closed
; form solutions.
;
; Itasca Consulting Group, Inc.
;
; ========================================================================
fish def define_meas_parameters(nmeas, pcntg)
; ARGUMENTS:
;  nmeas : integer value that indicates the number of randomlyplaced
; probing measurement circles that will be created
; Vary this parameter to see the effects on the final result
;  pcntg : 'percentage' argument of the MEASURE function
;
global p = pcntg
global n = nmeas
end
fish def compute_meas_center(r)
;Randomly computes the coordinates of the measurement sphere center
;ARGUMENTS:
;  r : radius of the measurement sphere
;
global xmin, xmax, ymin, ymax, zmin, zmax
local xcmin = xmin+r
local xcmax = xmaxr
local ycmin = ymin+r
local ycmax = ymaxr
local zcmin = zmin+r
local zcmax = zmaxr
global xc = xcmin + (xcmax  xcmin)*math.random.uniform
global yc = ycmin + (ycmax  ycmin)*math.random.uniform
global zc = zcmin + (zcmax  zcmin)*math.random.uniform
end
fish define create_meas_regions
measradius = array.create(n)
local maxmeasradius = math.min((0.99*(zmaxzmin)), ...
math.min((0.99*(xmaxxmin)), ...
(0.99*(ymaxymin))))
loop ii(1,n)
measradius(ii) = math.max(maxmeasradius/15,...
math.random.uniform*(0.5*(maxmeasradius)))
end_loop
loop ii(1,n)
command
[compute_meas_center(measradius(ii))]
measure create posx [xc] posy [yc] posz [zc] ...
radius [measradius(ii)] tolerance [p]
end_command
end_loop
end
fish define check_error(val,t,averageradius)
; Compares porosities computed in the measurement circles
; with analytical values
; ARGUMENTS:
;  val: exact value of the porosity
;  t: name of the table with the final numerical results
;
global result = table.create(t)
local it = 0
local averageerror = 0
loop foreach local mp measure.list
it = it + 1
local res = measure.porosity(mp)
local diff = 100.0 * (res  val) / val
averageerror = averageerror + math.abs(diff)
table(t,math.abs(diff))=measradius(it)/averageradius
endloop
global check_error = averageerror/it
end
;=========================================================================
; EoF: MeasureLogic.dat (3D)
PolydisperseAssembly.dat (3D)
; fname: PolydisperseAssembly.dat (3D)
; Build a polydisperse granular assembly and compute porosity using
; the measure logic in PFC
; ========================================================================
;
model new
program call 'MeasureLogic'
fish def poly_assembly(size)
; Create a polydisperse packing in a square domain
; ARGUMENTS:
;  size : size of the domain
;
global xmin = 0.0
global ymin = 0.0
global xmax = size
global ymax = size
global zmin = 0.0
global zmax = size
global width = xmaxxmin ; width of box
global height = ymaxymin ; height of box
global depth = zmaxzmin ; depth of the box
global tot_vol = (width)*(height)*(depth)
local num = 10000 ; number of particles
local rlo = size/100
local rhi = 1.25*rlo
local oversize = rhi
command
model domain extent ([0.0oversize],[width+oversize]) ...
([0.0oversize],[height+oversize]) ...
([0.0oversize],[depth+oversize])
wall generate box (0.0, [width]) (0.0, [height]) (0.0, [depth])
end_command
; generate the balls and give them their properties
command
ball generate radius [rlo] [rhi] gauss ...
box ([width*0.02], [width*0.98]) ...
([height*0.02], [height*0.98]) ...
([depth*0.02], [depth*0.98]) ...
id=1,[num]
end_command
end
fish def averageradius
local rad = 0
loop foreach local bp ball.list
rad = rad + ball.radius(bp)
endloop
global averageradius = rad/ball.num
end
fish def get_poros ;computes the exact value of porosity
local sum = 0.0
loop foreach local bp ball.list
sum = sum + 1.33333* math.pi * ball.radius(bp)^3
end_loop
global get_poros = 1.0  sum / tot_vol
end
[poly_assembly(5.0)]
[get_poros]
[define_meas_parameters(5000, 0.0)]
[create_meas_regions]
fish def display_error(porosity, s)
global error = check_error(porosity,s,averageradius)
local oo = io.out(string.build("Average Error = %1%)",error))
end
[display_error(get_poros,'Polydisperse')]
model save 'PolydisperseAssembly'
; ========================================================================
; EoF: PolydisperseAssembly.dat (3D)
RegularAssembly.dat (3D)
; fname: RegularAssembly.dat (3D)
; Build a regular granular assembly and compute porosity
; using the measure logic in PFC
; ========================================================================
;
model new
program call 'MeasureLogic'
fish define regular_assembly(size, nres,itype)
; Create a regular packing of monosized balls in a cubic domain
; ARGUMENTS:
;  size : size of the domain
;  nres : resolution of the packing
;  itype: type of packing (0:simple cubic; #0: hexagonal)
;
global xmin = 0.0
global ymin = 0.0
global xmax = size
global ymax = size
global zmin = 0.0
global zmax = size
global radius = size / (2.0*nres)
local numballs = nres^2
command
model domain extent [xmin] [xmax] [ymin] [ymax] [zmin] [zmax]
end_command
if itype = 0 then
command
ball generate rad [radius] number [numballs] cubic
end_command
else
command
ball generate rad [radius] number [numballs] hexagonal
end_command
endif
end
fish def display_error(porosity, s)
global error = check_error(porosity,s,radius)
local oo = io.out(string.build("Average Error = %1%)",error))
end
[define_meas_parameters(5000,0.0)] ;define measurement parameters
; Simple Cubic packing  low resolution
[regular_assembly(18.0, 9, 0)]
[create_meas_regions]
[display_error(0.4764,'cubic_lores')]
model save 'cubic_lores'
;
; Simple Cubic packing  high resolution
ball delete
[regular_assembly(18.0, 30, 0)]
[display_error(0.4764,'cubic_hires')]
model save 'cubic_hires'
;
; Hexagonal packing  low resolution
ball delete
[regular_assembly(18.0, 9, 1)]
[display_error(0.2595,'hexa_lores')]
model save 'hexa_lores'
;
; Hexagonal packing  high resolution
ball delete
[regular_assembly(18.0, 30, 1)]
[display_error(0.2595,'hexa_hires')]
model save 'hexa_hires'
model save 'RegularAssembly'
; ========================================================================
; EoF: RegularAssembly.dat (3D)
MeasureLogic.dat (2D)
; fname: MeasureLogic.dat (2D)
;
; Exercise the Mechanical Ball integration to the Measure Logic with PFC.
; Porosity is computed using measurement circles with random radii.
; Results for different model resolutions are compared against closed
; form solutions.
;
; Itasca Consulting Group, Inc.
;
; ========================================================================
fish def define_meas_parameters(nmeas, pcntg)
; ARGUMENTS:
;  nmeas : integer value that indicates the number of randomlyplaced
; probing measurement circles that will be created
; Vary this parameter to see the effects on the final result
;  pcntg : 'percentage' argument of the MEASURE function
;
global p = pcntg
global n = nmeas
end
fish def compute_meas_center(r)
;Randomly computes the coordinates of the measurement circle center
;ARGUMENTS:
;  r : radius of the measurement sphere
;
global xmin, xmax, ymin, ymax
local xcmin = xmin+r
local xcmax = xmaxr
local ycmin = ymin+r
local ycmax = ymaxr
global xc = xcmin + (xcmax  xcmin)*math.random.uniform
global yc = ycmin + (ycmax  ycmin)*math.random.uniform
end
fish define create_meas_regions
global measradius = array.create(n)
local maxmeasradius = math.min(0.99*(xmaxxmin),0.99*(ymaxymin))
loop ii(1,n)
measradius(ii)= math.max(maxmeasradius/15, ...
math.random.uniform*(0.5*(maxmeasradius)))
end_loop
loop ii(1,n)
command
[compute_meas_center(measradius(ii))]
measure create id=[ii] positionx [xc] positiony [yc] ...
radius [measradius(ii)] tolerance [p]
end_command
end_loop
end
fish define check_error(val,t,average_radius)
; Compares porosities computed in the measurement circles
; with analytical values
; ARGUMENTS:
;  val: exact value of the porosity
;  t: name of the table with the final numerical results
;
global result = table.create(t)
local it = 0
local averageerror = 0
loop foreach local mp measure.list
it = it + 1
local res = measure.porosity(mp)
local diff = 100.0 * (res  val) / val
averageerror = averageerror + math.abs(diff)
table(t,math.abs(diff))=measradius(it)/average_radius
endloop
global check_error = averageerror/it
end
; ========================================================================
; EoF: MeasureLogic.dat (2D)
PolydisperseAssembly.dat (2D)
; fname: PolydisperseAssembly.dat (2D)
; Build a polydisperse granular assembly and compute porosity
; using the measure logic in PFC
; ========================================================================
;
model new
program call 'MeasureLogic'
fish def poly_assembly(size)
; Create a polydisperse packing in a square domain
; ARGUMENTS:
;  size : size of the domain
;
global xmin = 0.0
global ymin = 0.0
global xmax = size
global ymax = size
global width = xmaxxmin ; width of box
global height = ymaxymin ; height of box
global tot_surface = (width)*(height)
local num = 3500 ; number of particles
local rlo = size/200
local rhi = 1.25*rlo
local oversize = rhi
command
model domain extent ([0.0oversize],[width+oversize]) ...
([0.0oversize],[height+oversize])
wall generate box (0.0, [width]) (0.0, [height])
end_command
; generate the balls and give them their properties
command
ball generate radius [rlo] [rhi] gauss ...
box ([width*0.01], [width*0.99]) ...
([height*0.01], [height*0.99]) ...
id=1,[num]
end_command
end
fish def averageradius
local rad = 0
loop foreach local bp ball.list
rad = rad + ball.radius(bp)
endloop
global averageradius = rad/ball.num
end
fish def get_poros ;computes the exact value of porosity
local sum = 0.0
loop foreach local bp ball.list
sum = sum + math.pi * ball.radius(bp)^2
end_loop
global get_poros = 1.0  sum / tot_surface
end
[poly_assembly(20.0)]
[get_poros]
[define_meas_parameters(5000, 0.00)]
[create_meas_regions]
fish def display_error(porosity, s)
global error = check_error(porosity,s,averageradius)
local oo = io.out(string.build("Average Error = %1%)",error))
end
[display_error(get_poros,'Polydisperse')]
model save 'PolydisperseAssembly'
; ========================================================================
; EoF: PolydisperseAssembly.dat (2D)
RegularAssembly.dat (2D)
; fname: RegularAssembly.dat (2D)
; Build a regular granular assembly and compute porosity
; using the measure logic in PFC
; ========================================================================
;
model new
program call 'MeasureLogic'
fish define regular_assembly(size, nres,itype)
; Create a regular packing of monosized balls in a square domain
; ARGUMENTS:
;  size : size of the domain
;  nres : resolution of the packing
;  itype: type of packing (0:simple cubic; #0: hexagonal)
;
global xmin = 0.0
global ymin = 0.0
global xmax = size
global ymax = size
global radius = size / (2.0*nres)
local numballs = nres^2
command
model domain extent [xmin] [xmax] [ymin] [ymax]
end_command
if itype = 0 then
command
ball generate rad [radius] number [numballs] cubic
end_command
else
command
ball generate rad [radius] number [numballs] hexagonal
end_command
endif
end
fish def display_error(porosity, tablename)
global error = check_error(porosity,tablename,radius)
local oo = io.out(string.build("Average Error = %1%)",error))
end
[define_meas_parameters(5000, 0.00)] ;define measurement parameters
; Simple Cubic packing  low resolution
[regular_assembly(18.0, 9, 0)]
[create_meas_regions]
[display_error(0.2146,'cubic_lores')]
model save 'cubic_lores'
;
; Simple Cubic packing  high resolution
ball delete
[regular_assembly(18.0, 30, 0)]
[display_error(0.2146,'cubic_hires')]
model save 'cubic_hires'
;
; Hexagonal packing  low resolution
ball delete
[regular_assembly(18.0, 9, 1)]
[display_error(0.0931,'hexa_lores')]
model save 'hexa_lores'
;
; Hexagonal packing  high resolution
ball delete
[regular_assembly(18.0, 30, 1)]
[display_error(0.0931,'hexa_hires')]
model save 'hexa_hires'
model save 'RegularAssembly'
; ========================================================================
; EoF: RegularAssembly.dat (2D)
Endnotes
Was this helpful? ...  Itasca Software © 2023, Itasca  Updated: Jan 18, 2024 