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 mono-sized balls organized in crystal-like 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 crystal-like packings are considered. Both packings consist of equal-radius 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 two-dimensional or three-dimensional 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

\[n = 1 - D = 1 - \frac{A_s}{A}\]

For the simple cubic packing, we have, in two dimensions:

(1)\[n = 1 - \frac{\pi R^2}{(2R)^2}=1-\frac{\pi}{4} \approx 0.2146\]

and, in three dimensions:

(2)\[n = 1 - \frac{4/3 \pi R^3}{(2R)^3}=1-\frac{\pi}{6} \approx 0.4764\]

For the hexagonal packing, in two dimensions, we have:

(3)\[n = 1 - \frac{\pi R^2}{(\sqrt{3}R)(2R)}=1-\frac{\pi}{2 \sqrt{3}} \approx 0.0931\]

and, in three dimensions:

(4)\[n = 1 - \frac{8 \pi R^3}{(6\sqrt{3}R^2)(4R\sqrt{2/3})}=1-\frac{\pi}{3 \sqrt{2}} \approx 0.2595\]

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.

(5)\[n = 1-\frac{S_{balls}}{S} \quad (2D) \quad \quad n = 1-\frac{V_{balls}}{V} \quad (3D)\]

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 post-processing 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.

Figure 1: Cubic regular packing, low resolution.

Figure 2: Hexagonal regular packing, low resolution.

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}.

Figure 3: Polydisperse packing, low resolution.

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.

Figure 4: Cubic packing and measurement regions.

Figure 5: Polydisperse packing and measurement regions.

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:

Table 1: Porosity Comparison Regular Assemblies - 2D Model

	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:

Table 2: Porosity Comparison Regular Assemblies - 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.

Figure 6: Porosity measurement, regular assemblies - average error [% - log] vs. circle diameter.

We do the same comparison in the case of the polydisperse assembly. The error is plotted in Fig. 7.

Figure 7: Porosity measurement, polydisperse assembly - average error [% - log] vs. ircle diameter.

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

[Deresiewicz1958a]

Deresiewicz, H. “Mechanics of Granular Matter”, in Advances in Applied Mechanics, Vol. 5, pp. 233-306. 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 randomly-placed
;             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 = xmax-r
    local ycmin = ymin+r
    local ycmax = ymax-r
    local zcmin = zmin+r
    local zcmax = zmax-r
    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*(zmax-zmin)), ...
                        math.min((0.99*(xmax-xmin)), ...
                        (0.99*(ymax-ymin))))
  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 pos-x [xc] pos-y [yc] pos-z [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      = xmax-xmin  ; width of box
  global height     = ymax-ymin  ; height of box
  global depth      = zmax-zmin  ; 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.0-oversize],[width+oversize])  ...
                        ([0.0-oversize],[height+oversize]) ...
                        ([0.0-oversize],[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 mono-sized 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 randomly-placed
;             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 = xmax-r
    local ycmin = ymin+r
    local ycmax = ymax-r
    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*(xmax-xmin),0.99*(ymax-ymin))
    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] position-x [xc] position-y [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      = xmax-xmin  ; width of box
  global height     = ymax-ymin  ; 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.0-oversize],[width+oversize]) ...
                        ([0.0-oversize],[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 mono-sized 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

[1]	To view this project in PFC, use the program menu. Help ▼ Examples…   ⮡   PFC     ⮡   Verifications       ⮡   MeasureLogic         ⮡   MeasurePorosity.prj