## Abstract

Prismatic closed cells, i.e., honeycomb structures, are often used as infill in additive manufacturing (AM) for providing physical stability to the skin and mechanical integrity to the object. These cells are periodic in nature and uniform in density. In this research, a new fabrication pattern for honeycomb infill is proposed for material deposition-based additive manufacturing applications. The proposed pattern uniformly distributes the material within the cell and can accommodate a controllable variational honeycomb infill while maintaining continuity with relative ease. First, the honeycomb unit cell geometry is defined for uniform and non-uniform voxel sizes. A continuous toolpath scheme is then designed to achieve the honeycomb structure with uniform wall thickness. Unlike traditional honeycomb cells, the aspect ratio of the proposed cell type is not restricted, which helps to introduce variational honeycomb architecture in the infill. Additionally, the proposed cell type is four-time smaller than the traditional cell, which increases the unit cell packing density for the same R3 space. The proposed infill structures are fabricated with both uniform and variational patterns, which are then compared with the traditional honeycomb pattern with compression testing. In comparison to the traditional samples, the proposed uniform and variational infill patterns have achieved higher elastic modulus, collapse strength, and absorbed more specific energy along the X-direction. However, the values measured for both proposed patterns are lower along the Y-direction. Similar results are achieved for two different materials (PLA and TPU), which indicates the consistency of our findings.

## Introduction

In additively manufactured parts, infills are primarily used to support the shell/skin of the part from collapsing and to provide a print base for the skin hanging over the interior of the part. Infills are also used to reduce the material volume by hollowing the interior of the object [1]. Additive manufacturing (AM) offers the capability of making prismatic lattice infills without the use of any special tooling, resulting in lightweight objects with a higher strength to weight ratio. In many instances, the use of infill becomes inevitable because the object skin requires interior support for printing. Furthermore, infills can be tailored to achieve desired part functionality such as guided and localized energy absorption, controlled heat transfer or insulation, part floatation, acoustic wave propagation, or even controlling the center of gravity of the part [2]. A variety of traditional porous infill patterns, such as 0 deg–90 deg zigzag, rectilinear, line, regular hexagonal honeycomb, triangular, cubic, grid, concentric, etc., are available on different existing slicing and 3D printing software [3,4].2

The infill meso-structure, which is in many times designed considering factors such as build time, surface finish, deposition continuity, and support material, has a strong correlation with the physical properties of the object [5]. Both intrinsic and extrinsic properties of the object can be controlled through the infill meso-structure. For instance, the material layout patterns are investigated for structural integrity [6] and other physical attributes including pore size and geometry [7]. Jin et al. [8] proposed a mixed toolpath generation algorithm for contour and zigzag toolpaths to reduce build time while preserving geometrical accuracy. A multipatch toolpath planning concept was proposed where the design domain was discretized into multiple patches, and the unidirectional zigzag toolpath path was created inside each patch to ensure deposition continuity [9]. Adaptive layout patterns [10] have been utilized to achieve desired infill porosity. Recently, a parametric function-based variational deposition pattern [11] has been proposed for tissue scaffolds, which mimics the heterogeneous topology of the native tissue.

Wu et al. [12] proposed a topology optimization approach to generating bone-like porous infill structures. The material distribution across the design domain was optimized under external loading conditions to maximize stiffness, resulting in irregular variational density porous infill. However, at the same time, such complex architectures experience some critical manufacturability issues, such as discontinuous toolpath, overhanging features requiring support material, and trapped voids. Jiang et al. [13] developed a parametric level-set-based topology optimization framework for stereolithography (SLA) AM process to generate dense boundary and periodic internal lattice structure. Cheng et al. [14] proposed an asymptotic homogenization-based topology optimization methodology to design variable density periodic lattice infill structure having improved natural frequency. Heuristic functions and stress or strain fields, i.e., linear interpolation of stress values derived from finite element analysis [15,16] are also considered for generating functional infill patterns. Roger and Krawczak [17] used the stress field to generate discrete varying density infill. However, the discrete interfaces between the high-density and low-density infill areas generated from these methodologies introduce deposition discontinuity and impact their manufacturability. To improve the manufacturability of the topology optimized parts, Mhapsekar et al. [18] integrated two manufacturing constraints, minimization of the number of thin features and volume of support structures, within the density-based topology optimization process. Furthermore, Wang et al. [19] constrained the part boundary slopes in the density-based topology optimization process to control self-support and surface roughness. However, both approaches did not consider the internal architecture of the parts and thus are not applicable to design part lattice infill structures.

Among various infill patterns, 0 deg–90 deg zigzag and regular hexagonal honeycomb with uniform density are the two most commonly used infill patterns. Between these two patterns, the zigzag infill is relatively simple to design with continuity, and printing is straightforward following the vector. A plethora of work has been performed to optimize the zigzag infill for different object geometries by controlling the process parameters [1,20,21]. However, object with zigzag infill demonstrates a higher degree of anisotropy. In contrast, superior structural performance has been reported by the hexagonal honeycomb infill structure [22]. Domínguez-Rodríguez et al. [23] studied the effect of printing orientation and relative density on compressive properties including stiffness of the 3D-printed samples with honeycomb and rectangular infill patterns. They found honeycomb infill patterns to be stiffer and stronger when compared with rectangular infill. However, the honeycomb infill took almost twice the printing time of the rectangular infill. Fernandez-Vicente et al. [24] investigated the mechanical behavior of rectilinear, honeycomb, and line infill patterns for different relative densities. The experimental results revealed that the honeycomb pattern had a better tensile strength [24,25] and impact strength [26] when compared to its counterpart for the same relative density. Using the response surface method, Moradi et al. [27] found that the infill relative density is one of the major parameters influencing the mechanical performance of the structure. With the experimental results, it was claimed that honeycomb was adequate as the infill structure for lightweight parts. Thus, the authors ruled out solid parts as their optimized honeycomb part possessed superior characteristics in terms of mechanical properties, weight, and build time. Additionally, toolpath continuity for honeycomb patterns can ensure better manufacturability. Thus, periodic honeycomb infill patterns have become an attractive option for infill structure in AM.

On the other hand, many researchers are lured by aperiodic or heterogeneous infill due to its potential for near-optimal functionality. Iovenitti et al. [28] 3D-printed polymeric uniform honeycomb structures with different cell wall thicknesses and studied their in-plane compressive behavior and energy-absorbing capability. The printing toolpath was obtained using a commercially available slicing software, and no variational density was incorporated in the structures. Conversely, Bates et al. [29] printed flexible polymer honeycomb structures with graded density in one direction to study energy-absorbing characteristics. The graded density was introduced through the varying wall thickness. Therefore, thicker cell walls requiring multiple passes, which introduced discrete toolpath segments with an increased amount of interruptions and chances of fabrication defects. A Voronoi-based irregular honeycomb-like porous infill proposed by Lu et al. [30] considered part strength–to–weight optimization taking into account the stress working across the part interior. This structure incorporated varying cell size as well as varying cell wall thickness, which can introduce discrete toolpath segments with an increased amount of interruptions resulting in thin features, and fabrication defects. Similarly, the self-supporting rhombic cell lattice infill structures developed by Wu et al. [31] are the macro-scale segmentation of the infill volume. Although minimum printable wall thickness was taken into design consideration, their manufacturability could be a challenge. Rahman et al. [32] investigated the mechanical performance of hierarchical honeycomb metamaterials having printing defects. They found that the presence of printing defects significantly impaired the stiffness of the hierarchical honeycomb structures compared with that of regular honeycomb structures. Cellular metamaterials with an increase in the order of hierarchy become more prone to printing defects due to the exponentially increased number of thin and tiny sections/walls. Therefore, a carefully designed continuous toolpath is an ultimate necessity to alleviate such manufacturing complexity and resulting defects.

Some non-conventional infill patterns, such as skin-frame structures [33] and shape balancing with the center of gravity [2], have been pronounced in literature as well. However, the strut-node frame nature of such structures creates challenges in printing due to tiny layer contour geometries and the resulting increased tool start-stops and air travel. Aremu et al. [34] demonstrated a voxel-based method of generating conformal lattice infills with complete unit cell. In this process, the design domain was first voxelized which is enclosed with net skin. Lattice cells are embedded inside the voxels conform to the arbitrary external geometry via tessellation and trimming. Verma et al. [35] proposed a combinatorial stitching strategy to create variational lattice topology by connecting the non-uniform radii element. Compton et al. [36] presented a variational honeycomb 3D printing pattern that utilizes two-dimensional parameters creating conformal voxelization. However, the manufacturing constraints inherent in some AM techniques, i.e., material extrusion may see deposition discontinuity, poor printing quality, and higher build time. Additionally, such designs are difficult to realize without the use of supports.

Although additive manufacturing offers an immense capability to manufacture any geometry/shape, uniform infill structures have been investigated by many researchers to reduce material consumption, support structure usage, trapped voids, discontinuous toolpath, and build time while achieving enhanced part functionality. However, aperiodic lattice architecture may have the key solution for multi-functionality, but their realization is a challenge due to the design and manufacturing gap. In this paper, a new fabrication pattern for honeycomb infill is proposed for material deposition-based additive manufacturing applications. The proposed pattern can accommodate a controllable variational honeycomb infill without changing the wall thickness while maintaining continuity with relative ease. First, the honeycomb unit cells geometry is defined for uniform and non-uniform voxel size. A continuous toolpath scheme is then designed to achieve the honeycomb structure with uniform wall thickness. The toolpath is characterized and compared with the traditional 3D printing honeycomb toolpath for which the wall thickness is not uniform. The infill structures are fabricated with both uniform and variational patterns, which are then compared with the traditional toolpath pattern with compression testing. The results show that the proposed design demonstrates uniform densification under compression and performs better while absorbing more energy.

## Continuous Toolpath for Honeycomb Lattice Infill

### Definition of Honeycomb Unit Cell.

The traditional uniform hexagonal honeycomb infill structures in AM are constructed with three consecutive layers following a 0 deg–60 deg–120 deg layout pattern. The unit cell construction is given in Fig. 1 demonstrates that each cell wall consists of two adjacent filaments in one layer and one filament in the other two layers creating a non-uniform wall thickness. This results in 2–1–1, 1–2–1, and 1–1–2 formation sequences of filaments making each side of the honeycomb unit cell. However, filaments in two sides of the cell are shared with adjacent cells resulting in 1–.5–.5 sequence of filaments making these two sides as shown in Fig. 1(c).

Fig. 1
Fig. 1
Close modal
To generate gradient as well as uniform density infill and reduce the toolpath overlaps and tool air travel, we are proposing a continuous toolpath scheme for hexagonal honeycomb infill. To model a generic hexagonal honeycomb unit cell, a rectilinear representative unit deemed as a voxel is first defined with parameters Lx and Ly in the standard Cartesian coordinate system, where Lx and Ly are length and width of the voxel along X- and Y-axes, respectively, as shown in Fig. 2(a). It should be mentioned here that both traditional and proposed honeycomb cells are constructed with a stack of three layers which is considered as the height of the voxel fitting the cell. Now, this voxel is discretized along Lx and Ly the honeycomb cell parameters a, b, and c, where c = 2a are derived, as shown in Figs. 2(a) and 2(b). This unit cell is composed of three layers to make continuous toolpath resulting in a 3D unit cell, as shown in Figs. 2(d) and 2(e). The total filament length inside one 3D unit cell will be 4a + 8b + 4(c/2). Here, all the c struts of the unit cell are shared by the adjacent cells, and hence, halves of the c struts belong to the unit cell, which can be seen in Fig. 2(e). Now, if the voxel is discretized into 1/6 units of Lx and Ly, the unit cell parameters can be derived as a = 1/6Ly, $b=1/4Lx2+1/36Ly2$, and c = 2a = 1/3Ly. When the filament diameter is d, the total material volume, Vlattice, inside the 3D unit cell can be calculated using Eq. (1).
$Vlattice=π4d2(4a+8b+2c)=π3d2(Ly+9Lx2+Ly2)$
(1)
Fig. 2
Fig. 2
Close modal
The voxel parameters Lx and Ly may have any of the following relations: Lx = Ly, Lx > Ly, or Lx < Ly. If we consider non-square bounding voxels such that $Lx=1/3Ly$, the resulting honeycomb cell becomes the regular hexagon with all equal sides and 120 deg interior angles between sides, as shown in Fig. 2(b). Thus, the relative density of the cell $ρ*$, which is the ratio of lattice cell material volume Vlattice to solid cell (voxel) material volume, Vsolid, can be represented in terms of Ly (Eq. (2)).
$ρ*=VlatticeVsolid=πd3Ly$
(2)
Therefore, for a given infill density $ρ*$ of a layer, the voxel parameters for the regular hexagonal honeycomb cell can be determined using Eqs. (3) and (4). Then, the infill layer can be voxelized, and honeycomb cells can be fitted to generate the layer-by-layer continuous tool-paths as represented in the “Infill Space Voxelization” section:
$Lx=πd3ρ*$
(3)
$Ly=πd3ρ*$
(4)

In the proposed honeycomb cells, two filament segments appear along each unit cell wall in the three-layer arrangement (see Fig. 2). The unit cell construction shown in Fig. 2(e) demonstrates 1–0–1, 1–1–0, and 0–1–1 formation sequence of filaments making the sides of the honeycomb unit cell in the three-layer arrangement.

In contrast, the traditional 3D-printed honeycomb cells are regular hexagon with all equal sides and 120 deg interior angles between sides. Thus, analogous to our proposed honeycomb cell definition, the traditional honeycomb infill can also be represented with voxels parameterized by $Lx=1/3Ly$, a = 1/6Ly, and b = c = 1/3Ly as shown in Fig. 3. The 3D unit cell can again be easily decomposed into three layers as shown in Figs. 3(d) and 3(e). The total filament length inside one 3D unit cell will be 8a + 16b + 4(c/2) + 2c. Here, the c struts in the top two layers (blue and green colored layers in Figs. 3(d) and 3(e)) of the unit cell are shared by the adjacent cells. Thus, half of each c strut in the top two layers and a whole of each c strut in the bottom layer belong to the unit cell. Considering the filament diameter as d, the total material volume, Vlattice, inside the 3D unit cell can be calculated using Eq. (5).
$Vlattice=π4d2(8a+16b+4c)=2πd2Ly$
(5)
Fig. 3
Fig. 3
Close modal
Thus, the relative density of the cell $ρ*$ can be expressed in terms of Ly by Eq. (6).
$ρ*=VlatticeVsolid=2πd3Ly$
(6)
Therefore, for a given infill density $ρ*$ of a layer, the voxel parameters for the traditional regular hexagonal honeycomb cell can be determined using Eqs. (7) and (8).
$Ly=2πd3ρ*$
(7)
$Lx=2πd3ρ*$
(8)

Therefore, Eqs. (2) and (6) indicate that the relative density of the traditional honeycomb unit cell is twice the relative density of the proposed 1/6 division honeycomb unit cell for the same size of a voxel. Consequently, for a given relative density $ρ*$, the proposed cell size parameters (Lx and Ly) will be half the size of traditional cell parameters. This unit cell size and density relationship for both proposed and traditional honeycomb cells are further illustrated with Fig. 4. Thus, the proposed honeycomb infill design accommodates smaller cells (half of its counterpart) compared to the existing regular hexagonal honeycomb infills for the same relative density. This feature of the infill can provide better support to the object skin resulting in enhanced surface quality and wall strength.

Fig. 4
Fig. 4
Close modal
For a square voxel (Lx = Ly = L) discretized into 1/8 units of Lx and Ly, if the cell parameters are considered as a = 1/8Ly, $b=1/44Lx2+Ly2$, and c = 2a = 1/4Ly (see Fig. 5), the total material volume (Vlattice) inside the 3D unit cell can be obtained from Eq. (9).
$Vlattice=π4(1+25)d2L$
(9)
Fig. 5
Fig. 5
Close modal
In a similar fashion, the relative density of the cell $ρ*$ can be represented by Eq. (10). Consequently, the voxel parameters can be determined using Eq. (11).
$ρ*=(1+25)πd12L$
(10)
$L=(1+25)πd12ρ*$
(11)

From Fig. 6, the proposed and traditional cell patterns can be compared based on how many cells can be fitted in the same given space. Compared to the traditional honeycomb cell, the proposed cell types (both 1/6 and 1/8 division types) occupy one-fourth of the space for the same relative density. For this same reason, a substantially higher number of unit cells can be packed in the same space, which provides the opportunity to voxelized the design space with higher resolution, as shown in Fig 6. The density versus cell number plot in Fig. 6 shows that the measured density reaches more than 100% with 16 traditional type cell packing whereas our proposed cell types require 64 cells packing. Although this is numerically possible and can be demonstrated in wireframe models given in Fig. 6, achieving density over 100% during manufacturing is practically infeasible since the volume of material deposited inside a cell will be greater than the cell volume.

Fig. 6
Fig. 6
Close modal

### Infill Space Voxelization.

For a given relative density of each infill layer/slice, the voxel parameters (size) Lx and Ly can be derived from either Eqs. (3) and (4) or Eqs. (10) and (11). Then, the infill layer can be voxelized using the calculated voxel parameters as shown in Fig. 7. Voxelization is performed with parallel grid lines of spacing Lx and Ly along X and Y directions, respectively, in the standard coordinate system. The start points along both directions are the lower extreme values of X and Y. The intersection points of the grid lines determine the voxels. After voxelization, honeycomb parameters a, b, and c are determined and honeycomb cells are fitted inside the voxels as represented in Fig. 7. For honeycomb structures with variational density distribution, the relation between Lx and Ly will vary throughout the structure in accordance with the density distribution. Therefore, elongated/flattened hexagonal cells and resulting in irregular interstitial hexagonal cells will appear as depicted in Fig. 7.

Fig. 7
Fig. 7
Close modal

For variational infill, the relative density may follow continuous distributions $ρx*$ and $ρy*$ along the X and Y-directions, respectively. In this case, the parallel grid lines along both X and Y-directions are determined by plugging the locations of the immediate previous X and Y grid lines, respectively, into the continuous density functions. Hence, like other unit cell-based cellular structures, this voxelization represents a discrete approximation of the continuous density distributions. Furthermore, the resultant density along the structure is determined with Lx and Ly grid lines that are controlled independently. Therefore, although the voxelization may influence the density of the neighboring cells/voxels due to the rectilinear nature of the grid lines, this process can minimize overall density approximation error.

### Continuous Honeycomb Toolpath Scheme.

Our proposed honeycomb cells can be decomposed into three distinct layers of parallel kinked lines, as shown in Fig. 8. The parallel lines in each of the three layers are kinked in a certain pattern following the honeycomb cell geometry. Each layer's parallel kinked lines can be connected in a continuous zigzag pattern. Stacking these three sublayers results in complete honeycomb cells, as shown in Fig. 8. Repetitively putting this tri-layer one upon another will generate 3D honeycomb lattice. Thus, each layer in this lattice renders continuous toolpath for 3D printing.

Fig. 8
Fig. 8
Close modal

## Sample Parts Fabrication

The proposed methodology was implemented with a visual basic-based scripting language and requires only a few minutes to process and generate the patterns. In order to perform mechanical tests, variational and uniform toolpaths of a series of sample parts with dimensions of 40 × 40 × 12.6 mm3 and relative densities of 0.23 were designed and fabricated. The generated continuous honeycomb toolpaths were printed with an extrusion-based AM process (Creality Ender 3 Pro 3D Printer) using PLA as well as TPU materials. Traditional honeycomb toolpaths generated by the open-source slic3r software3 were fabricated. All the samples were modeled with the same dimensions and the same relative density for the compression test. All fabricated test samples were printed with no skin; therefore, they all had rough sides where the honeycomb was cut off. A sine function $(ρx,y*(φ)=1−sinφ,0.25π≤$$φ≤0.75π)$ was used to generate the density gradient along with the sample structures, as shown in Fig. 9. Equation (7) was used to determine the parallel grid lines for voxelization. Figure 10 depicts the proposed and traditional honeycomb toolpaths. For all the fabricated samples, the layer thickness used was 0.2 mm. Pictures of the fabricated samples are shown in Fig. 11.

Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal

## Mechanical Tests

Honeycomb structures are used for load-bearing/absorbing applications since they provide superior performance under compression loading [17]. Numerous research has reported the compressive/crushing behavior and energy absorption capability of different type of 3D-printed lattice structures, for instance, regular honeycombs [28,29,37,38], Ti-6Al-4V regular diamond lattice [39], BCC and BCC-Z lattices [40,41], BCC and BCC-Z, and micro-lattices [42,43], gyroid lattice [44], etc. To characterize the mechanical behavior of the proposed honeycomb pattern, we investigate their stress–strain behavior, and the energy-absorbing capabilities under in-plane compression loading for a rigid thermoplastic (PLA) material as well as a flexible thermoplastic (TPU) material.

The in-plane compression tests were performed on the fabricated honeycomb structures with an Instron 5567 universal testing machine utilizing 30 kN load cell. The crosshead speeds applied for PLA and TPU samples were 4 mm/min and 13 mm/min, respectively. The tests were controlled, and the data were recorded using bluehill software for Instron. All the tests were recorded with a high-speed camera and pictures at different strain levels were snapped from the video to illustrate the deformation and progressive failure mechanism of the lattice structures. In all cases, the Z-direction indicates the build direction for the designed and fabricated structures. Therefore, in-plane compression loads were applied to the proposed honeycomb structures along the X and Y-directions to capture the mechanical anisotropy. The traditional uniform honeycombs generated by Slic3r are expected to have in-plane isotropic properties [23]. Although in our experiment, the orientation of the traditional cells is inclined compared with that in the proposed cell, we anticipate that this is not going to influence the compression result of the traditional ones. Hence, they were tested along one in-plane direction. For each type of structure, i.e., uniform, variation, and traditional, three samples were tested along each of the two in-plane directions (X and Y). Also, all the sample structures were designed with the same overall relative density.

The nominal stress (σlattice) of the lattice structures is measured by dividing the applied load by the initial cross-sectional area of the structures perpendicular to the in-plane load. The nominal lattice strains (ɛlattice) were measured from the deflection of the interface between the structures and the compression platens. The stress–strain curves obtained after averaging the measured data from three replicates of each lattice type were used to study the compression behavior and energy-absorbing capability of the sample structures.

## Results and Discussion

Pictures in Fig. 12 demonstrate the failure process of the proposed toolpath pattern (variational and uniform) as well as traditional PLA honeycomb samples, which are subject to compression loading along the X-direction, at the overall lattice strains (ɛlattice) of 0%, 20%, 40%, and 60%. The failure of regular hexagonal honeycomb cells subject to a similar in-plane X-directional loading typically propagates through deformation localization along two opposite oblique (diagonal) bands intersecting each other [28]. At $εlattice=20%$, it can also be observed for the proposed variational honeycomb that its deformation occurred through the failure of the diagonal cells forming two opposite bands resembling an “X” shape. Additionally, the deformation further localized through the plastic collapse of larger cells in the mid-region where the two opposite bands intersected. As we increase the load $(εlattice=40%)$, cells in the mid-region start collapsing along the Y-direction. The lower density mid-region with larger cells was basically responsible for this failure. The deformation of the proposed uniform honeycomb (middle row in Fig. 12) started evenly through the collapse of mostly diagonal cells forming almost perfect “X” shaped bands like the regular hexagonal honeycomb cells. However, at $εlattice=40%$, the failure is more prominent in the bottom half of the structure. This can be attributed to a relatively rough bottom edge of that structure that occurred during its fabrication. The traditional uniform honeycomb started collapsing asymmetrically and randomly from the bottom side and then progressively collapsed through the neighboring cells with the shearing mode of deformation.

Fig. 12
Fig. 12
Close modal

Figure 13 shows the failure processes of the PLA honeycomb samples identical to those discussed above, which are now subject to compression loading along the Y-direction. In case of the proposed variational honeycomb structure, cells in the mid-region were primarily deformed along the X-direction which dominated the entire failure of the proposed variational honeycomb. The larger cells in the mid-region guided the failure process. For the proposed uniform honeycomb, although yielding started through diagonal cells, the structure progressively failed through an asymmetric shearing mode of deformation in horizontal rows of cells along the mid-region. In both proposed structures, deformations of the cells away from the mid-region were quite symmetric and homogeneous, such as regular hexagonal honeycombs [45,46].

Fig. 13
Fig. 13
Close modal

Despite being regular hexagon, the failure of traditional honeycomb samples did not predominantly occur through either failure of diagonal cells or localized deformation along mid-region in both X and Y-direction loading conditions, respectively. Instead, the deformation began with buckling of cell walls and then localized with the asymmetric shearing mode in more than one place throughout the structure. Furthermore, uniform densification can be observed in all the proposed toolpath samples compared to the traditional toolpath model. This can be attributed to the non-uniform wall thickness of the traditional honeycomb at their design phase. The thicker part of the wall segment acts as reinforcement causing random densification.

The stress (σlattice) − strain (ɛlattice) behavior of the PLA honeycombs can be observed in Fig. 14(a). The linear elasticity of all the samples spans over a very small strain range. The lattice elastic moduli Elattice are determined from this strain range and are listed in Table 1. The linear elasticity ends through plastic failure at the plastic strain $εlatticepl$ giving the plastic collapse strength $σlatticepl$. After plastic failure starts, the honeycombs show long plastic plateau regions until densification starts at the densification strain $εlatticeD$. Table 1 lists the values of the plastic collapse strengths, plastic strains, and elastic moduli of the PLA samples. It can be observed that the proposed variational structures have the largest modulus and collapse strength along the X-direction and start collapsing at the lowest strain. Also, both variational and uniform structures proposed here have similar collapse strength along the Y-direction. Overall, both proposed design structures (variational and uniform) demonstrate superiority in terms of strength and modulus compared to the traditional pattern structure.

Fig. 14
Fig. 14
Close modal
Table 1

Plateau stress, plateau strain, and elastic modulus PLA honeycomb lattices

Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$σlatticepl$ (MPa)2.1821.1811.8421.1511.588
$εlatticepl$ (%)4.5887.1684.72211.1104.975
Elattice (MPa)67.323.858.719.443.2
Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$σlatticepl$ (MPa)2.1821.1811.8421.1511.588
$εlatticepl$ (%)4.5887.1684.72211.1104.975
Elattice (MPa)67.323.858.719.443.2

It is also noticeable in Fig. 14 and Table 1 that the mechanical properties of the proposed variational and uniform structures along the X-direction differ from those along the Y-direction. Although an identical density distribution is followed along both directions for these structures, this mechanical property variation can be attributed to the elongation of honeycomb cells along the X-direction in the design (see Figs. 9(d)10(a) and 10(b)). More specifically, the elongated cells contributed to an increase in both strength and modulus along the X-direction and a reduction in plastic strain. This in-plane anisotropy in the proposed structures also results in dissimilar failure processes for X and Y-direction loadings as observed in Figs. 12 and 13, respectively.

The compression behavior along X and Y-directions of the proposed variational and traditional honeycombs made with TPU material were also observed. Figure 15 is showing only the compression behavior of variational and traditional TPU structures under the X-direction loading. Again, the failure of both variational and traditional structures predominantly occurred in the mid-region of the samples. The stress (σlattice) − strain (ɛlattice) curves of the TPU honeycombs given in Fig. 14(b) demonstrate that they also have low-strain elastic regions followed by long plastic plateau up to the densification strain. Uniform densification is also observed in the proposed honeycomb patterns.

Fig. 15
Fig. 15
Close modal

Table 2 lists the values of the plastic collapse strengths, plastic strains, and elastic moduli of the TPU samples. Compared to traditional honeycomb samples, the proposed variational as well as uniform TPU honeycombs have slightly higher plastic strength and elastic modulus and start collapsing at an almost similar strain level under the X-direction loading. The plastic collapse strengths of variational and uniform structures are also comparable with that of traditional one along the Y-direction. Thus, for TPU materials, both proposed design structures (variational and uniform) also demonstrate superiority in terms of strength and modulus compared with the traditional pattern structure.

Table 2

Plateau stress, plateau strain, and elastic modulus of TPU honeycomb lattices

Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$σlatticepl$ (MPa)0.2470.1290.2160.1520.157
$εlatticepl$ (%)6.1338.6826.96113.8807.062
Elattice (MPa)5.92.04.91.33.6
Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$σlatticepl$ (MPa)0.2470.1290.2160.1520.157
$εlatticepl$ (%)6.1338.6826.96113.8807.062
Elattice (MPa)5.92.04.91.33.6

The energy-absorbing capabilities of the honeycomb structures were determined through the numerical integration of the area under the stress (σlattice) − strain (ɛlattice) curves. This can also be called specific energy absorbed or the cumulative energies absorbed per unit volume (EV). Figure 16 demonstrates the specific energies absorbed by PLA and TPU honeycombs. The strains $(εlatticeD)$ where densification starts and the total energies absorbed per unit volume $(EVD)$ of the PLA and TPU structures were determined and are listed in Tables 3 and 4, respectively. The linear portions of the curves starting from the plastic strain $(εlatticepl)$ moving to densification strain $(εlatticeD)$ where the specific energy is proportional to the lattice strain refers to the long plastic plateau of the honeycombs. The densification of the proposed variational and uniform honeycombs, and the traditional honeycombs started in the range of 55–65% strain.

Fig. 16
Fig. 16
Close modal
Table 3

Densification strains and total specific energies absorbed at densification of PLA honeycombs

Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$εlatticeD$ (%)55.0050.6860.3152.9451.1
$EVD$ (J/m3)0.902 × 1060.578 × 1061.135 × 1060.611 × 1060.671 × 106
Proposed variationalProposed uniformTraditional uniform
X-directionY-directionX-directionY-direction
$εlatticeD$ (%)55.0050.6860.3152.9451.1
$EVD$ (J/m3)0.902 × 1060.578 × 1061.135 × 1060.611 × 1060.671 × 106
Table 4

Densification strains and energies absorbed per unit volume at densification of TPU honeycombs

Proposed variationalProposed uniformTraditional uniform
X directionY directionX directionY direction
$εlatticeD$ (%)59.3254.1557.4651.9857.16
$EVD$ (J/m3)1.336 × 1050.876 × 1051.343 × 1050.884 × 1051.126 × 105
Proposed variationalProposed uniformTraditional uniform
X directionY directionX directionY direction
$εlatticeD$ (%)59.3254.1557.4651.9857.16
$EVD$ (J/m3)1.336 × 1050.876 × 1051.343 × 1050.884 × 1051.126 × 105

For PLA material, both variational and uniform patterns proposed in this paper absorbed 34% and 69% more specific energy $(EVD)$ along the X-direction than the traditional samples. On the other hand, densifications of the proposed variational and uniform honeycombs along the Y-direction started nearly at the same strain level, but they have relatively lower (less than 14%) total specific energy absorbed $(EVD)$ compared with the traditional honeycomb samples.

Similarly, for TPU material, both variational and uniform patterns proposed in this paper absorbed about 19% more specific energy $(EVD)$ along the X-direction than the traditional samples. However, densifications of the proposed patterns along the Y-direction started earlier, and they have around 21% lower total specific energy absorbed $(EVD)$ compared to the traditional samples.

The bison example is the demonstration of our proposed method being applicable to the freeform 3D shape.

For the demonstration of our proposed method being applicable to free-form 3D shape, we used a Bison model as an example shown in Fig. 17. The boundary contours were first obtained through slicing the model with a given slice thickness. Then, following the method discussed in previous sections, the layers were voxelized and variational infill was generated. For the demonstration purpose, we have shown a tri-layer continuous toolpath generation process for a Bison layer in Fig. 17.

Fig. 17
Fig. 17
Close modal

## Conclusion

Continuity-based toolpaths for both variational and uniform honeycomb infill structures of 3D-printed objects are proposed and characterized in this paper. The proposed and traditional honeycomb infill fabrication patterns were fabricated using a commercially available 3D printer. Mechanical tests were conducted to investigate the in-plane compressive deformation behavior. The study of the deformation process reveals that the proposed honeycomb patterns possess relatively higher strength, stiffness, and total specific energy absorbed at densification along the one (X) of the two orthogonal directions on the build (XY) plane compared to the traditional equivalent relative density honeycomb infills. Unlike traditional honeycomb cells, the aspect ratio of the proposed cell type is not restricted, which can produce longer or wider cell type fitting irregular geometry. This helps to introduce variational honeycomb architecture in the infill but may bring anisotropic mechanical behavior.

Furthermore, the proposed honeycomb infill design accommodates smaller cells compared to the traditional regular hexagonal honeycomb infills for the same relative density. This feature of the infill can provide higher cell-packing and better support to the object skin, which may result in enhanced surface quality. That may have helped avoid sagging in our samples, but for larger unit cell, sagging may become unavoidable. The variational honeycomb pattern introduced a unique failure mechanism with uniform densification, which could be studied further. Thus, the controlled variational nature of the honeycomb can be utilized for very specific applications such as guided and localized energy absorption, heat transfer, acoustic wave propagation, etc.

3

See Note 2.

## Acknowledgment

The authors thank Stephen Abbadessa and UG student Brandon Johnstone from UMaine for their help in printing and testing some of the specimens. Financial supports provided by the National Science Foundation Grant # OIA-1355466 and by the US-DOT # 693JK31850009CAAP are acknowledged.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. No data, models, or codes were generated or used for this paper.

## References

1.
Ahsan
,
A. M. M. N.
,
Habib
,
M. A.
, and
Khoda
,
B.
,
2015
, “
Resource Based Process Planning for Additive Manufacturing
,”
Comput.-Aided Des.
,
69
, pp.
112
125
. 10.1016/j.cad.2015.03.006
2.
Prevost
,
R.
,
Whiting
,
E.
,
Lefebvre
,
S.
, and
Sorkine-Hornung
,
O.
,
2013
, “
Make it Stand: Balancing Shapes for 3D Fabrication
,”
ACM Trans. Graph.
,
32
(
4
), pp.
1
10
. 10.1145/2461912.2461957
3.
Ultimaker
. CURA version 4.4. [cited 2020 January]. Available from https://ultimaker.com/software/ultimaker-cura
4.
Stratasys
. Insight version 5.1. [cited 2015 December]. Available from https://support.stratasys.com/en/resources/software-download
5.
Hansen
,
C. J.
,
Saksena
,
R.
,
Kolesky
,
D. B.
,
Vericella
,
J. J.
,
Kranz
,
S. J.
,
Muldowney
,
G. P.
,
Christensen
,
K. T.
, and
Lewis
,
J. A.
,
2013
, “
Inkjet Printing: High-Throughput Printing via Microvascular Multinozzle Arrays (Adv. Mater. 1/2013)
,”
Adv. Mater.
,
25
(
1
), p.
2
.
6.
Es-Said
,
O. S.
,
Foyos
,
J.
,
Noorani
,
R.
,
Mendelson
,
M.
,
Marloth
,
R.
, and
Pregger
,
B. A.
,
2000
, “
Effect of Layer Orientation on Mechanical Properties of Rapid Prototyped Samples
,”
Mater. Manuf. Processes
,
15
(
1
), pp.
107
122
. 10.1080/10426910008912976
7.
Domingos
,
M.
,
Chiellini
,
F.
,
Cometa
,
S.
,
Giglio
,
E. D.
,
Grillo-Fernandes
,
E.
,
Bartolo
,
P.
, and
Chiellini
,
E.
,
2011
, “
Evaluation of in Vitro Degradation of PCL Scaffolds Fabricated via BioExtrusion—Part 2: Influence of Pore Size and Geometry
,”
Virtual Phys. Prototyping
,
6
(
3
), pp.
157
165
. 10.1080/17452759.2011.605839
8.
Jin
,
G. Q.
,
Li
,
W. D.
, and
Gao
,
L.
,
2013
, “
An Adaptive Process Planning Approach of Rapid Prototyping and Manufacturing
,”
Robot. Comput. Integr. Manuf.
,
29
(
1
), pp.
23
38
. 10.1016/j.rcim.2012.07.001
9.
Yu
,
H.
,
Hong
,
H.
,
Cao
,
S.
, and
Ahmad
,
R.
,
2020
, “
Topology Optimization for Multipatch Fused Deposition Modeling 3D Printing
,”
Appl. Sci.
,
10
(
3
), p.
943
. 10.3390/app10030943
10.
Khoda
,
A. K. M. B.
, and
Koc
,
B.
,
2013
, “
Functionally Heterogeneous Porous Scaffold Design for Tissue Engineering
,”
Comput.-Aided Des.
,
45
(
11
), pp.
1276
1293
. 10.1016/j.cad.2013.05.005
11.
Ahsan
,
A. M. M.
,
Xie
,
R. N.
, and
Khoda
,
B.
,
2018
, “
Heterogeneous Topology Design and Voxel-Based bio-Printing
,”
Rapid Prototyp. J.
,
24
(
7
), pp.
1142
1154
. 10.1108/RPJ-05-2017-0076
12.
Wu
,
J.
,
Aage
,
N.
,
Westermann
,
R.
, and
Sigmund
,
O.
,
2018
, “
Infill Optimization for Additive Manufacturing—Approaching Bone-Like Porous Structures
,”
IEEE Trans. Vis. Comput. Graph.
,
24
(
2
), pp.
1127
1140
. 10.1109/TVCG.2017.2655523
13.
Jiang
,
L.
,
Ye
,
H.
,
Zhou
,
C.
, and
Chen
,
S.
,
2019
, “
Parametric Topology Optimization Toward Rational Design and Efficient Prefabrication for Additive Manufacturing
,”
ASME J. Manuf. Sci. Eng.
,
141
(
4
), p.
041007
. 10.1115/1.4042580
14.
Cheng
,
L.
,
Liang
,
X.
,
Belski
,
E.
,
Wang
,
X.
,
Sietins
,
J. M.
,
Ludwick
,
S.
, and
To
,
A.
,
2018
, “
Natural Frequency Optimization of Variable-Density Additive Manufactured Lattice Structure: Theory and Experimental Validation
,”
ASME J. Manuf. Sci. Eng.
,
140
(
10
), p.
105002
. 10.1115/1.4040622
15.
Steuben
,
J. C.
,
Iliopoulos
,
A. P.
, and
Michopoulos
,
J. G.
,
2016
, “
Implicit Slicing for Functionally Tailored Additive Manufacturing
,”
Comput. Aided Des.
,
77
, pp.
107
119
. 10.1016/j.cad.2016.04.003
16.
Adams
,
D.
, and
Turner
,
C. J.
,
2018
, “
An Implicit Slicing Method for Additive Manufacturing Processes
,”
Virtual Phys. Prototyping
,
13
(
1
), pp.
2
7
. 10.1080/17452759.2017.1392684
17.
Roger
,
F.
, and
Krawczak
,
P.
,
2015
, “
3D-printing of Thermoplastic Structures by FDM Using Heterogeneous Infill and Multi-Materials: An Integrated Design-Advanced Manufacturing Approach for Factories of the Future
,”
22ème Congrès Français de Mécanique 2015 (CFM 2015)
,
Lyon, France
.
18.
Mhapsekar
,
K.
,
McConaha
,
M.
, and
Anand
,
S.
,
2018
, “
Additive Manufacturing Constraints in Topology Optimization for Improved Manufacturability
,”
ASME J. Manuf. Sci. Eng.
,
140
(
5
), p.
051017
. 10.1115/1.4039198
19.
Wang
,
C.
,
Qian
,
X.
,
Gerstler
,
W. D.
, and
Shubrooks
,
J.
,
2019
, “
Boundary Slope Control in Topology Optimization for Additive Manufacturing: for Self-Support and Surface Roughness
,”
ASME J. Manuf. Sci. Eng.
,
141
(
9
), p.
091001
. 10.1115/1.4043978
20.
Habib
,
A.
,
Ahsan
,
N.
, and
Khoda
,
B.
,
2015
, “
Optimizing Material Deposition Direction for Functional Internal Architecture in Additive Manufacturing Processes
,”
Procedia Manuf.
,
1
, pp.
378
392
. 10.1016/j.promfg.2015.09.045
21.
Ahsan
,
N.
, and
Khoda
,
B.
,
2016
, “
AM Optimization Framework for Part and Process Attributes Through Geometric Analysis
,”
Addit. Manuf.
,
11
, pp.
85
96
.
22.
Ćwikła
,
G.
,
Grabowik
,
C.
,
Kalinowski
,
K.
,
Paprocka
,
I.
, and
Ociepka
,
P.
,
2017
, “
The Influence of Printing Parameters on Selected Mechanical Properties of FDM/FFF 3D-Printed Parts
,”
Mater. Sci. Eng.
,
227
, p.
012033
. 10.1088/1757-899X/227/1/012033
23.
Domínguez-Rodríguez
,
G.
,
Ku-Herrera
,
J. J.
, and
Hernández-Pérez
,
A.
,
2018
, “
An Assessment of the Effect of Printing Orientation, Density, and Filler Pattern on the Compressive Performance of 3D Printed ABS Structures by Fuse Deposition
,”
Int. J. Adv. Manuf. Technol.
,
95
(
5
), pp.
1685
1695
. 10.1007/s00170-017-1314-x
24.
Fernandez-Vicente
,
M.
,
Calle
,
W.
,
Ferrandiz
,
S.
, and
Conejero
,
A.
,
2016
, “
Effect of Infill Parameters on Tensile Mechanical Behavior in Desktop 3D Printing
,”
3D Print. Addit. Manuf.
,
3
(
3
), pp.
183
192
. 10.1089/3dp.2015.0036
25.
Bartolai
,
J.
,
Wilson-Heid
,
A. E.
,
Kruse
,
J. R.
,
Beese
,
A. M.
, and
Simpson
,
T. W.
,
2019
, “
Full Field Strain Measurement of Material Extrusion Additive Manufacturing Parts With Solid and Sparse Infill Geometries
,”
JOM
,
71
(
3
), pp.
871
879
. 10.1007/s11837-018-3217-1
26.
Ramkumar
,
P. L.
,
2019
, “
Investigation on the Effect of Process Parameters on Impact Strength of Fused Deposition Modelling Specimens
,”
IOP Conference Series: Materials Science and Engineering
,
491
, p.
012026
.
27.
Moradi
,
M.
,
Meiabadi
,
S.
, and
Kaplan
,
A.
,
2019
, “
3D Printed Parts With Honeycomb Internal Pattern by Fused Deposition Modelling; Experimental Characterization and Production Optimization
,”
Metals Mater. Int.
,
25
(
5
), pp.
1312
1325
. https://doi.org/10.1007/s12540-019-00272-9
28.
Iovenitti
,
P.
,
Masood
,
S. H.
, and
Nikzad
,
M.
,
2017
, “
In-Plane Energy Absorption Evaluation of 3D Printed Polymeric Honeycombs
,”
Virtual Phys. Prototyping
,
12
(
2
), pp.
117
131
. 10.1080/17452759.2017.1291354
29.
Bates
,
S. R. G.
,
Farrow
,
I. R.
, and
Trask
,
R. S.
,
2019
, “
Compressive Behaviour of 3D Printed Thermoplastic Polyurethane Honeycombs With Graded Densities
,”
Mater. Des.
,
162
, pp.
130
142
. 10.1016/j.matdes.2018.11.019
30.
Lu
,
L.
,
Sharf
,
A.
,
Zhao
,
H.
,
Wei
,
Y.
,
Fan
,
Q.
,
Chen
,
X.
,
Savoye
,
Y.
,
Tu
,
C.
,
Cohen-Or
,
D.
, and
Chen
,
B.
,
2014
, “
Build-to-Last: Strength to Weight 3D Printed Objects
,”
ACM Trans. Graph.
,
33
(
4
), pp.
1
10
.
31.
Wu
,
J.
,
Wang
,
C. C. L.
,
Zhang
,
X.
, and
Westermann
,
R.
,
2016
, “
Self-Supporting Rhombic Infill Structures for Additive Manufacturing
,”
Comput.-Aided Des.
,
80
, pp.
32
42
. 10.1016/j.cad.2016.07.006
32.
Rahman
,
K. M.
,
Hu
,
Z.
, and
Letcher
,
T.
,
2018
, “
In-Plane Stiffness of Additively Manufactured Hierarchical Honeycomb Metamaterials With Defects
,”
ASME J. Manuf. Sci. Eng.
,
140
(
1
), p.
011007
. 10.1115/1.4038205
33.
Wang
,
W.
,
Wang
,
T. Y.
,
Yang
,
Z.
,
Liu
,
L.
,
Tong
,
X.
,
Tong
,
W.
,
Deng
,
J.
,
Chen
,
F.
, and
Liu
,
X.
,
2013
, “
Cost-Effective Printing of 3D Objects With Skin-Frame Structures
,”
ACM Trans. Graph.
,
32
(
6
), pp.
1
10
.
34.
Aremu
,
A. O.
,
Brennan-Craddock
,
J. P. J.
,
Panesar
,
A.
,
Ashcroft
,
I. A.
,
Hague
,
R. J. M.
,
Wildman
,
R. D.
, and
Tuck
,
C.
,
2017
, “
A Voxel-Based Method of Constructing and Skinning Conformal and Functionally Graded Lattice Structures Suitable for Additive Manufacturing
,”
Addit. Manuf.
,
13
, pp.
1
13
.
35.
Verma
,
C. S.
,
Rankouhi
,
B.
, and
Suresh
,
K.
,
2020
, “
A Combinatorial Approach for Constructing Lattice Structures
,”
ASME J. Mech. Des.
,
142
(
4
), p.
041404
. 10.1115/1.4044521
36.
Compton
,
B. G.
,
Goin
,
M. J.
,
Shams
,
H.
, and
Shylo
,
O.
,
2019
, Methods of Producing a Cellular Structure and Articles Produced Thereby, University of Tennessee Research Foundation.
37.
Bates
,
S. R. G.
,
Farrow
,
I. R.
, and
Trask
,
R. S.
,
2016
, “
3D Printed Polyurethane Honeycombs for Repeated Tailored Energy Absorption
,”
Mater. Des.
,
112
, pp.
172
183
. 10.1016/j.matdes.2016.08.062
38.
Pollard
,
D.
,
Ward
,
C.
,
Herrmann
,
G.
, and
Etches
,
J.
,
2017
, “
The Manufacture of Honeycomb Cores Using Fused Deposition Modeling
,”
Adv. Manuf.: Polym. Compos. Sci.
,
3
(
1
), pp.
21
31
. 10.1080/20550340.2017.1306337
39.
Van Grunsven
,
W.
,
Hernandez-Nava
,
E.
,
Reilly
,
G. C.
, and
Goodall
,
R.
,
2014
, “
Fabrication and Mechanical Characterisation of Titanium Lattices With Graded Porosity
,”
Metals
,
4
(
3
), pp.
401
409
. 10.3390/met4030401
40.
Maskery
,
I.
,
Hussey
,
A.
,
Panesar
,
A.
,
Aremu
,
A.
,
Tuck
,
C.
,
Ashcroft
,
I.
, and
Hague
,
R.
,
2017
, “
An Investigation Into Reinforced and Functionally Graded Lattice Structures
,”
J. Cell. Plast.
,
53
(
2
), pp.
151
165
. 10.1177/0021955X16639035
41.
Azzouz
,
L.
,
Chen
,
Y.
,
Zarrelli
,
M.
,
Pearce
,
J. M.
,
Mitchell
,
L.
,
Ren
,
G.
, and
Grasso
,
M.
,
2019
, “
Mechanical Properties of 3-D Printed Truss-Like Lattice Biopolymer non-Stochastic Structures for Sandwich Panels With Natural Fibre Composite Skins
,”
Compos. Struc.
,
213
, pp.
220
230
. 10.1016/j.compstruct.2019.01.103
42.
Gumruk
,
R.
, and
Mines
,
R. A. W.
,
2013
, “
Compressive Behaviour of Stainless Steel Micro-Lattice Structures
,”
Int. J. Mech. Sci.
,
68
, pp.
125
139
. 10.1016/j.ijmecsci.2013.01.006
43.
Gümrük
,
R.
,
Mines
,
R. A. W.
, and
Karadeniz
,
S.
,
2013
, “
Static Mechanical Behaviours of Stainless Steel Micro-Lattice Structures Under Different Loading Conditions
,”
Mater. Sci. Eng., A
,
586
, pp.
392
406
. 10.1016/j.msea.2013.07.070
44.
Yan
,
C. Z.
,
Hao
,
L.
,
Hussein
,
A.
, and
Raymont
,
D.
,
2012
, “
Evaluations of Cellular Lattice Structures Manufactured Using Selective Laser Melting
,”
Int. J. Mach. Tools Manuf.
,
62
, pp.
32
38
. 10.1016/j.ijmachtools.2012.06.002
45.
Gibson
,
L. J.
, and
Ashby
,
M. F.
,
1999
,
Cellular Solids: Structure and Properties
,
Cambridge University Press
,
Cambridge, UK
.
46.
Papka
,
S. D.
, and
Kyriakides
,
S.
,
1994
, “
In-Plane Compressive Response and Crushing of Honeycomb
,”
J. Mech. Phys. Solids
,
42
(
10
), pp.
1499
1532
. 10.1016/0022-5096(94)90085-X