VOLCO: A predictive model for 3D printed microarchitecture

Material extrusion additive manufacturing is widely used for porous sca ﬀ olds in which polymer ﬁ laments are extruded in the form of log-pile structures. These structures are typically designed with the assumption that ﬁ laments have a continuous cylindrical pro ﬁ le. However, as a ﬁ lament is extruded, it interacts with previously printed ﬁ laments (e.g. on lower 3D printed layers) and its geometry varies from the cylindrical form. No models currently exist that can predict this critical variation, which impacts ﬁ lament geometry, pore size and me- chanical properties. Therefore, expensive time-consuming trial-and-error approaches to sca ﬀ old design are currently necessary. Multiphysics models for material extrusion are extremely computationally-demanding and not feasible for the size-scales involved in sca ﬀ old structures. This paper presents a new computationally-e ﬃ cient method, called the VOLume COnserving model for 3D printing (VOLCO). The VOLCO model simulates material extrusion during manufacturing and generates a voxelised 3D-geometry-model of the predicted microarchitecture. The extrusion-deposition process is simulated in 3D as a ﬁ lament that elongates in the direction that the print-head travels. For each simulation step in the model, a set volume of new material is simulated at the end of the ﬁ lament. When previously 3D printed ﬁ laments obstruct the deposition of this new material, it is deposited into the nearest neighbouring voxels according to a minimum distance criterion. This leads to ﬁ lament spreading and widening. Experimental validation demonstrates the ability of VOLCO to simulate the geometry of 3D printed ﬁ laments. In addition, ﬁ nite element analysis (FEA) simulations utilising 3D-geometry-models generated by VOLCO de- monstrate its value and applicability for predicting mechanical properties. The presented method enables structures to be validated and optimised prior to manufacture. Potential future adaptations of the model and integration into 3D printing software are discussed.


Introduction
The material extrusion additive manufacture of complex structures such as lattices and porous scaffolds is an active research field [1,2]. 3D printing of tissue engineering scaffolds has reached a level of maturity where cell-laden hydrogels can be extruded together with biodegradable polymers for tissue regeneration in any shape [3,4]. Computer modelling and simulation can be used to enhance and optimise the 3D printing process. In 2009, Mironov et al. [5] identified a need for new computer models in their biofabrication review paper and stated, "We strongly believe that biofabrication from the start must be a predictable technology and built on predictable models and measurable parameters." But this issue has still not been addressed since Pati et al. [6], Paulsen and Miller [7] and Tang et al. [8] still highlight a need for mathematical modelling and computer simulations to support the design of tissue engineering constructs. Similarly, a 2017 review of additive manufacturing of lattice structures [2] identified several aspects of 3D printing that require new modelling capabilities.
There is growing interest in the use of computational modelling and simulation for biofabrication [9]. Therefore, porous scaffolds are used for demonstrations in this study. An important application of computational models is to enable validation and optimisation of scaffold microarchitecture a priori, as opposed to through a traditional experimental trial-and-error approach. Through a priori optimisations, Giannitelli et al. suggest that "finite element analysis has played a major role in the reduction of in vitro and in vivo experimental efforts" [9]. Indeed, a large number of studies have used FEA to validate or optimise scaffold structures including for mechanical properties [10,11], oxygen diffusion [12] and cell responses to external loads [13]. A key limitation of the above simulations is that the accuracy strongly depends on the geometric accuracy of the 3D-geometry-model being evaluated. Other complex structures such as lattices also need computational simulation of mechanical properties due to the impracticality of iterative experimentation in terms of cost and time [2].
At present, discrepancies between the designed 3D geometry and final printed structure are poorly understood [2,9]. For example, in many scaffolds fabricated by material extrusion additive manufacturing, filaments widen as they cross previously deposited perpendicular filaments [14][15][16][17]. This impacts the scaffold geometry and mechanical properties but no modelling techniques currently exists that can simulate this geometrical variation. Thus, computational analyses and optimisation of scaffolds prior to 3D printing is currently limited by a lack of accurate predictive geometric modelling techniques. Additional causes of variation to filament geometry include: changes in print-toolpath direction; filament interactions with the build platform; and filament interactions due to deposition directly on top of existing parallel filaments. A recent review of additive manufacturing of lattice structures [2] concluded that "geometrical discrepancy [between the designed and manufactured geometry] is a critical issue" and that "the manufacturing influence of additive manufacturing processes cannot be neglected by designers." They also report that it is "imperative to simulate the mechanical properties of lattice structures" and that "[manufacturing factors] need to be incorporated in to the geometric models." This clearly indicates the immediate need for new models to predict the microscale geometry of 3D printed structures.
Multiphysics simulations can give highly detailed insights into manufacturing processes. However, the following examples demonstrate the high computational requirements of simulations that consider heat transfer/convection, liquid/solid phase changes and material flow: • Even when simplifying the simulation to be 2D, a recent multiphysics study of heat transfer during welding required several days of high performance computing time [18].
• A 3D welding simulation study [19] required 1-3 days of computation time even with an order of magnitude coarser resolution (0.2-3.2 mm) than may be required for tissue engineering scaffold simulations (note that the number of 3D elements increases with resolution to the power of 3).
• A recent study for heat transfer in laser-based additive manufacture utilised a resolution of approximately 1 mm [20], over an order of magnitude more coarse than is required to effectively simulate a tissue engineering scaffold.
It is unfeasible to use multiphysics models to simulate filament extrusion of 3D printed scaffolds due to the high computational demands. This applies to any structure with geometric features spanning several size-scales. Therefore, a simplified modelling approach is required to enable feasible computational requirements for complex structures.
Simplified models for material extrusion additive manufacturing are reviewed in Table 1, but none can predict the 3D microarchitecture for complex structures with multi-directional interactions between individual filaments: • Studies with a resolution > 500 μm [21][22][23][24] consider geometry variation over a distance of several millimetres to centimetres. They allow useful simulation of large structures, but the variation of 3D geometry for individual filaments is out of scope -their fundamental modelling concepts are, by design, not applicable to individual filaments.
• Studies with a finer resolution (typically ≤25 μm) [35][36][37][38][39][40][41] focus on specific aspects of the material extrusion additive manufacturing process in order to ensure computational feasibility. Models have been developed for dynamic bond formation between parallel fused filaments [35,36], but they are not applicable to tissue engineering scaffolds because filaments in scaffolds have non-parallel orientations. Other studies considered individual stages of the material extrusion additive manufacturing process, including the feeding of input material [37], liquefied polymer dynamics in the melt chamber [37][38][39] and thermal history of an extruded filament [35,40]. Several models were integrated together in the thesis of Bellini [41] to simulate melt flow, nozzle extrusion and filament swelling/deposition. However, the computation-demands of their  [21] ≥ 500 μm Model to simulate variations in strut radii of lattice structures due to manufacturing process uncertainties [22] Adaptation of [21] to more accurately model structures over a larger scale [23] Model for the effect of variation in strut diameter on mechanical properties [24] Model for the effect of heating and cooling during manufacture on overall part distortion [25] ≈ 50 -500 μm Model air gaps and stair-step geometries of lattice structures [26] Print-toolpath strategies to optimise structural properties [27,28] Print-toolpath generation methods to improve the overall product quality (e.g. by avoiding large gaps between filaments) [29] Analytical model to predict the effects of the print-toolpath and filament orientation on mechanical properties [30][31][32][33][34] Analytic models for the surface geometry of 3D printed parts [35] ≤ 25 μm Dynamic model for bond formation between aligned neighbouring filaments [36] Improved [35] to include temperature-dependent material properties [37] Model considers the feeding of input material for material extrusion additive manufacturing [37][38][39] Models for liquefied polymer dynamics in the melt chamber [40] Model the thermal history of an extruded filament [41] Multi-part model including melt flow, nozzle extrusion and filament swelling/deposition but only suitable for 2D simulation due to computational demands THIS STUDY Model to simulate the 3D geometry of individual extruded filaments A. Gleadall et al.
This study presents a new VOLume COnserving model for extrusionbased 3D printing (VOLCO). The complex interaction between multiple extruded filaments is modelled based on a simple principle of conservation of volume. An accurate 3D-geometry-model of the predicted as-fabricated microarchitecture is generated by the model. No other models currently exist to achieve such predictions. The novel predictive model establishes the first step in bridging the gap between multiphysics models (which strive to capture the full range of factors that influence the fundamental 3D printing process) and computer programs currently available to 3D printer end-users (limited or no predictive capabilities). An application of porous scaffolds is used to demonstrate the model capabilities here because the individual extruded filaments in such scaffolds have highly variable microscale geometries. The applicability of the model to other fields is discussed.
The structure of the rest of this article is as follows: in the next section, the VOLCO model is described along with the materials and methods used in simulations and experiments; Section 3 presents the results of experimental validation and demonstrates predictive capabilities of the model; and Section 4 discusses the value, limitations and future opportunities for VOLCO.

Voxel modelling concept
The concept of the VOLCO model is to simulate the deposition of 3D printed filaments in a virtual 3D voxel environment. The model simulates a virtual filament that elongates in the direction that the print head moves in a series of steps (Fig. 1a). The key strength of VOLCO is Fig. 1. Overview of the VOLCO modelling concept: (a) shows a schematic representation of how filaments are simulated to elongate in the direction of the printhead travel in several simulation steps; (b) shows that when filament elongation is obstructed, the obstructed volume is reallocated according to a minimum distance criterion; (c) shows a 3D voxel model (used for simulations in this study) for the situation shown in (b); and (d) shows how voxels are assigned to contain either air or filament material during each simulation step.
to simulate how the newly 3D printed filament interacts with previously printed filaments -at present, there are no methods available for the simulation of filament interactions (except continuous linear stacking of parallel filaments). Fundamentally, VOLCO is based on a simple assumption of conservation of volume: if new material deposited at the end of the filament is blocked by previously deposited filaments (top part of Fig. 1b), the material is positioned as close as possible to the end of the filament in the virtual environment according to a minimum distance criterion (bottom part of Fig. 1b). This concept removes the need to consider complex effects of temperature-dependent viscoelasticity, fluid mechanics and heat transfer, which would render simulations of the size-scales undertaken in this paper computationally-unfeasible. Assumptions and limitations of this fundamental concept are considered in Section 2.3 and in the discussion. A 3D representation of the voxel modelling concept is shown in Fig. 1c, in which a filament is blocked by a previously 3D printed filament and therefore expands.

Model parameters
Parameters used in VOLCO are given in Table 2 and described in this section. The simulation step distance (d s ) is the length of elongation of the filament versus the previous simulation step, as shown in Fig. 1a. This is equal to the distanced travelled by the printhead between simulation steps. As d s approaches zero, the model approaches simulation of continuous extrusion (as opposed to discrete steps). As d s increases, computational demands reduce but the value must be set small enough to avoid the introduction of artifacts in the 3D-geometry-model. Practically, to ensure simulation accuracy, the value of the d s parameter is reduced during simulations until further reduction has negligible impact on model results. Opportunities for automated recommendation of appropriate parameter values are considered in the discussion.
The deposition volume parameter (V d ) refers to the volume of material that is deposited by the nozzle during each simulation step, as shown in Fig. 1a. For each simulation step, the filament elongates and occupies new previously-empty voxels. The total volume of newly-occupied voxels must equal the deposition volume (for conservation of volume). In order to identify which volumes are occupied with new material when the filament elongates, the deposition radius parameter (r d ) is used. Any voxels that lie within this radius from the point at the end of the filament are considered to be occupied by polymer, as shown in Fig. 1(b) and (d). During simulations, r d iteratively increases from zero in small steps (step size stated in Table 2), encompassing an everincreasing number of voxels, until the volume of newly-occupied voxels equals the required deposition volume. This gradually broadening search for new voxels to be occupied by the filament naturally achieves the modelling concept of positioning newly-deposited material according to a minimum distance criterion (when previously printed filaments obstruct the new filament).
The value of V d is directly derived from the 3D printer setup according to the 3D printer extrusion rate (v e ), which indicates the volume of material extruded per millimetre of printhead travel. The calculation of v e depends on the design of 3D printer and is typically detailed in the printer manual: • For printers that are fed material from a filament-reel, the extrusion rate is controlled by setting the rate at which the filament-reel is fed into the melt chamber.
• For screw-extrusion printers, the extrusion rate is controlled by setting the rate at which the screw mechanism turns.
• For syringe-based 3D printers with mechanical syringe-piston displacement control, the extrusion rate is controlled by setting the rate at which the piston is depressed.
• For syringe-based 3D printers with pneumatic syringe-piston control, the rate of material extrusion depends on many factors including needle size and viscosity. Pneumatic-controlled systems are out of scope for VOLCO, although the model can be potentially extended to consider such systems (see Discussion section).
Once the rate of material extrusion has been identified, the deposition volume for each simulation step is simply calculated as the rate of extrusion multiplied by the distance travelled by the nozzle during a simulation step, as given in Eq. (1): The 3D position of the printhead is derived from the print-toolpath or machine control code (GCODE). The filament is simulated immediately below the nozzle and elongation of the filament follows an identical path to the printhead. Voxels that are blocked by the printer nozzle in each simulation step are not considered for deposition since they are already occupied (by the nozzle). Similarly, no extrusion was permitted in voxels below the build platform.
None of the parameters in VOLCO need to be adjusted to fit experimental data. They are all either directly derived from the 3D printer setup or adjusted purely to ensure numerical accuracy of the model. The model was deliberately kept as fundamental and simple as possible here to allow flexible future adaptation; more advanced, enhanced versions of the model are possible as considered in the discussion. A voxel side length (l v ) of 25 μm was utilised for the simulations in this study, although smaller voxels with l v = 12.5 μm were also tested to check that simulation results remained similar, as discussed in Section 3.2.

Model assumptions
Key simplification and assumptions of the modelling concept are described below. They ensure feasible computational demands and avoid the need for experimentally-calibrated parameters. In this study, we present the fundamental, simplest version of the model and Table 2 Parameters used in the new modelling concept. No parameters were adjusted to fit experimental data. All parameters were either directly derived from the 3D printer setup or varied purely for purposes of ensuring computation accuracy and efficiency. demonstrate that even the simplest form of the model has value in simulating 3D printed structures. Many of the simplifications and assumptions listed below may be addressed in extended versions of the model as required for specific 3D printing applications. For example, a more advanced version of the model specifically for pneumatic 3D printing hydrogels at room temperature may include the capability to relate extrusion rate to syringe pressure but has no need for thermal considerations. Simplification and assumptions: • Gravity is not considered ○ Justification: gravitational forces are low compared to those that can be exerted by the nozzle; filaments on lower layers evidently surpass gravitational forces to prevent scaffold collapse; the results section of this study show several-fold increases in compressive modulus due to filaments interacting with each other, outweighing any minor effects of gravity.
• Material does not flow after it has been positioned in the model • Justification: this avoids unfeasible computation demands associated with multiphysics material flow simulations -without this simplification the model may need to simulate temperature, viscoelasticity and phase changes, which have extremely high computational demands.
• Thermal contraction is not modelled • Justification: VOLCO effectively considers all thermal contraction to be completed by the end of each simulation step. Thermal contraction has a minor effect on geometry (several percent) by comparison to filaments obstructing each other (up to 66% change in experimental filament width in Section 3.1); thermal simulations including convection have extremely high computational demands.
• Material is extruded from the nozzle at a constant rate • Justification: Molten polymer is considered to be an incompressible fluid; mechanically-driven extrusion forces can overcome forces resisting polymer extrusion.
• Obstructed material is relocated according to a minimum distance criterion • Justification: An alternative assumption was initially considered whereby filaments are represented as a series of cylindrical segments of varying radius. However, such as assumption runs into complications when the printhead changes direction (jagged edges would be simulated). The minimum distance criterion was considered to be a flexible fundamental approach.

Software simulation routine
The simulation routine is presented in the flowchart in Fig. 2. The following steps were followed in MATLAB 2016b to generate 3D-geometry-models of the predicted polymer scaffold structure: 1 Import a list of printer nozzle coordinates (the print-toolpath) and import the simulation setup parameters in Table 2 2 Identify the Cartesian location of the nozzle, and therefore the direction of elongation of the filament, for each simulation step 3 Simulate elongation of the polymer filaments in the voxel 3D matrix for the print-toolpath. Each simulation step includes the following sub-steps: a Identify the coordinates considered to be the end of the filament b Set r d to equal the nominal filament radius and set voxels within this radius to contain polymer c Determine the increase in the total volume of polymer-voxels in the modelling environment d If this increase in volume is less than V d , some of the targeted voxels already contained polymer (e.g. a previous simulated filament). Therefore, iteratively increase r d and set voxels within this radius to contain polymer until the target deposition volume is achieved. 4 Output a 3D-geometry-model of the simulated structure (STL file) The following steps were taken to model scaffold compressive modulus 1 Import the 3D model into the FEA software 2 Setup boundary conditions and other FEA parameters 3 Run simulation 4 Determine compressive modulus from reaction forces

Experimental scaffold design and fabrication
Tissue engineering scaffolds are frequently fabricated by 3D printing layers of polymer filaments in "log-pile" structures, as shown in the idealised CAD representation in Fig. 3. The layer thickness (LT) indicates the vertical distance moved up by the nozzle between layers. It is often varied when 3D printing scaffolds because it can be used to control mechanical properties and porosity. Therefore, in this paper, the model is validated by varying LT and comparing modelled and experimental results for the effects the variation has on compressive modulus, filament geometry and pore fraction.
All scaffolds were 3D printed on an Orion Delta 3D printer in poly lactide (NatureWorks ® polylactide 4043D) using a 300 μm nozzle diameter at a temperature of 210°C. The rate of filament extrusion was set to achieve a nominal 300 μm filament diameter and the nozzle travel rate was 600 mm min −1 . The height of the nozzle above the build platform for the first printed layer of all scaffolds was set to 175 μm. The 3D printer control code (GCODE) for each scaffold was generated using a custom-developed Visual Basic routine. Lay-down patterns of 0/90°and 0/45/90/135°have both been studied in the literature [15,42] and are therefore used to validate the model in this work as detailed in Table 3. To test the effect of layer thickness on pore fraction, 12 mm wide × 12 mm deep × 20-layer scaffolds were 3D printed with LT = 75, 100, and 125 μm and a 0/90°l ay-down pattern. Small layer thicknesses, in relation to nozzle diameter, were used to ensure a high interaction between filaments stacked on top of each other; this forced the filaments to spread sideways (to have a flatter, wider cross section as opposed to circular) and therefore provided an effective challenge for the model. To test the effect of LT on compressive modulus, scaffolds were printed as 18 mm wide × 18 mm deep × 9 mm tall scaffolds with LT = 100, 150, 200 and 250 μm and a 0/45/90/135°lay-down pattern. For each printed scaffold, nine 4.5 mm × 4.5 mm × 9 mm samples were cut from the central section with a knife blade in a press, using a new blade for each sample.

Microscopy for porosity and filament width
A Zeiss Stemi 2000-C microscope with a Schott S40-10D ring light was used for microscopy. ImageJ 1.49v (National Institutes of Health) was used to analyse pore fraction and filament width. Filament width refers to the width of a filament from a top-down view and must increase when filaments are more closely stacked on top of each other because the filaments are impeded from taking a cylindrical form and must spread and widen due to interactions with filaments on lower layers. This is discussed in more detail in Section 3.1. To determine pore fraction, images with a pixel size of 5.4 μm were binarised according to greyscale value (0-255) using a threshold of 105. The four central internal pores were analysed for each scaffold to determine average volume fraction. Filament width was measured from the microscope images, midway between crossover points with perpendicular filaments. Four filament width measurements were taken for each scaffold.

Mechanical compression testing
Compressive tests were conducted on an Instron 5969 machine with a 5 kN load cell according to ASTM standard D695. Unconstrained samples were compressed between two steel plates at a rate of 1 mm min −1 to 50% strain. Five samples (4.5 mm × 4.5 mm × 9 mm) were tested for each layer thickness. Scaffolds were compressed in the build-direction and the initial linear-elastic compression phase was used to calculate compressive modulus.

Finite element analysis
Linear elastic FEA simulations were completed using the commercial software package ABAQUS/CAE 6.12-3 (Dassault Systèmes). Fig. 4 shows the mesh and boundary conditions, which were as follows: the lower surface of the 3D model was kept planar and prevented from zdisplacement; the top surface of the 3D model was kept planar and a negative z-displacement was applied to achieve 2% strain; the sides of the model were not constrained since they directly represent the unconstrained external surface of the experimental sample. Compressive modulus was calculated based on the z-reaction forces and unit cell area. The elastic modulus parameter was set to 2.45 GPa, which was determined through calibration with the experimental sample for LT = 250 μm, as discussed in the results section. This value is reasonable in comparison to values in the literature that range from < 1 GPa to > 7 GPa [43][44][45][46][47]. This single calibrated value was used in simulations for all values of LT. Poisson's ratio was set to 0.36 [48,49].
in which a t is the experimental measurement, f t is the model prediction and n is the number of samples. Fig. 3. CAD representation of the idealised log-pile structure illustrates a 0°/90°l ay-down pattern and the LT parameter that was varied in scaffolds for experimental validation of the model. Note that in this idealised representation, interactions between filaments (e.g. filament spreading/bulging) are not indicated -the lack of any existing method to predict such interactions gives justification to and highlights the importance of the presented model.

Results
The results are split into two sections. The first section validates the ability of VOLCO to simulate a scaffold in which 3D printed filaments are obstructed by previously printed filaments on lower layers. Due to this obstruction, filaments become wider and flatter. For validation, the 3D-geometry-models generated by VOLCO are compared to experimental results. The second section demonstrates the ability of VOLCO to support the prediction of scaffold mechanical properties. The 3Dgeometry-models generated by VOLCO are imported into FEA software and used to predict compressive modulus. Experimental samples are 3D printed and tested in order to validate the FEA predictions.

Modelling filament geometry
Twenty-layer scaffolds with square pores from top-to-bottom were 3D printed (Fig. 5a) to test the ability of the VOLCO modelling concept to simulate the 3D geometry (Fig. 5b). Three different printing setups were used in which LT = 75, 100 and 125 μm. The top-view "Experimental macrograph" images in Fig. 5c show that filaments widened (indicated by the arrows to either side of hashed regions) as layer thickness reduced from left to right in the figure. This widening was due to filaments being obstructed by previously-printed filaments (on lower layers) as shown in the idealised schematic at the top of Fig. 5c. As LT reduced, the filaments were more obstructed and therefore widened to a greater degree. The VOLCO model effectively captured this trend, as can be seen in the top-view images of 3D-geometry-models at the bottom of Fig. 5c. Filament width was measured from the experimental and model images and the results are given in Fig. 6a. The mean percentage error between the model predictions and experimental values was 5.0% (3 sets of n = 4). This is a good level of predictive accuracy given that experimental filament width increased by 66%, from 305 to 505 μm experimentally, as LT reduced from 125 to 75 μm. For the sample with LT = 75 μm, filaments were ≈505 μm wide by 150 μm tall (150 μm = 2 × 75 μm layer thickness because stacked filaments had the same orientation every two layers, as shown in Fig. 3), which highlights the ability of VOLCO to simulate filaments with a non-circular cross-section. Fig. 6b shows the pore fraction measured experimentally (65.2%-77.0%) and from modelling results (68.9%-79.1%). The model simulated the trend effectively but predicted a slightly greater pore fraction for all samples (MPE = 4.6%; 3 sets of n = 4). This may have been due to experimental variation between 3D printed layers: any misalignment between filaments on different layers reduced the visible pore area in top-view images.

Compressive modulus prediction
This section validates the ability of VOLCO to support FEA simulations of scaffold compressive modulus. Four different scaffolds were 3D printed with LT = 100, 150, 200 and 250 μm and tested for compressive modulus. The sample shown in Fig. 7(a) and (b) had LT = 200 μm. VOLCO was used to general 3D-geometry-models of the scaffolds (Fig. 7d). These 3D models were compared to the experimental scaffold geometries and were used in FEA simulations to predict compressive moduli. The predicted and experimental values of compressive moduli were compared in order to validate the model. To demonstrate the value of the VOLCO model, FEA simulations were also conducted using an idealised CAD model of the scaffold, in which filaments had a constant diameter equal to the nozzle diameter (Fig. 7c)this is a simple approach that researchers can currently use to generate 3D models of a scaffolds in standard 3D CAD packages. The key difference between the VOLCO model geometries and the idealised CAD model geometries was that the idealised CAD models did not and could not have the scope to consider filament widening due to obstruction by previously printed filaments. Fig. 8a shows that filaments had a constant diameter in the experimental sample with LT = 250 μm. This indicates that the filaments were not greatly obstructed by other filaments (on lower layers), which is expected because the nozzle raised by 250 μm between layers: only 50 μm less than the nominal filament diameter (300 μm). In contrast, in the sample with LT = 100 μm, the nozzle only raised by 33% of the nominal filament diameter between layers and therefore the filaments were obstructed and widened at cross-over points with lower-layer filaments, as can be seen in Fig. 8b. The widening of filaments at crossover points resulted in compressive modulus increasing more than 5.5× as layer thickness reduced from 250 to 100 μm, as shown in Fig. 9 (black triangles). The scatter of the individual experimental results (n = 5) was < ± 8% for the samples with LT = 150-250 μm and ± 20.5% for the 100 μm samples.

Proposed model FEA results
The VOLCO model effectively simulated filaments that either maintained a constant width (Fig. 8a) or widened at the crossover points (Fig. 8b), as can be seen in the comparisons between the modelled and experimental geometries. FEA simulations were conducted to predict the effect of this filament widening on compressive modulus. The 3D-geometry-models generated by VOLCO for four different layer thicknesses were used in FEA simulations to determine the compressive modulus of each scaffold. The values of compressive moduli determined through FEA simulations are shown in Fig. 9 (hollow squares) and followed a similar trend to experimental results (solid black triangles). Fig. 10 shows the FEA von Mises stress distribution for the four scaffold models. It can be seen that compressive forces were sustained by areas of highly stressed polymer at points where filaments crossed over. The stresses were greatest in these areas because there was a continuous column of polymer from top-to-bottom of the scaffold (in the direction of compression) that resisted compression. In contrast, isolated sections of filaments elsewhere in the scaffold were surrounded by air and could therefore be displaced relatively easily and experienced low levels of stress (as can be seen in the Fig. 10). As LT reduced, filaments widened to a greater degree at crossover points (visible in Figs. 8b and Figure 10) and therefore compressive modulus increased because a wider column of polymer resisted compression compared to samples with a larger layer thickness.
The bulk material elastic modulus parameter in FEA simulations was calibrated to the sample with LT = 250 μm because the value is highly dependent on experimental factors such as testing strain rate, stress-strain-curve analysis procedure, polymer molecular weight, compressive versus tensile testing, and many other factors. A calibrated elastic modulus value of 2.45 GPa, which is reasonable for poly lactide [43][44][45][46][47], was used in all four FEA simulations. Even with the FEA calibration step, the experimental results present a fair and difficult challenge for the model since there was a large 5.5x increase in experimental compressive modulus for samples with 100 μm versus 250 μm layer thickness. This change was accurately predicted in the FEA simulations as 5.9× (MPE = 6.6%; n = 5). The mean percentage error versus experimental compressive modulus for all samples was 11% (3 sets of n = 5; 250 μm sample was excluded since it was used for calibration).

Idealised CAD model FEA results
Four idealised CAD models were generated in SolidWorks 2016 for the four different LT values and used in FEA simulations for comparison to the VOLCO model. In contrast to 3D-geometry-models generated by the VOLCO model, the idealised CAD models did not consider the effects of filaments widening, therefore FEA predictions obtained using them ( Fig. 9 hollow circles) did not effectively capture the experimental trend of increasing compressive modulus as layer thickness reduced ( Fig. 9 solid black triangles). The mean percentage error versus experimental compressive modulus for all samples was 42.9% (3 sets of n = 5; 250 μm sample was excluded since it was used for calibration). These results highlight the importance of the VOLCO modelling concept being able to simulate filaments with varying cross-sectional geometries. The calibrated elastic modulus value used in FEA for all idealised CAD models was 2.64 GPa.
Each voxel simulation of polymer deposition took less than 5 min and each FEA simulation less than 10 min on a desktop computer. FEA simulations were repeated with a smaller unit cell (one ninth of the size) and planar boundary conditions in the X and Y directions. These simulation results remained within 0.2%-2.8% of those for the original unit cell whilst simulations took just 1 min. Additionally, to test the sensitivity of the model, the voxel side length was halved to l v = 12.5 μm. The simulations took an order of magnitude longer but the change had little impact on compressive modulus (< 1% change).

Validated capabilities of the VOLCO model
Comparisons of the simulated filament geometries to experimental filament geometries in Figs. 5, 6 and 8 demonstrate the ability of the VOLCO model to simulate substantial changes in filament geometry ranging from a nominally circular cross-section of 300 μm diameter to wide flat cross-sections approximately 500 μm wide by 150 μm tall. For stacked filaments oriented in the same direction, it is possibly to analytically determine the approximate width of filaments based on their constrained height (to match the cross-sectional area of a circular filament). However, analytical methods become complicated for interactions between just two filaments crossing perpendicularly or at an angle of 45°, and analytic modelling of interactions for a whole 3D printed structure is unfeasible. The fundamental, basic nature of VOLCO means it can simulate a large range of filament geometries and therefore has applicability to a wide range of different structures: it is demonstrated for log-pile structures but is anticipated to be applicable to other types, such as commonly 3D printed shell + infill parts. The authors are unaware of any models or any previous publications presenting methods to simulate the varying geometry of individual 3D printed filaments (aside from models for simple interactions between isolated filaments, which have vastly different scope and relevance in comparison to the method presented in this study for full 3D printed structures).
The FEA simulation results demonstrate the value of the 3D-geometry-models generated by VOLCO. These 3D-geometry-models are predictions of the as-fabricated 3D printed microarchitecture. In the present study, these models were used in FEA simulations to predict compressive modulus of the scaffolds. However, the 3D-geometrymodels could be used for a wide range of other simulations or analyses, for example to simulate fluid flow or molecular diffusion through a scaffold or to analyse the distribution of pores within the scaffold. By changing the layer thickness of the 3D printed scaffold, compressive modulus was shown to increase 5.5x for experimental samples. The VOLCO model simulated an increase of 5.9x through FEA, which demonstrates a remarkably accurate prediction for such a large change in compressive modulus with no calibration or fitting of any model parameters (the calibration of the elastic modulus FEA parameter in Section 3.2.2 bears no impact on these 5.5× and 5.9× increases since they are relative increases). No other models exist that can predict the relationship between compressive modulus and layer thickness demonstrated here. The results for a simple CAD representation of the structure, which did not account for filament widening, show that it could not capture the experimental trends. This highlights the importance of VOLCO's ability to simulate variation in filament cross-  section. Although layer thickness was varied in this study to validate the model, it could also simulate a wide range of other factors that have been studied experimentally including nozzle size, filament spacing and lay down pattern.

Computational demands and 3D printing software integration
The modelling simulations in this study took 3-50 min. Simulations with greater interactions between filaments (smaller layer thicknesses) took longer because a greater number of iterations (increasing r d ) were required in each simulation step to achieve wider filaments. The authors anticipate that simulation times could be reduced to under a minute with program optimisations for speed. If the method was adapted to consider a periodic unit cell (particularly relevant given the repeating nature of many 3D printed scaffolds/lattices), it is likely that effectively-instantaneous predictions of 3D-geometry could be integrated into 3D printer software.
In most 3D printing software, users are currently presented with visualisations of the print-toolpath before the part is 3D printed. This enables the user to identify problems such as undesired gaps, which they can correct by either revising the 3D geometric design or changing the print-toolpath-generation parameters. If VOLCO was integrated into 3D printing software, the 3D printing process could be simulated in more detail to identify potential issues such as voids or regions of As LT reduced, filaments widened to a greater degree at crossover points with previously printed filaments due to greater obstruction. This caused larger areas of polymer to sustain the compressive force and resulted in an increase in compressive modulus. excessive filament overlapping. Furthermore, automated predictions of pore size, pore fraction and mechanical properties could be generated (some FEA simulations took < 1 min in the present study).
As an alternative to presenting designers with modelling results, VOLCO could be implemented within 3D printer software to automatically optimise the print-toolpath or print parameters. To achieve this, feedback of the model results within the 3D printer program would be utilised to iteratively update the printing setup. In such an approach, the user may specify design requirements (e.g. pore size, compressive modulus, pore fraction) as opposed to printer setup parameters (e.g. extrusion temperature, filament spacing, layer thickness).
Although some model parameters were iteratively adjusted in this study to ensure numerical accuracy and avoid model artefacts, it would be possible to develop algorithms for automatic identification of model parameters. Preliminary assessment of the print-toolpath and 3D printer setup would determine the model size, the level of interaction between filaments and the extent of filament widening. This would enable software to identify effective parameter values and present the user with a choice of simulation accuracies and respective simulation time estimates.
Lattice generation software can automatically create a 3D lattice structure optimised for requirements such as mechanical properties [50]. The software generates a 3D-geometry-model of the intended structure, which is imported into 3D printing software for print-toolpath generation. During generation of the print-toolpath and during fabrication, discrepancies between the intended geometry and the final 3D printed geometry are introduced. VOLCO could be integrated into lattice generation software to optimise the final fabricated lattice geometry, rather than the conceptual 3D geometry. Park and Rosen [25] recently proposed a method to simulate as-fabricated lattice geometries predicted for theoretical designs. Significant improvements in predictions for mechanical properties were achieved which demonstrates the importance of modelling discrepancies between designed and manufactured 3D geometries. In their method, filaments were assumed to have a constant elliptic-rectangular cross section. Therefore, it is a completely different concept to VOLCO and cannot replicate results in the present study because, similar to the idealised CAD model in Section 3.2.3, it does not consider geometric variation of individual filaments.

Model limitations and potential extension
The VOLCO model is based on a simple assumption that as a filament is 3D printed, the rate of material deposition remains constant. In other words, if filament extrusion is partly obstructed by previously printed filaments, material is still extruded from the nozzle. In the model, the newly extruded material for each simulation step is positioned according to a minimum distance criterion from a point at the end of the growing filament. This geometric modelling concept is intended to be as fundamental and simple as possible; it allows for a wide range of future adaptations. The modelling concept can be readily extended to consider extra factors for specific deposition conditions or for more detailed analysis. For example: gravitational effects may be added for scaffolds where filaments bridge large gaps; process imperfections may be considered such as fluctuations in deposition rate or nozzle position; and melt flow or cooling effects could be incorporated for microstructure analysis.
In many cases, extensions to VOLCO would have negligible impact on computational demands; for example, to simulate the effect of experimental fluctuations in deposition rate or nozzle position. In terms of implementation, many extensions would require very little reprogramming of the physical software code; for example, to consider a random misalignment error for each layer, just a single extra line of program code would be required to offset the coordinate list for each new layer.
A limitation of VOLCO is that the extrusion rate of material was considered to be constant, which makes the model applicable for 3D printers with a mechanically-driven syringe-plungers (e.g. 3DYNAMIC SYSTEMS OMEGA) or screw-extruders (e.g. RegenHU 3DDiscovery), or 3D printers fed by a polymer filament reel (e.g. REGEMAT3D). For pneumatic printers (e.g. BioBots), it would be interesting to extend the model to include a term to represent pressure in the syringe and relate it to a varying extrusion rate of new material. It may be possible to achieve such an extension with minor increases in computation-demands.
Practically, it may be useful to use VOLCO in two phases: 1) uncalibrated -for quick analyses of scaffold geometries and the relative effects of design choices on factors such as compressive modulus; 2) calibrated -(e.g. elastic modulus FEA parameter calibration) for predictive use with a higher level of confidence or to optimise a design for specific requirements.

Applicability to tissue engineering
For tissue engineering scaffolds, porosity has been shown to be an important factor for cell ingrowth [51], vascularisation [52], nutrient and oxygen exchange [12], and mechanical properties [14]. The ability to predict pore size and pore fraction enables scaffold designs and printer settings to be optimised for desired characteristics with reduced experimental trials. It addresses the gap that currently exists between designed and manufactured scaffold microarchitectures and improves the link between experimental and computational models that is currently lacking [7,9].
As discussed above, two key potential uses of the model are as follows: • Validation and optimisation of a scaffold design prior to printing • Automated optimisation of printer parameters and the print-toolpath through feedback of model results within 3D printer software Validation of a scaffold design prior to printing enables better design optimisation during the design process and is therefore hugely beneficial to reduce experimental iteration [6,7]. Currently, characterisation of scaffold design is often achieved through 3D printing > microCT scanning > microCT 3D-geometry-model generation > microCT 3D-geometry-model analyses. This is a time-consuming process that requires several days of expensive resources. In contrast, the VOLCO model can be used a priori to optimise or validate a construct design in a few minutes. Although final scaffold geometries should still be characterised experimentally, the use of the model may save days or weeks of trial-and-error printing. Furthermore, the 3Dgeometry-models generated by VOLCO enable a very wide range of potential analyses including pore fraction, pore shape/size, open/closed porosity, oxygen diffusion, mechanical properties and computational fluid dynamics. Future work could also translate the model to biomaterials including cell-laden hydrogels.
In situations where scaffold designs must be optimised for several factors (e.g. mechanical properties, degradation profile and biological response), the number of experiments required for trial-and-error optimisation may be impractical. New design strategies to determine the optimal trade-off between the conflicting requirements are necessary [9]. The 3D geometry predictions of the VOLCO model may enable topological optimisation of multiple parameters due to low computational demands. The improved predictions achieved by the model may enable more accurate computational modelling of biomechanical forces in cells or biological responses, and thus support important efforts to reduce the number of animals required for in vivo studies [8].
With effective software-implementation, model simulations and analysis can be conducted with little or no modelling expertise. As clinical readiness increases for tissue engineering scaffolds, VOLCO may have potential value in the regulatory process as a validation tool for quality control checks of medical implants prior to printing; to identify any regions of closed porosity, for example.

Applicability outside of tissue engineering
Although this study demonstrated VOLCO for tissue engineering scaffolds, it is potentially suitable to a wide range of 3D printing applications. For most 3D printed parts, filaments are printed closer together than in log-pile tissue engineering scaffolds. The validation results in Section 3.1 show that the modelling concept can simulate the geometry of filaments that are closely stacked. Even coarse simulations with very low computational demands (e.g. a large voxel side length of l v = 50 μm) could prove highly useful in previewing parts prior to printing. In the authors' opinion, such simulations could potentially be conducted in less time than that taken by print-toolpath generation algorithms. Many aspects discussed in relation to tissue engineering above directly translate to other applications; for example where infill is used, the new modelling concept could enable the user to specify infill mechanical properties rather than infill density. Furthermore, the model may be translated to other processes; particularly those based on volumetric deposition such as inkjet 3D printing, for which the geometry of each droplet would be simulated.

Conclusions
This paper presented a new modelling technique, called VOLCO, to predict the microarchitecture of tissue engineering scaffolds by considering interactions between 3D printed filaments. It simulated filaments being 3D printed and generated 3D-geometry-models of the predicted structures. The voxel-based technique successfully simulated filament geometries ranging from nominally circular cross-sections (300 μm diameter) to wide flat cross-sections (505 μm wide × 150 μm tall). Filament width and pore fraction were predicted with mean percentage errors of 5.0% and 4.6%, respectively. It was also shown that the 3D-geometry-models generated by the VOLCO model can be used to predict compressive modulus through FEA simulations. During validation experiments, compressive modulus was found to increase 5.5 fold experimentally as the 3D printer layer thickness reduced from 250 to 100 μm. VOLCO predicted this increase to be 5.9 fold, which is extremely accurate for such a large modulus range and given that the prediction was achieved with no fitting of any model parameters. No other models or methods have been reported that can generate such predictions.
VOLCO has low computation demands. Generation of 3D-geometrymodels that were used in FEA simulations took less than 5 min. For the repeating structures often used in tissue engineering scaffolds, effectively-instantaneous predictions are a realistic prospect if a repeating unit cell were utilised in simulations. The model is therefore highly suited for integration into 3D printer software, both for tissue engineering and potentially more general 3D printing applications.
The new modelling technique is a powerful tool that can be used to validate the geometry and properties of tissue engineering scaffolds before manufacture. It is hoped that it forms the foundation for a wide range of future adaptations that extend modelling capabilities for more detailed, application-specific analyses. Further research into computational modelling and simulation is critical to support future advancement of the 3D printing field, in particular to further understand or predict the as-fabricated structure and properties of 3D printed constructs.