**2. Methodology**

#### *2.1. Materials and Manufacturing*

The large-format composite structures discussed in this research are made of a blend of recycled High-Density Polyethylene (HDPE), recycled Glass Fiber Polypropylene (GFPP), ADC (chemical blowing agent), and carbon black. This mixture consists of ~20 % GFPP, ~80% HDPE, ~0.5% ADC, and carbon black by mass. This blend was mixed, melted, and pressurized into a mold cavity. The exterior of the mold was cooled using room temperature water. The cooling affected the outer shell first and allowed the blowing agen<sup>t</sup> to react in the core forming a dense outer core and a foamed inner core, with varying cell size distribution and spatially varying density. Three samples from different crossties were sectioned using an industrial chop saw, and debris from the cutting process was removed. The cross-section of the crosstie was then photographed in a low light environment with the flash from the COTS camera, yielding an image that was approximately 3000 × 2300 pixels. The dimensions of the core for each sample were 0.19 m × 0.14 m (7.38 in × 5.38 in) as shown in Figure 2a.

**Figure 2.** Cross-section of the presented samples. (**a**) Images of sample 1, 2, and 3 from the commercial off-the-shelf (COTS) camera. (**b**) Binary plots of samples 1, 2, and 3.

#### *2.2. Image Processing to Generate Cells*

The photographed samples were imported into the MATLAB environment as an array representing the image color density as a function of pixel location. The uploaded images taken from three different cross-sections are shown in Figure 2a. Pixels were converted from a grayscale value to a binary value using color contrast differentiation. It was found that a yellow tint using the automatic enhance option included in the Windows photo editor provided reasonable pre-processing contrast variations in differentiating the solid region from the darker cells. A threshold value of 25% was selected for the images, and the noise within the binary image was removed using a morphological structure operation by ignoring individual false single pixels. As shown in Figure 2b, the white pixels denote cells, and the neighboring solid surface region is denoted in black.

This method then considers the effective area of the cells shown in white pixels and draws a circle around it. The dimensions of this circle are defined as the effective diameter using the automated MATLAB subroutine which is much faster compared to highly sophisticated and manual use of an optical microscope. During image processing, false positives must be identified and properly addressed. As seen in Figure 3a,b, it can be noted that in the picture taken with the COTS camera, there are color contrast variations within the cell cavity itself. The first false positive, termed a donut, is a cell region that contains black pixels resembling a solid shell. For the overall sample size of 7.38 × 5.38 inches, the image is approximately 3000 × 2300 pixels, and the size of each pixel is 63.5 μm. For the composite structure discussed, the radius of the fibers was measured to be 8.56 μm ± 0.57 μm [39] and thus the fibers are not visible within the image. Optical microscopy was performed on the surface of the cells, and the fibers did not protrude from the cell wall. A second false positive, termed a double circle, is a shell region that appears to cover a portion of the cell due to flaws introduced from the saw blade used in the sectioning process.

The first false positives were corrected using a subroutine that automatically takes into account the overall area by circles drawn around each cell, including the black pixels within the cell region, as shown in Figure 3c. The second false-positive requires an additional processing step whereby any captured features that reside within the domain of another feature, such as those highlighted in Figure 3c, were removed from consideration. If the double circles were not removed, the area of the cells within the region would be overestimated. Figure 3d shows the cell distribution requiring substantially less archived data as only the centroid and effective radius is stored and used in subsequent processing. This down-select process allows the homogenization algorithm to process more quickly and retains the effective features of the various cells for micromechanical predictions.

**Figure 3.** Image processing false positives. (**a**) Image with a yellow tint showing color contrasts within cells. (**b**) Pixel formation of cell donuts based on color contrast. (**c**) Cells with double circles. (**d**) Filtering the image of double circles.

It is important to note that the cross-sections shown are cut at a random plane along the longitudinal axis of the composite. As this cutting process may not necessarily pass through the centroid of the cell, it will tend to underestimate the actual cell spherical equivalent volume. To account for the underestimation, a correction factor is used to estimate the largest radius of the 3D spheroidal cell instead of the chord length as seen in a 2D cell. This correction factor is taken from the work of Pinto et al. [40] and the radius of each cell is adjusted as

$$r\_{3D} = 1.273 \, r\_{2D} \tag{1}$$

#### *2.3. 2D Area Validation*

To validate the cell area obtained using image thresholding and the COTS camera, a detailed study was conducted on the subregion of the COTS image shown in Figure 4a for sample 2. This region contains 36 identifiable cells of varying sizes, and each cell is compared to the area measured manually using the high-resolution 3D Keyence 3100 microscope. For example, a single cell was selected and is shown in Figure 4b. Its effective pixel diameter was measured using the image thresholding algorithm discussed above, and the effective area is compared to the area measured using the Keyence digital microscope that includes the non-circular perimeter of the cell in the area calculation. Each of the 36 cells were measured manually using the Keyence digital microscope and their regions are highlighted in Figure 4c. The areas from the aforementioned imaging algorithm, shown in Figure 4d, were then compared to each of manually identified areas from the high-resolution microscopy results.

**Figure 4.** Cell area validation using a Keyence digital microscope. (**a**) Measuring diameter as seen from the Keyence microscope. (**b**) Measuring the effective diameter of the same cell using the new proposed method. (**c**) Manually measuring the area of each cell using the Keyence microscope. (**d**) Automated measuring of the effective diameters of all the cells.

The total area for the region of one cell is highlighted in Figure 4b. The total area for the cells in the image was calculated with image thresholding approach using the image from the COTS camera as 170.69 mm2, and the total area for all 36 cells was measured using high-resolution microscopy was 166 mm2, yielding an error of 3%. The error was calculated using the equation,

$$Err = \frac{|A\_{Microsopy} - A\_{COTS\ Imging}|}{A\_{Microsopy}} \times 100\tag{2}$$

#### *2.4. Image Mirroring and Density Homogenization*

To calculate the homogenized density of the surface for use in the finite element model, a grid of 400 × 250 points was mapped onto the surface of each respective sectioned sample. The sample image was then mirrored about the edges to avoid biasing in any region near the surface and to allow continuous regions for homogenization at all points within the image. A region for a single point was depicted in Figure 5. A fixed pixel radius of *R* = 300 was taken for each of the 3000 × 2300 pixel images, and every cell was identified within the radius *R*.

**Figure 5.** Mirrored image preventing image biasing near the boundaries. A single grid point is highlighted, showing the region of consideration used for local density homogenization.

The cells identified within the circle were averaged to calculate the local density using an exponential decaying weight function. This allowed the areas of the cells closest to the grid point to have the highest contribution to the effective density, whereas the cells farther from the center have less of a contribution to the effective density. This is also denoted by the color variations, as seen in Figure 5. The cells closest to the center are colored red and provide the highest contribution to the density homogenization, whereas the cells farthest from the center are shown in colors that transition to yellow and then to blue as their contribution decreases. The exponential decay function is defined as,

$$\widetilde{\rho}(\mathbf{x},y) = \int\_0^{2\pi} \int\_0^{\infty} \Psi(\mathbf{r},\theta) \varepsilon^{-(r(\mathbf{x},y)/R)^2} r(\mathbf{x},y) dr d\theta \simeq \int\_0^{2\pi} \int\_0^{3R} \Psi(\mathbf{r},\theta) \varepsilon^{-(r(\mathbf{x},y)/R)^2} r(\mathbf{x},y) dr d\theta \tag{3}$$

where *r* is the local radius from the center of the region of radius *R*, *ρ*(*<sup>x</sup>*, *y*) is the localized density, and *ψ*(*<sup>r</sup>*, *θ*) is the probability of finding a cell at *r*, *θ*. The upper limit on the radial integral is approximated as being upper bounded by 3*R*, as this will take into account 99.9% of the area within the distribution while avoiding unnecessary calculations over all cells captured within the image. As observed, the cell sizes are much smaller than the region *R*. Therefore, for discrete data sets, Equation (3) can be approximated as,

$$\widetilde{\rho}(\mathbf{x}, \mathbf{y}) \simeq \sum\_{i=1}^{N} A\_i \mathbf{e}^{-(r(\mathbf{x}, \mathbf{y})/R)^2} \tag{4}$$

where *Ai* is the area of the *ith* cell and *N* is the total number of cells in the entire image. Again, Equation (4) as in Equation (3), will disregard cells that have a centroid more than 3*R* from the grid point being analyzed as their contribution will not significantly alter the final solution due to the pre-multiplication of the exponential decaying function. To normalize over the entire area of the image, the density *ρ*(*<sup>x</sup>*, *y*) is calculated as,

$$\rho(\mathbf{x}, y) = \frac{\tilde{\rho}(\mathbf{x}, y) \sum\_{i=1}^{N} A\_i}{\int \int \tilde{\rho}(\mathbf{x}, y) d\mathbf{x} dy} \tag{5}$$

where all the cell areas within the image are summed from *i* = 1 to *N* and *x* and *y* are the grid points selected.

#### *2.5. Material Properties Prediction*

Based on the previous research [39], material properties were predicted using only constitutive properties. This was accomplished using the inputs of raw materials, including as elastic moduli and Poisson's ratios of HDPE and GFPP, glass fiber aspect ratio, and foamed core density. Since modeling over each individual fiber in a 2.62 m long composite structure is computationally impractical, Advani and Tucker [41] proposed a model which uses the average orientation of the fibers within a region, termed the orientation tensors. Specifically of interest are the second order *aij* and the fourth order *aijkl* tensors, defined as,

$$a\_{ij} = \int\_{\mathbb{S}} p\_i p\_j \psi(\mathbf{p}) d\mathbb{S}, \quad a\_{ijkl} = \int\_{\mathbb{S}} p\_i p\_j p\_k p\_l \psi(\mathbf{p}) d\mathbb{S} \tag{6}$$

where, S is the surface of a unit sphere and **p** is the unit vector representing an individual fiber. The orientation tensor is symmetric, *aij* = *aji*, and *aijkl* = *aklij* = *ajikl* = *ailkj* = ···. Based on the property that the integral of the probability density function over the entire unit sphere must be 1, the trace of second-order orientation tensor is similarly equal to 1. For example, the second order orientation tensor of randomly distributed fibers has components of *a*11 = *a*22 = *a*33 = 1/3 with all other components being equal to zero.

Using the second-order orientation tensor *aij*, the fourth-order orientation *aijkl* was approximated using the orthotropic fitted (ORT) closure method [42]. Along with the orientation tensors, the unidirectional composite stiffness tensor *Cijkl* was used in the homogenization approach to calculate the anisotropic stiffness of the composite. The underlying unidirectional stiffness tensor of the composite was calculated using the method given by Tandon and Weng [29], which uses the matrix and fiber material properties and the fiber aspect ratio as shown in Table 1. Within the fiber aspect ratio calculations, the weighted average and number average fiber lengths were calculated and the weight average fiber length was used for the aspect ratio based on the discussion as shown in [39].

**Table 1.** Micromechanics input material properties, taken from [39].


From the calculated unidirectional stiffness tensor *Cijkl*, the effective anisotropic stiffness tensor < *Cijkl* > was computed using the equation,

$$\begin{aligned} \left< \mathcal{L}\_{ijkl} \right> &= \left. B\_1 a\_{ijkl} + B\_2 (a\_{ij} \delta\_{kl} + a\_{kl} \delta\_{ij}) + B\_3 (a\_{ik} \delta\_{jl} + a\_{il} \delta\_{jk} + a\_{jl} \delta\_{ik} \\ &+ a\_{jk} \delta\_{il}) + B\_4 (\delta\_{ij} \delta\_{kl}) + B\_5 (\delta\_{ik} \delta\_{jl} + \delta\_{il} \delta\_{jk}) \end{aligned} \right.$$

where *δij* is the Kronecker delta, and the *Bi* terms are shown in terms of the calculated unidirectional stiffness tensor as (see, e.g., [41]),

$$\begin{array}{rclcrcl}B\_{1} &=& \overline{\mathsf{C}}\_{1111} + \overline{\mathsf{C}}\_{2222} - 2\overline{\mathsf{C}}\_{1122} - 4\overline{\mathsf{C}}\_{1212} & & B\_{4} &=& \overline{\mathsf{C}}\_{2233} \\ B\_{2} &=& \overline{\mathsf{C}}\_{1122} - \overline{\mathsf{C}}\_{2233} & & B\_{5} &=& \frac{1}{2}(\overline{\mathsf{C}}\_{2222} - \overline{\mathsf{C}}\_{2233}) \\ B\_{3} &=& \overline{\mathsf{C}}\_{1212} + \frac{1}{2}(\overline{\mathsf{C}}\_{2233} - 2\overline{\mathsf{C}}\_{2222}) \end{array} \tag{8}$$

From the calculated anisotropic stiffness matrix < *Cijkl* >, the fourth-order compliance tensor was calculated using < *Sijkl* > = < *Cijkl* >−1. Using this method, Young's modulus of the isotropic solid shell region can be predicted using the equation *E*11 = 1/*S*1111 and was calculated as 1.73 GPa. The Young's modulus of the foamed core depends on varying densities calculated using Equation (5). As the density of the cells locally increases, the Young's modulus decreases. Based on the comparison study of existing micromechanics models discussed by Zhang et al. [20], the empirical power law equation given by Moore [26] was used in the present study and is given as

$$E\_f = E\_\mathbf{c} \left( 1 - f \right)^n \tag{9}$$

where *Ec* is the Young's Modulus of the solid composite, *Ef* is the Young's modulus of the localized cell region, *f* is the calculated density, and *n* = 2 is given by the square power law relationship. This calculation was performed over each individual grid point mapped onto the surface of the cross-sectioned image and yielded an array of Young's moduli values of *Ef* at 400 × 250 points.

#### *2.6. Finite Element Analysis*

In a previous work [39], the aforementioned micromechanical modeling on coupon samples was experimentally validated for the shell and the core region. The earlier work was extended to include a full non-linear finite element solution, with inputs of the constitutive properties of the neat polymers and fibers. The previous modeling work was extended in the present study to allow for spatial variations of the foam density as indicated in Figure 2 using the image processing technique discussed in the previous sections. The Young's modulus grid developed using Equation (9) was exported into the finite element software package COMSOL (Version 5.4, COMSOL Inc., Burlington, MA, USA) using a custom subroutine written in MATLAB that takes the spatial location as inputs and returns the effective local anisotropic stiffness tensor. The finite element domain is that of a 0.178 m × 0.229 m (7 in × 9 in) tie that is 2.62 m (8.6 ft) in length as shown in Table 2, which is placed in a four-point bend testing configuration conforming to the ASTM D6109-13 dimensions.

**Table 2.** Finite Element Analysis model dimensions.


The domain was given a point cloud of over 10 million points, and the spatially varying local stiffness was computed at each point using the results from the cell homogenization process discussed above. The lookup table data set contains the spatial *x*, *y*, *z* location of the point cloud along with the Young's modulus value. This data set was then exported to a single lookup table for later access by a subroutine within COMSOL using a gridded interpolation technique. This pre-processing work was done to avoid a constant shifting between MATLAB and COMSOL that would significantly increase the computational time. The use of a lookup table resulted in reducing the total computational time of the FEA model to the order of several hours, whereas the prior solution process required more computational time than the authors were willing to consider, and solutions were terminated after a week of computation.

To make use of the available computing power and reduce the computation time, symmetry was used to model the four-point bend test. A quadrilateral mesh was used to model this geometry and is shown in Figure 6. This model was then run with 3.5 million degrees of freedom, and the presented results were compared to results obtained at 8 million degrees of freedom, which was the threshold of the computational resources available on the workstation utilized for the study, an Intel i9-7940X CPU at 3.10 GHz, 14 cores, 28 logical processors, and 64 GB memory. The presented results are graphically indistinguishable from those from the higher resolution model.

**Figure 6.** Quadrilateral Mesh for various surfaces of the symmetric FEA model.

The lookup table was read by the non-linear FEA model to study the load-deflection relationship under a four-point bend test that conforms to ASTM-D6109-13. Load-deflection data samples at regular intervals from 0 to 111 kN (0 to 25 kip) were collected to generate results. To account for the overall deflection, the (*<sup>x</sup>*, *y*, *z*) coordinates of the undeformed structure were mapped to the displaced points as (*<sup>x</sup>*, *y*, *z* − *<sup>w</sup>*), where *w* = *<sup>w</sup>*(*<sup>x</sup>*, *y*, *z*) represents the *z*-deflection in COMSOL. This allowed each point to retain its individual spatially varying stiffness behavior due to cell density variations.

To analyze the importance of including the non-linearity of the composite and spatial variations of the density within the core, the material properties from the uniform core structure were used from the identical material system from the authors' previous work [39] and are given in Table 3 for completeness. This model uses the Ramberg–Osgood [43] model, which is a non-linear curve fitting method as given in COMSOL Multiphysics to analyze the stress *<sup>σ</sup>*(*x*), and strain *<sup>ε</sup>*(*x*), based on Youngs's modulus *E*, a nonlinear reference stress *<sup>σ</sup>ref* , and a nonlinear reference strain *<sup>ε</sup>ref* , and is shown as,

$$\varepsilon(\mathbf{x}) = \frac{\sigma(\mathbf{x})}{E} + \varepsilon\_{ref} (\frac{\sigma(\mathbf{x})}{\sigma\_{ref}})^n \tag{10}$$


**Table 3.** Material inputs for the non-linear FEA model.
