Area of the Intersection between a Sphere and a Cylindrical Plane

Abstract

A proper understanding of the porous structure of packed beds of spheres is imperative in the analysis and design of the processes involving fluid flow and heat and mass transfer. The radial variation in porosity is of specific interest. When the positions and sizes of the spheres are known, the radial variation in porosity can be determined using volume-based, area-based, or line-based approaches. Here, the focus is on the area-based methods which employ the intersections between the spheres and selected cylindrical planes to determine the radial variation in porosity, focusing specifically on the calculation of the area of the curved elliptic intersection between a sphere and a cylindrical plane. Using geometrical considerations, analytical integral expressions have been derived based on the axial direction, angular direction, or the radial direction as independent variables. The integral expressions cannot be integrated analytically and have been evaluated using approximations or numerical integration. However, only indirect validation of the calculation of the intersection area has been provided by comparing the radial porosity profiles obtained with experimental data. This study provides direct validation of the calculated area through refined numerical integration of the primary integral expressions and the evaluation of the area employing computer-aided design software.

1. Introduction

Cylindrical packed beds consisting of spherical particles are encountered in industrial applications such as in pebble bed reactors in the nuclear industry, in operations in chemical processes requiring interphase contact (for example, separation and heterogeneous catalysis), and practical applications such as packed columns for chromatography and regenerative heat exchangers [1,2,3]. A proper understanding of the characteristics of the porous structure of cylindrical packed beds is imperative in the analysis and design of the processes involving fluid flow, and heat and mass transfer [1]. It has been found that the container wall has a significant effect on the packing structure in the near-wall region, which affects the flow distribution and heat and mass transfer [2,4]. The most important characteristic is the porosity or void fraction and the radial variation in porosity is of specific interest [5]. The radial variation in porosity has been determined directly from experimental measurements using various approaches [2,5,6,7]. When the positions and sizes of the spheres are known, the radial variation in porosity can be determined using numerical techniques. The positions and diameters of the spheres have been determined experimentally [1,8,9] as well as numerically using packing algorithms [10,11] and the discrete element method [12,13]. The radial variation in porosity can be determined numerically using volume-based [8,14,15], area-based [16,17,18,19], or line-based [20,21,22] approaches. In this study, the focus is on the area-based methods of Mariani et al. [16], Mueller [18], and Du Toit [17,19], which employ the intersection between the spheres and selected cylindrical planes to determine the radial variation in porosity. The focus is specifically on the calculation of the area of the curved elliptic intersection between a sphere and a cylindrical plane.

Using geometrical considerations, Mariani et al. [16] derived an analytical integral expression in terms of the axial direction based on the analytical expression describing the lengths of the in-plane arcs as a function of axial position, defining the intersection. The integral cannot be evaluated analytically and was transformed to expressions involving Legendre complete elliptic integrals of the first and second kind. Mariani et al. [16] evaluated the elliptic integrals using an efficient numerical procedure to calculate the area. Du Toit [19], following a similar approach, calculated the area by integrating the derived analytical integral expression numerically.

Using the same geometrical considerations, Du Toit [17,19] derived an integral expression in terms of the circumferential or angular direction based on the analytical expression describing the in-plane vertical or axial lines, as function of angular position, defining the intersection. Since the integral cannot be evaluated analytically, Du Toit [17,19] integrated the analytical integral expression numerically to determine the area.

Mueller [18] derived an integral expression in terms of the radial direction based on the analytical expression for the axial height of the area, as a function of the radial position, obtained from the equations describing the surfaces of the sphere and the cylindrical plane. This integral can also not be evaluated analytically, and Mueller [18] therefore evaluated it numerically to determine the intersection area.

Mariani et al. [16], Du Toit [17,19], and Mueller [18] only provided indirect validation of the calculation of the intersection area by comparing the radial porosity profiles they obtained with available corresponding experimental results or by calculating the number of spheres in the bed. This study, in the first place, provides direct validation of the calculation of the intersection area through the refined numerical integration employing Simpson’s rule [23] of the primary integral expressions of Mariani et al. [16], Du Toit [17,19], and Mueller [18] and the evaluation of the area employing computer-aided design software. Secondly, the characteristics and the performance of the respective approaches are compared in terms of the number of uniformly spaced integration points that are required to obtain an accurate result. Four test cases representing typical sphere—cylindrical plane configurations that can be encountered are used to perform the analysis.

2. Theoretical Overview

2.1. Intersection of Cylindrical Plane and Sphere

A proper understanding of the radial variation in the porosity in cylindrical packed beds is important due to the fact that the container wall has a significant effect on the porous structure in the near-wall region which affects the heat and mass transfer and the flow distribution [1,2,4]. When the positions and the diameters of the spheres are known, the radial variation in porosity can be obtained employing numerical techniques [14,15,16,17,18,19,20,21,22]. In the area-based methodologies of Mariani et al. [16], Mueller [18], and Du Toit [17,19], the radial variation in porosity is determined by considering the intersections between the spheres and selected cylindrical planes. The accurate calculation of the area of the curved elliptical intersection between a sphere and a cylindrical plane is critical for the success of these methodologies.

If we define the radius of the sphere to be $r_p$, the radial position of the sphere to be $r_s$, and the radial position of the cylindrical plane to be $r$, the cylindrical plane intersects the sphere; that is, it cuts through the sphere or is inside the sphere when:

$$\max \left(0,r_s-r_p\right)\le r\le r_s+r_p.$$

(1)

Figure 1a shows a top view of the case where the cylindrical plane cuts through the sphere and Figure 1b shows a side view depicting the resulting curved elliptical intersection surface.

The cylindrical plane cuts through the sphere when:

$$r_p\le r_s\le R_c-r_p \mathrm{or} 0<r_s<r_p\hspace{1em}\mathrm{and}\hspace{1em}{r}_p-r_s<r,$$

(2)

where $R_c$ is the radius of the cylindrical container. The angular limits of the elliptical intersection surface are defined by the intersections of the cylindrical plane and the circumference of the center plane of the sphere at $-{\theta }_s$ and $+{\theta }_s$. It is assumed that the axial position of the sphere center is at $z=0$. The upper limit of the elliptical surface is at $z_s=z_B$.

Figure 2a shows a top view when the cylindrical surface is inside the sphere and Figure 2b shows the resulting intersection surface. The cylindrical plane is inside the sphere when:

$$0<r_s<r_p\hspace{1em}\mathrm{and}\hspace{1em}r\le r_p-r_s.$$

(3)

The intersection plane lies between $-\pi \le \theta \le +\pi$ or as shown in Figure 2b between $0\le \theta \le 2\pi$. The lower limit of the curve defining the upper edge of the intersection surface is at $z_A$. In the case of Figure 1b we have that $z_A=0$.

2.2. Du Toit Angular Integration

Du Toit [17,19] has shown the area $A$ of the intersection surface can be obtained by integrating in the angular or tangential $\left(\theta \right)$ direction as indicated in Figure 1b. The area is given as:

$$A=4{\displaystyle {\int }_0^{{\theta }_s}{z}_{\theta }}rd\theta ,$$

(4)

where $z_{\theta }$ is axial height of the elliptical curve defining the upper edge of the intersection surface relative to $z=0$. The factor 4 accounts for the symmetry of the intersection surface around $z=0$ and $\theta =0$. The upper integration limit ${\theta }_s$ (intersection angle) is given as:

$${\theta }_s={\cos }^{-1}\left[\frac{r^2-r_p^2+r_s^2}{2 r r_s}\right],$$

(5)

when the cylindrical plane cuts through the sphere. In the case where the cylindrical plane remains inside the sphere the intersection angle is given as:

$${\theta }_s=\pi .$$

(6)

The axial height $z_{\theta }$ of the upper edge of the intersection surface at the angular position $\theta$ can be obtained as:

$$z_{\theta }=\sqrt{r_p^2-r^2-r_s^2+2 r r_s\cos \left(\theta \right)}.$$

(7)

The integral Equation (4) cannot be calculated analytically, and Du Toit [17,19] therefore integrated Equation (4) numerically using Simpson’s rule [23].

2.3. Du Toit Axial Integration

Du Toit [19] has shown the area $A$ of the intersection surface can also be obtained by integrating in the axial $\left(z\right)$ direction. The area is given as:

$$A=2{\displaystyle {\int }_0^{z_s}{S}_z}dz,$$

(8)

where $S_z$ is the in-plane arc at the axial position $z$ defining the intersection surface. The factor 2 accounts for the symmetry of the intersection surface around $z=0$. The upper integration limit $z_s$ (axial integration height) is given as:

$$z_s=\sqrt{r_p^2-{\left(r-r_s\right)}^2}.$$

(9)

The length of the in-plane arc $S_z$ is obtained as:

$$S_z=2{\theta }_{z}r,$$

(10)

where ${\theta }_z$ is the in-plane arc angle. When the cylindrical plane cuts through the sphere the in-plane arc angle ${\theta }_z$ can be obtained from:

$${\theta }_z={\cos }^{-1}\left[\frac{r_s^2+r^2+z^2-r_p^2}{2 r r_s}\right],$$

(11)

whilst when the cylindrical plane remains within the sphere, the in-plane arc angle is given as:

$${\theta }_z=\pi$$

(12)

.

The integral Equation (8) can also not be calculated analytically, and Du Toit [19] therefore also integrated Equation (8) numerically using Simpson’s rule [23].

2.4. Mariani Elliptical Integration

Mariani et al. [16] also adopted the axial integration approach discussed in Section 2.3, but chose to write the integral expression to calculate the area of the intersection plane as:

$$A=4r\left[\pi z_A+{\displaystyle {\int }_{z_A}^{z_B}{\cos }^{-1}\left(\frac{r_s^2+r^2+z^2-r_p^2}{2r{r}_s}\right)dz}\right],$$

(13)

where the upper integration limit is $z_B=z_s$ from Equation (9). The lower integration limit $z_A$ is given as:

$$z_A=\{\begin{array}{ll}\sqrt{r_p^2-{\left(r+r_s\right)}^2}\hfill & \mathrm{when} r<r_p-r_s\hfill \\ 0\hfill & \mathrm{when} r>r_p-r_s.\hfill \end{array}$$

(14)

The first condition in Equation (14) occurs when the cylindrical plane remains inside the sphere as shown in Figure 2 and the second condition occurs when the cylindrical plane cuts through the sphere as shown in Figure 1. To facilitate the integration of Equation (13) Mariani et al. [16] have shown that the integral expression can be transformed to:

$$A=(\begin{array}{ll}8r{z}_{B}E\left(k\right)\hfill & \mathrm{when} k\le 1\hfill \\ 8r{z}_B\left[\left(k^{-1}+k\right)K\left(k^{-1}\right)+kE\left(k^{-1}\right)\right]\hfill & \mathrm{when} k>1,\hfill \end{array}$$

(15)

where the factor $k$ is defined as:

$$k=\frac{2\sqrt{r{r}_s}}{z_B}.$$

(16)

Substituting Equation (9) in Equation (16) it can be shown that when Equation (2) is valid, then $k>1$ and when Equation (3) is valid, then $k\le 1$. $K\left(m\right)$ is the Legendre complete elliptic integral of the first kind and is defined as:

$$K\left(m\right)={\displaystyle {\int }_0^{\pi /2}\frac{d\theta }{\sqrt{1-m^2{\sin }^2\left(\theta \right)}}}\hspace{1em}\mathrm{for}\hspace{1em}0\le m\le 1,$$

(17)

whilst $E\left(m\right)$ is the Legendre complete elliptic integral of the second kind and is defined as:

$$E\left(m\right)={\displaystyle {\int }_0^{\pi /2}\sqrt{1-m^2{\sin }^2\left(\theta \right)} d\theta }\hspace{1em}\mathrm{for}\hspace{1em}0\le m\le 1$$

(18)

.

The Legendre complete elliptic integrals cannot be calculated analytically, and the values must be obtained from tables, series expansions, or by employing numerical integration [24,25,26]. Note the fact that $K\left(k\right)$ has a singular value at $k=k^{-1}=1$ does not affect the application of Equation (15). Mariani et al. [16] integrated the Legendre complete elliptic integrals using an efficient numerical procedure.

2.5. Mueller Radial Integration

Mueller [18] derived an integral expression in terms of the radial direction $x$ based on the analytical expression for the axial height $z_x$, at the radial position $x$, of the upper edge of the intersection surface above $z=0$ obtained from the equations describing the surfaces of the sphere and the cylindrical plane. The equation describing the surface of the sphere is given as:

$${\left(x-r_s\right)}^2+y^2+z^2=r_p^2,$$

(19)

whilst the equation for the cylindrical plane is given as:

$$x^2+y^2=r^2.$$

(20)

Substituting Equation (20) in Equation (19) gives the expression for the axial height $z_x$ of the edge of the intersection surface:

$$z_x=\sqrt{r_p^2-r^2+2x{r}_s-r_s^2}.$$

(21)

The in-plane arc length element $ds$ can be written, applying Equation (20), as:

$$ds=\sqrt{d{x}^2+d{y}^2}=\sqrt{1+{\left(\frac{dy}{dx}\right)}^2} dx=\frac{rdx}{\sqrt{r^2-x^2}}.$$

(22)

The integral expression to calculate the area of the intersection surface can then be written as:

$$A=C{\displaystyle {\int }_{LL}^r\left(\frac{z_x}{\sqrt{r^2-x^2}} r\right)dx},$$

(23)

where the lower integration limit $LL$ and the integration constant $C$ are defined as:

$$LL=\frac{r_s^2+r^2-r_p^2}{2r_s}\hspace{1em}\mathrm{and}\hspace{1em}C=4,$$

(24)

when the cylindrical plane cuts through the sphere and the conditions in Equation (2) are applicable. When the cylindrical plane remains within the sphere and the condition in Equation (3) are applicable, then the lower integration limit and the integration constant are:

$$LL=-r\hspace{1em}\mathrm{and}\hspace{1em}C=4$$

(25)

Lastly, when $r_s=0$ and $0\le r\le r_p$, then the lower integration limit and the integration constant are defined as:

$$LL=0\hspace{1em}\mathrm{and}\hspace{1em}C=8.$$

(26)

Note that $x=-r$ and $x=r$ are singular points. The integral Equation (23) cannot be calculated analytically and therefore needs to be evaluated using numerical integration or another suitable method. Mueller [18] does not describe the approach employed to integrate the integral expression numerically.

3. Results

3.1. Calculation of Intersection Areas

Four sphere–cylindrical plane test configurations were selected to evaluate the ability of the methodologies of Mariani et al. [16], Du Toit [17,19], and Mueller [18] to calculate the areas of the intersections between the spheres and the cylindrical planes accurately. The four test configurations were chosen to be representative of sphere–cylindrical plane configurations that typically occur in cylindrical packed beds consisting of spheres for which the radial variation in porosity must be obtained.

In this study, Equations (4), (8), (15), and (23) were evaluated numerically employing Simpson’s rule [23]. It was not the purpose of the study to compare the performance of different numerical integration methods. Simpson’s rule was, therefore, selected for the numerical integration because all the functions to be integrated are smooth, as well as the fact that Simpson’s rule is third-order accurate and simple to implement.

The first step in the analyses was to find the maximum integration interval size in each case for which the value obtained for the intersection area changed by less than $1\times {10}^{-4}{ \mathrm{mm}}^2$ between successive values selected for the integration interval size.

The second step in the analyses was to compare the intersection areas that were obtained for the different cases with the corresponding areas obtained using the finite element code COMSOL Multiphysics [27] and the computer aided design (CAD) codes SOLIDWORKS [28] and NX [29].

Finally, the performance of the methodologies was compared considering the number of integration points required in each case to obtain an accurate solution for the intersection area.

3.1.1. Test Configurations

The test configurations that were selected for the analyses are summarized in

Table 1. Test configurations for the analyses.

Parameter	Case 1c	Case 2c	Case 3c	Case 4c
$r_p \left(\mathrm{mm}\right)$	30.0	30.0	30.0	30.0
$r_s \left(\mathrm{mm}\right)$	7.5	15.0	22.5	60.0
$r \left(\mathrm{mm}\right)$	15.0	7.5	37.5	60.0
${\theta }_s \left(\mathrm{deg}\right)$	180.0	180.0	53.13010	28.95502
$z_s=z_B \left(\mathrm{mm}\right)$	29.04738	29.04738	25.98076	30.0
$z_A \left(\mathrm{mm}\right)$	19.84313	19.84313	0.0	0.0
$k$ (-)	0.730296743	0.730296743	2.236068	4.0
$LL \left(\mathrm{mm}\right)$	−15.0	−7.5	22.5	52.5

and shown in Figure 3.

For Cases 1c and 2c the cylindrical plane remained within the sphere with $r>r_s$ for Case 1c and $r<r_s$ for Case 2c. For Cases 3c and 4c the cylindrical plane cuts through the sphere with $r_s<r_p$ and $r>r_p-r_s$ for Case 3c and $r_s>r_p$ for Case 4c. The same sphere radius of 30 mm was used for all the cases reported in this study. In Table 1 the values of the limits and parameters ${\theta }_s$, $z_s=z_B$, $z_A$, $LL$, and $k$ occurring in Equations (4), (8), (15), and (23) are also tabulated.

The red lines in the CAD drawings of Cases 1c, 2c, 3c, and 4c shown in Figure 3 indicate the radial distance of the cylindrical plane from the centre line, the radial position of the centre of the sphere from the centre line and the radius of the sphere respectively.

3.1.2. Du Toit Angular Integration

To integrate Equation (4) numerically using Simpson’s rule [23], to obtain the intersection area, it is rewritten to become:

$$A=\frac{4{\displaystyle \sum \left(z_{\theta ,i-1}+4z_{\theta ,i}+z_{\theta ,i+1}\right)r\Delta \theta }}{3}\hspace{1em}\mathrm{for}\hspace{1em}i=2,4,6\dots ,n-1,$$

(27)

where

$$z_{\theta ,i}=\sqrt{r_p^2-r^2-r_s^2+2 r r_s\cos \left({\theta }_i\right)}\hspace{1em}\mathrm{for}\hspace{1em}0\le {\theta }_i\le {\theta }_s,$$

(28)

and $\Delta \theta$ the angular increment. $n$ is the number of angular integration points with $n\ge 3$ and an uneven number. The results for the analysis of the numerical integration of Equation (4) are summarized in

Table 2. Convergence analysis for Equation (4) [17,19].

	Case 1c	Case 2c	Case 3c	Case 4c
$\Delta \theta \left(\mathrm{deg}\right)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$
1.00 × 10−1	4648.7417	2324.3708	2811.8029	2849.9743
1.00 × 10−2	4648.7417	2324.3708	2811.8345	2850.0553
1.00 × 10−3	4648.7417	2324.3708	2811.8355	2850.0579
1.00 × 10−4			2811.8355	2850.0580
1.00 × 10−5				2850.0580

. The values for the angular increment $\Delta \theta$ are the nominal values in degrees that were specified. The code adapts the angular increment where necessary to obtain an even number of increments for the interval $0\le \theta \le {\theta }_s$ in order to apply Simpson’s rule. Note that in Equation (27) the curvature of the elliptic surface is represented exactly, whilst the edge of the surface is represented discretely.

It can be seen in Table 2 that an angular increment of $\Delta \theta =1\times {10}^{-3}$ deg can be considered as sufficient to obtain a converged solution for the integral.

3.1.3. Du Toit Axial Integration

To integrate Equation (8) numerically using Simpson’s rule [23], to obtain the intersection area, it is rewritten to become:

$$A=\frac{2{\displaystyle \sum \left(S_{z,i-1}+4S_{z,i}+S_{z,i+1}\right)\Delta z}}{3}\hspace{1em}\mathrm{for}\hspace{1em}i=2,4,6\dots ,n-1,$$

(29)

where

$$S_{z,i}=2{\theta }_{s,i}r,$$

(30)

with

$${\theta }_{z,i}={\cos }^{-1}\left[\frac{r_s^2+r^2+z_i^2-r_p^2}{2 r r_s}\right]\hspace{1em}\mathrm{for}\hspace{1em}0\le z_i\le z_s,$$

(31)

and $\Delta z$ the axial increment. $n$ is the number of axial integration points with $n\ge 3$ and an uneven number. The results for the analysis of the numerical integration of Equation (8) are summarized in

Table 3. Convergence analysis for Equation (8) [19].

	Case 1c	Case 2c	Case 3c	Case 4c
$\Delta z \left(1/d_p\right)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$
1.00 × 10−1	4546.4251	2273.2126	2761.4722	2812.8134
1.00 × 10−2	4645.6167	2322.8083	2810.4152	2848.8880
1.00 × 10−3	4648.6216	2324.3108	2811.7911	2850.0210
1.00 × 10−4	4648.7411	2324.3705	2811.8341	2850.0568
1.00 × 10−5	4648.7416	2324.3708	2811.8355	2850.0579
1.00 × 10−6	4648.7417	2324.3708	2811.8355	2850.0580
1.00 × 10−7	4648.7417			2850.0580

. The values for the axial increment $\Delta z$ are the nominal values as a fraction of the sphere diameter that were specified. The code also adapts the axial increment where necessary to obtain an even number of increments for the interval $0\le z\le z_s$ in order to apply Simpson’s rule. Note that in Equation (29), the curvature of the elliptic surface is represented exactly, whilst the edge of the surface is represented discretely.

It can been in Table 3 that an axial increment of $\Delta z=1\times {10}^{-6} (1/d_p)$ is required to obtain a converged solution for the integral. The fact that the numerical integration for Equation (8) requires a finer increment than the numerical integration of Equation (4) to obtain a converged solution can be attributed to the intersection surface being flatter at the top $\left(+z_s\right)$ and bottom $\left(-z_s\right)$ than at the left hand $\left(+{\theta }_s\right)$ and right hand $\left(-{\theta }_s\right)$ sides of the elliptical intersection plane. It should also be noted that the angular increment of $\Delta \theta =1\times {10}^{-3}$ deg in the case of the angular integration translates to $9\times {10}^{-6}<r\Delta \theta /d_p<3\times {10}^{-5}$ for the four cases considered.

3.1.4. Mariani Elliptic Integration

The ability of Equation (15) to calculate the intersection area accurately is dependent on the accurate numerical integration of the Legendre complete integrals Equations (17) and (18). Using Simpson’s rule [23], Equation (17) can be rewritten to become:

$$K\left(m\right)=\frac{{\displaystyle \sum \left(I_{K,i-1}+4I_{K,i}+I_{K,i+1}\right)\Delta \theta }}{3}\hspace{1em}\mathrm{for}\hspace{1em}i=2,4,6\dots ,n-1,$$

(32)

where

$$I_{K,i}=\frac{1}{\sqrt{1- m^2{\cos }^2\left({\theta }_i\right)}}\hspace{1em}\mathrm{for}\hspace{1em}0\le {\theta }_i\le \pi /s,$$

(33)

and $\Delta \theta$ the angular increment. $n$ is the number of angular integration points with $n\ge 3$ and an uneven number. Using Simpson’s rule [23], Equation (18) can be rewritten to become:

$$E\left(m\right)=\frac{{\displaystyle \sum \left(I_{E,i-1}+4I_{E,i}+I_{E,i+1}\right)\Delta \theta }}{3}\hspace{1em}\mathrm{for}\hspace{1em}i=2,4,6\dots ,n-1,$$

(34)

where

$$I_{E,i}=\sqrt{1- m^2{\sin }^2\left({\theta }_i\right)}\hspace{1em}\mathrm{for}\hspace{1em}0\le {\theta }_i\le \pi /s,$$

(35)

and $\Delta \theta$ the angular increment. $n$ is the number of angular integration points with $n\ge 3$ and an uneven number. The results for the analysis for the numerical integration of Equations (17) and (18) are summarized in

Table 4. Validation of numerical integration of Legendre complete elliptic integrals.

	Simpson’s Rule		Milne-Thomson [30]
$m^2$	$E(m^2)$	$K(m^2)$	$E(m^2)$	$K(m^2)$
0.1	1.530757636	1.612441348720	1.530757636	1.612441348720
0.5	1.350643881	1.854074677301	1.350643881	1.854074677301
0.9	1.104774733	2.578092113348	1.104774733	2.578092113348

and compared with the corresponding values tabulated by Milne-Thomson [30]. It can be seen in Table 4 that the values obtained using Simpson’s rule are in exact agreement with corresponding values given by Milne-Thomson [30]. The results were obtained for integration increment of $\Delta \theta =1\times {10}^{-2}$ rad.

The results for the analysis for the numerical integration of Equation (15) are summarized in

Table 5. Convergence analysis for Equation (15) [16].

	Case 1c	Case 2c	Case 3c	Case 4c
$\Delta \theta \left(\mathrm{rad}\right)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$
1.00 × 10−1	4648.7417	2324.3708	2811.8355	2850.0580
1.00 × 10−2	4648.7417	2324.3708	2811.8355	2850.0580
1.00 × 10−3	4648.7417	2324.3708	2811.8355	2850.0580

. The values for the independent variable increment $\Delta \theta$ are the nominal values in radians that were specified. The code adapts the independent variable increment where necessary to obtain an even number of increments for the interval $0\le \theta \le \pi /2$ in order to apply Simpson’s rule.

It can been in Table 5 that an independent variable increment of $\Delta \theta =1\times {10}^{-1}$ rad is sufficient to obtain a converged solution of the integral. Note that in Equation (15), the curvature of the elliptic surface and the edge of the surface are represented exactly.

3.1.5. Mueller Radial Integration

To integrate Equation (23) numerically using Simpson’s rule [23], to obtain the intersection area, it is rewritten to become:

$$A=\frac{C{\displaystyle \sum \left(I_{M,i-1}+4I_{M,i}+I_{M,i+1}\right)r\Delta z}}{3}\hspace{1em}\mathrm{for}\hspace{1em}i=2,4,6\dots ,n-1,$$

(36)

where

$$I_{M,i}=\frac{\sqrt{r_p^2-r^2+2 x r_s-r_s^2}}{\sqrt{r^2-x^2}}\hspace{1em}\mathrm{for}\hspace{1em}LL\le x_i\le r,$$

(37)

and $\Delta x$ the radial increment. $n$ is the number of radial integration points with $n\ge 3$ and an uneven number. The singularities at $x=-r$ and $x=+r$ were circumvented by setting the value of the integrand for $x_1$ equal to the value of the integrand for $x_2$ and by setting the value of the integrand for $x_n$ equal to the value of the integrand for $x_{n-1}$.

The results for the convergence analysis for the numerical integration of Equation (23) are summarized in

Table 6. Convergence analysis for Equation (23) [18].

	Case 1c	Case 2c	Case 3c	Case 4c
$\Delta x \left(1/d_p\right)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$
1.00 × 10−1	3357.7813	1492.4812	1837.1173	1828.8018
1.00 × 10−2	4266.0634	2049.1817	2485.3647	2375.8312
1.00 × 10−3	4527.3808	2238.5827	2709.9127	2700.6604
1.00 × 10−4	4610.3532	2297.2269	2779.5819	2802.9507
1.00 × 10−5	4636.6018	2315.7867	2801.6353	2835.1597
1.00 × 10−6	4644.9027	2321.6563	2808.6099	2845.3467
1.00 × 10−7	4647.5277	2323.5124	2810.8155	2848.5681
1.00 × 10−8	4648.3578	2324.0994	2811.5129	2849.5868
1.00 × 10−9	4648.6203	2324.2850	2811.7335	2849.9090
3.00 × 10−10	4648.6752	2324.3238	2811.7899	2849.9764

.

The values for the radial increment $\Delta x$ are the nominal values as a fraction of the sphere diameter that that were specified. The code adapts the independent variable increment where necessary to obtain an even number of increments for the interval $LL\le x\le r$ in order to apply Simpson’s rule. It can be seen in Table 6 that the solutions cannot be considered as fully converged for the smallest radial increment. Note that in Equation (23), the curvature of the elliptic surface, as well as the edge of the surface, are represented discretely.

3.2. Comparison of Areas and Performance

Table 7. Summary of results for intersection areas.

	Case 1c	Case 2c	Case 3c	Case 4c
	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$	$A{ (\mathrm{mm}}^2)$
Equation (4) [17]	4648.7417	2324.3708	2811.8355	2850.0580
Equation (8) [19]	4648.7417	2324.3708	2811.8355	2850.0580
Equation (15) [16]	4648.7417	2324.3708	2811.8355	2850.0580
Equation (23) [18]	4648.6752	2324.3238	2811.7899	2849.9764
COMSOL [27]	4648.7410	2324.3619	2811.8342	2850.0255
SOLIDWORKS [28]	4648.5749	2324.3166	2811.5967	2849.9363
NX [29]	4647.0269	2323.6438	2809.5591	2848.5306

is a summary of the converged numerical solutions for the areas of the intersection surfaces for cases 1c to 4c obtained using Equations (4), (8), (15), and (23) obtained from Table 2, Table 3, Table 5, and Table 6. Included in Table 7 also are the areas of the intersection surfaces obtained using the Measuring Geometry Tool in COMSOL Multiphysics [27], and the CAD packages SOLIDWORKS [28] and NX [29].

In Table 7 it can be observed that the agreement between the results obtained by Du Toit [17] (Equation (4)), Du Toit [19] (Equation (8)), Mariani et al. [16] (Equation (15)), Mueller [18] (Equation (23)), and COMSOL [27], SolidWorks [28], and NX [29] is very good. The maximum relative difference between the numerical results and the COMSOL results is 0.001%, the maximum difference between the numerical results and the SOLIDWORKS results is 0.008%, and the maximum difference between the numerical results and the NX results is 0.081%. It can, therefore, be concluded that the methodologies of Mariani et al. [16], Du Toit [17,19], and Mueller [18] calculated the area of the intersection between a sphere and a cylindrical plane correctly.

It has been noted that in the numerical integration of Equations (4) and (8) the curvature of the curved elliptical surface is represented exactly, but the edge of the surface discretely, whilst in the numerical integration of Equation (23) both the curvature of the curved elliptic surface and the edge of the surface are represented discretely. Compared to that in the numerical integration of Equation (15), both the curvature of the curved elliptic surface and the edge of the surface are represented exactly. The values of the respective integration increments are a reflection of these characteristics. Due to the differences in the physical meaning of the integration increments, they cannot be directly used as a measure of the computational effort, which is of interest for the practical implementation of the methodologies to determine the radial variation in porosity of a packed bed. A truer reflection of the computational effort is the number of integration points in each case required to obtain an accurate numerical solution.

Table 8. Summary of number of integration points.

	Case 1c	Case 2c	Case 3c	Case 4c
Equation (15) [16]	159	159	159	159
Equation (4) [17,19]	18,001	18,001	53,131	28,957
Equation (8) [19]	96,825	96,825	86,603	100,001
Equation (23) [18]	1,666,666,667	833,333,335	1,250,000,001	416,666,667

gives a summary of the number of integration points required in each of the cases listed in Table 7.

It can be seen in Table 8 that the numerical integration of Equation (15) requires the least number of integration points followed by Equation (4), then Equation (8), and lastly Equation (23). It can further be observed in Table 8 that the number of integration points required in the case of Equation (15) is independent of the sphere–cylindrical plane configuration. The reason for this is the fact that the numerical integration of Equation (15) is only dependent on the accurate numerical integration of the Legendre complete elliptic integral Equations (17) and (18). It can, therefore, be concluded that the methodology of Mariani et al. [16] is the most effective approach as measured by the number of integrations points that are required to calculate the area of the intersection between a sphere and a cylindrical plane accurately.

4. Conclusions

Cylindrical packed beds consisting of spherical particles are found in various industrial and practical applications [1,2,3]. The container wall has a significant influence on the packing structure in the near-wall region affecting the flow distribution and heat and mass transfer [2,4]. This requires an understanding of the characteristics of the radial variation in in the porosity [1,5]. When the positions and the diameters of the spheres forming the packed bed are known, numerical techniques [14,15,16,17,18,19,20,21,22] can be used to obtain the radial variation in porosity. Employing area-based methodologies, Mariani et al. [16], Mueller [18], and Du Toit [17,19] determine the radial variation in porosity by considering the areas of the intersections between the spheres in the packed bed and selected cylindrical planes. The accurate calculation of the area of the curved elliptical intersection between a sphere and a cylindrical plane is of critical importance for the success of these methodologies.

This study endeavoured to provide a direct validation of the calculation of the intersection through the refined numerical integral using Simpson’s Rule [23] of the primary integral expressions of Equation (4) [17,19], Equation (8) [19], Equation (15) [16], and Equation (23) [18]. Four representative sphere–cylindrical plane configurations were chosen to evaluate the ability of the methodologies to calculate the intersection area. The first step in the analyses was to find the maximum size of the integration interval for each case that gave a converged value for the intersection area. The second step in the analyses was to compare the intersection areas that were obtained for the different cases with the corresponding areas obtained using the finite element code COMSOL Multiphysics [27] and the computer aided design (CAD) codes SOLIDWORKS [28] and NX [29]. It was found that the corresponding intersection areas obtained by the methodologies of Mariani et al. [16], Du Toit [17,19], and Mueller [18] are in excellent agreement and also in very good agreement with the corresponding areas obtained using the Measuring Geometry Tool in COMSOL Multiphysics, and the CAD packages SOLIDWORKS and NX.

The study also considered the performance of the methodologies of Mariani et al. [16], Du Toit [17,19], and Mueller [18] by comparing the number of integration points required in each case to obtain an accurate solution for the intersection area. The number of integration points is considered to be a representative reflection of the computational effort and of interest for the practical implementation of the methodologies. It was found that the numerical integration of Equation (15) requires the least number of integration points followed by Equation (4), then Equation (8), and lastly Equation (23). It was further observed that the number of integration points required in the case of Equation (15) is independent of the sphere–cylindrical plane configuration.

It can thus be stated that in this study the calculation of the area of the intersection surface between a cylindrical plane and a sphere using the approaches of Mariani et al. [16], Du Toit [17,19], and Mueller [18] have been validated. It can also be concluded that the procedure of Mariani et al. [16], compared to the methodologies of Du Toit [17,19] and Mueller [18], requires the least computational effort to obtain an accurate solution for the intersection area based on Simpson’s rule [23], which was employed to perform the numerical integration of the relevant integral expressions.

As a natural extension of the study, it is recommended that the computational effort of more advanced numerical integration methods [23,31,32] be evaluated in view of the fact that it has been shown that the methodologies of Mariani et al. [16], Du Toit [17,19], and Mueller [19] give the correct results for the area of the intersection between a sphere and a cylindrical plane. However, it is recommended that this be done in the context of the calculation of the radial variation in porosity for a selection of cylindrical packed beds consisting of varying numbers of spheres. The study should include the pebble bed model consisting of 450,000 spheres generated by Suikkanen et al. [13].