A discrete lattice model for assessment of buildability performance of 3D‐printed concrete

In this work, the lattice model is applied to study the printing process and quantify the buildability (i.e., the maximum height that can be printed) for 3D concrete printing (3DCP). The model simulates structural failure by incorporating an element birth technique, time‐dependent stiffness and strength, printing velocity, non‐uniform gravitational load, localized damage, and spatial variation of the printed object. The model can reproduce the plastic collapse failure modes reported in the literature. In this research, three main contributions for 3DCP modeling work can be found. A new failure criterion is proposed and adopted to improve the estimation of critical printing height; the element birth technique is utilized to mimic the continuous printing process and study the impact of non‐uniform gravitational load; variability of a printed structure is modeled through the inclusion of disorder during mesh generation and Gaussian distributions of material properties. Using this model, parametric analyses on non‐uniform gravitational load and material variation are conducted to assess their impact on the failure–deformation response and the critical printing height. Finally, the model is validated by comparison with two 3D printing experiments from the literature. The proposed lattice model can reproduce the correct failure‐deformation modes of two types of structures commonly used for buildability quantification: A 3D‐printed hollow cylinder and a square wall layout. Lattice modeling of the square structure yields a relative difference of around 10% with the experimental printing height. For the cylinder structure, the predicted radial deformation and corresponding height show good agreement with the experimental data; the model yields a 41.38% overprediction of the total number of printing layers, compared with the experimental data. Possible reasons for the quantitative discrepancy are discussed.

. This additive manufacturing technology can transform digital data from a computer model to a physical product (Gosselin et al., 2016). Not only can 3DCP do away with the need for conventional molds, thereby reducing postconstruction waste, it can also considerably accelerate the manufacturing process and yield geometrically complex, non-standard structural elements instead of rectilinear shapes (Bos et al., 2016;Paul et al., 2018;Tay et al., 2017).
A layer-wise extrusion approach is generally adopted for 3DCP. Specifically, the printing materials are extruded from the nozzle to build a computer-designed object using a robotic arm. After deposition, the material should be strong enough to retain the shape and avoid collapse under self-weight and the gravitational load from subsequent printing layers. This is defined as buildability performance.
The buildability performance is codetermined by the printing process and material characteristics. Experimental investigations on various materials, including sulfur-based cement (Khalil et al., 2017), limestonecalcined clay cement (Chen et al., 2019;Chen, Figueiredo, et al., 2020), calcium aluminate cement (Shakor et al., 2017), geopolymer D. -W. Zhang, Wang, et al., 2018), and magnesium potassium phosphate cement blended with fly ash (Weng et al., 2019), have been carried out to assess their feasibility for 3DCP. Furthermore, some researchers have also studied the influence of environmental conditions (e.g., relative humidity and temperature) and printing parameters (e.g., time interval, printing nozzle height and head speed) on the final quality and interlayer bond strength of the printed object (Panda, Mohamed, et al., 2019;Sanjayan et al., 2018;Tay et al., 2019;. In spite of the rapid development of 3D printing technology, not many investigations concern the relationship between the buildability performance and printing parameters, such as the printing velocity, time-dependent material stiffness and strength, and randomly distributed material properties.
In general, most printing parameters are determined through a series of trial-and-error experiments, whereby it remains unclear whether optimized parameters have been derived for the target object. Besides, experiments are time and resource-intensive, especially for large structures such as those used in real-life engineering and architectural applications. Compared with experiments, numerical modeling can provide a cheaper, faster, and more easily controlled alternative.
However, due to the novelty and complexity of 3DCP, developing a modeling approach for the printing process presents numerous challenges. In general, the major steps for 3DCP can be divided into two processes: (1) transport of printable material to the nozzle; (2) deposition of the printable material, accompanied by self-weight and gravitational load from the subsequent printing layers. A schematic diagram can be found in Figure 1. During the printing process, two performance metrics, namely, pumpability and buildability, are utilized for quantifying the printing process. Regarding the pumpability, the material should flow like a fluid to be transferred in the hose and avoid blockage. The semi-solid characteristics such as apparent viscosity and dynamic yield stress determine the pumpability. Thus, the rheology theory or fluid mechanics approaches such as computational fluid dynamics are usually adopted to quantify the pumpability (Jeong et al., 2019;Khan et al., 2020;Perrot et al., 2016;Roussel, 2018;Wangler et al., 2016) or predict the cross-sectional shape of 3D printed segments (Comminal et al., 2020). On the other hand, solid mechanics approaches and rheology theory are adopted for buildability quantification as described in Figure 1.
Among the limited approaches for 3DCP (listed in Table 1), quantifying the buildability performance is mostly done based on the material rheology (Kruger, Cho, et al., 2019;Perrot et al., 2016;Roussel, 2018;Wangler et al., 2016). Rheological models accounting for flocculationinduced thixotropy and chemical phenomena generally predict buildability performance based on the static yield stress. Notable contributions include the analytical model of Roussel (2018), the experimentally validated lower bound analytical model of , and empirical models considering linear or exponentially decaying curing function (Perrot et al., 2016;Wangler et al., 2016). These models consider the printable material as a non-Newtonian fluid and evaluate the buildability performance based on rheological theory. However, in reality, the printable material should have a high yield stress to retain its shape and avoid collapse after extrusion. Thus, printable materials are at rest during most of the printing process. Consequently, their elastic-plastic behavior is more significant than their visco-plastic behavior (Roussel, 2018). Therefore, solid mechanics might be more suitable for buildability quantification. The most advanced examples are the mechanistic model of Suiker et al. (2018) and the finite element model (FEM) of Wolfs et al. (2018). The mechanistic model incorporating time-dependent material properties, printing velocity, and geometrical features can consider the two possible failure mechanisms of a wall structure in 3DCP: elastic buckling and plastic collapse. However, this mechanistic model is limited to the wallshaped geometry. In order to further investigate the structural failure in the experiment, an FEM model based on an age-dependent Mohr-Coulomb failure criterion and linear strain-stress relation was developed by Wolfs et al. (2018) using ABAQUS. This numerical model is able to qualitatively reproduce a correct failure-deformation mode for F I G U R E 1 Processing steps and analytical methods for 3D concrete printing (3DCP) in different stages the hollow cylinder structure with plastic collapse failure mode, but the quantitative agreement with experimental data still needs to be improved. More importantly, this model simulates the printing process through a layer-bylayer extrusion process, and the system state is regarded as having failed as soon as any point in the printed object reaches the material yield strength. Therefore, this numerical model might underestimate the loading capacity of the printed system since the material strength is exceeded only locally, while in reality, stress redistribution may occur. Furthermore, 3DCP is a continuous printing process, and the gravitational load is applied to the system in a gradual fashion. This kind of loading leads to the non-uniform stress distribution and threatens the structure buildability. In addition, the printing material undergoes an extrusion process, which may heighten material heterogeneity. Combined with the non-uniform gravitational load, this heterogeneity may also play a significant role in the occurrence of localized damage, thereby affecting the structural buildability. However, at present, no approaches considering these effects have been put forward in the published literature.
To investigate the failure mechanism in 3DCP, developing reliable tools is essential for buildability quantification. Given the aforementioned reasons, the solid mechanics method is adopted for buildability quantification in this study. Since the printable materials generally have the yield stresses as high as a couple thousand pascal (Pa) after deposition, the gradual increasing self-weight can ultimately result in cracking or localized damage similar to that observed in fragile solids (Roussel, 2006(Roussel, , 2018. This localized damage affects the stress redistribution and influences printing characteristics. Considering the need to account for localized damage, a lattice model that is suitable for fracture analysis is extended to quantify the buildability for 3DCP. In this numerical model, a new failure criterion is proposed and adopted for buildability quantification. Considering the continuous printing process, the element birth technique is adopted and utilized for investigating the influence of non-uniform gravitational load. Compared to the previous homogeneous models for 3DCP, variability of a printed structure is simulated through the inclusion of disorder through the mesh generation and using Gaussian distributions to assign the breaking parameters of the lattice elements. Through this proposed model, the buildability performance of a 3D printed object is quantified and the printing characteristics, namely, the critical printing height, radial deformation and corresponding height are predicted. This article is organized as follows. Section 2.1 gives an overview of model development. The model discretization, the failure criterion, loading condition and time-dependent material strength and stiffness, the plastic collapse typical failure mode and the model applicability are described in Sections 2.2 to 2.10. After the model synopsis, validation efforts are conducted by employing two 3D printing experiments with different geometries. In Section 3.1, computational uniaxial compression tests are first performed to calibrate the time-dependent material properties. To further investigate the influence of material and printing parameters, parametric analyses including the non-uniform gravitational load and randomly distributed material properties are performed to study their impact on printing characteristics. Then, based on the calibrated linear stress-strain relations and parametric analyses, two experimental campaigns are simulated through this proposed model for buildability quantification in Section 3.3. Finally, conclusions are drawn in Section 4.

Model overview
An overview of the numerical approach based on the lattice model is presented in Figure 2, which comprises six branches. In the flowchart, the k and K stand for the local and global stiffness matrix, respectively. T is the transformation matrix used to relate the local and global coordinate systems. L and D are the global load vector and displacement vector, respectively. F I G U R E 2 A general procedure for simulating the concrete 3D printing process using the lattice model Branch A: Establish a numerical model for 3DCP. The details are given in Sections 2.2 and 2.3.
Branch B: Assemble the global stiffness matrix K, the load vector L and displacement vector D based on the current printing time. When reaching the predetermined printing time, the corresponding printing pieces are activated and the inner elements within these activated printing pieces are then assigned with age-dependent stiffness and strength.
Branch C: Calculate the nodal force. After the printing pieces activation, the self-weight of these parts and the boundary conditions (as indicated in Section 2.5) are applied to the current system. Then, a conjugate gradient solution procedure is utilized for displacement solution; parallel computation is employed to improve the comput-ing efficiency and enable the modeling of a large lattice system.
Branch D: Update the analysis model and check the system stability. Through the conjugate gradient solution, nodal displacements are found and used to determine stresses in the lattice elements. Afterward, all critical elements (i.e., those reaching the compressive/tensile strength) are removed from the lattice mesh. The removal of these elements stands for localized damage during the printing process. The solution procedure within this time step repeats until no additional elements rupture and the system stabilizes. Eventually, the system state is regarded as a failure when reaching the predetermined system failure criterion. The structural failure criterion can be found in Section 2.6. Branch E: Output the system displacement and localized damage.
Branch F: Describe the system failure criterion as stated in Section 2.7.

Model establishment
In the lattice model, the continuum object is represented by a lattice mesh. When establishing a numerical model such as that of a cylinder structure, the cubical domain is first divided into a number of cubic cells, representative of the continuum system. Then, sub-cells are generated in the middle of each cell. Nodes are then randomly positioned within each sub-cell using a pseudorandom number generator (Qian, 2012). The ratio between the length of the subcell and the cell is defined as the mesh randomness, a representative for the mesh disorder. Considering that printable cementitious materials are not completely isotropic, the randomly distributed material properties are mimicked through irregularities of the network geometry (Berton et al., 2006;Bolander et al., 2005;Qian, 2012;H. Zhang et al. 2018b;H. Zhang, et al., 2018a). Based on our previous research, a randomness of 0.2 is adopted in this study to avoid large variations in lengths of individual lattice elements . As soon as the lattice nodes are generated in the cubical domain (as shown in Figure 3(a)), the target printed object is then established in this domain as indicated in Figure 3(b). After that, the nodes within the printed object are retained for model discretization while the outside nodes are deleted as shown in Figure 3(c).

Model discretization
After the model establishment, the target object can be defined using a set of lattice nodes. Employing these lattice nodes, the printed structure is discretized as a network of beam elements. An example of model discretization is shown in Figure 4(a). Delaunay triangulation is then performed to establish the node connectivity (Yip et al., 2005). A 3D printed object is manufactured by a layer-by-layer printing process. In this numerical model, the printed object comprises a series of layers, representative of the continuous printing process (as shown in Figure 4(b)). To account for the non-uniform gravitational load that occurs during the printing of a single layer, each layer is subsequently divided into a series of printing pieces (as shown in Figure 4(c)). Regarding the interface between different printing pieces, tie constraint is utilized, which means that different printing pieces share the same nodes and there is no relative motion between two neighboring pieces.

Element birth technique
The continuous printing process is simulated using the element birth technique, which is based on the deactiv- The continuous printing process with element birth iting and reactivating elements that compose the printed object. Each element has an initial time t 0 less than zero, and this initial time is determined by the element location in printed system. This initial time for the specific element can be derived based on where Δd represents the placement distance to the starting point of the printed structure; and v is the printing velocity. During the printing process, the element time label t increases by an incremental time step. Once the current printing time t is more than zero for a given printing segment, the elements representing that segment are activated. This process is illustrated in Figure 5 for the placement of segments 1 through 4.

Properties of lattice elements
Since the ratio between the length and cross-section of the lattice beam is generally low in the lattice network, Timoshenko beam elements are utilized to account for the shear deformation. To simulate the development of the mechanical properties in time due to hydration, the lattice beams are assigned time-dependent stiffness and strength properties. At present, two types of relations have been proposed in the literature to account for the development of mechanical properties of printable concretes at a very early age: A linear relationship Suiker, 2018;Wolfs et al., 2018 and an exponentially decaying relationship between the strength/stiffness and curing time (Panda, Lim, et al., 2019;Billberg, 2003;Perrot et al., 2016;Suiker, 2018;Wangler et al., 2016). The actual development of mechanical properties can be derived from experimental programs such as compression tests, direct shear tests (Wolfs et al., 2018) and triaxial compression tests . In principle, both types of relations can be incorporated into the lattice model to define the stiffness and strength development in time.
F I G U R E 6 Schematic view of (a) 3D printing object with sequence printing, (b) 3D printing model in lattice

Non-uniform gravitational load
In the 3D printing process, the material strength competes with the gradual increasing load, codetermining the buildability performance. In principle, little is known about the actual loads acting on the object during printing . The most significant load is certainly the self-weight, the magnitude of which can be established relatively easily. However, the effects of the kinematic pressure from the filament deposition and non-uniform gravitational load caused by the deposition process have not been quantitatively studied.
In general, the non-uniform gravitational load is ascribed to the continuous printing process during the printing period of a single layer. Since printing is a continuous process, the deadweight loading for the newly printed layer is gradually transferred into the whole system, resulting in the non-uniform loading force in the vertical direction. Consequently, this non-uniform gravitational load may threaten the overall stability of the object being printed, allowing fewer layers to be printed than otherwise expected. However, this influence has not been studied in any of the previous prediction models (Kruger Cho, et al., 2019;Roussel et al., 2020;Suiker et al., 2020;Wolfs et al., 2018. To account for the continuous printing process, each printing layer is divided into several sequential printing pieces (as indicated in Figure 6). In the numerical analysis, the value and position of system gravity are associated with the number of activated printing pieces. Thus, the influence of nonuniform gravitational load can be investigated through the number of printing pieces as shown in Section 3.2.
In the lattice model, the spatial variability is taken into account from two aspects, namely, model discretization and material variability. The model discretization, as described in Sections 2.2 and 2.3, generates randomly distributed lattice nodes that are connected by beam elements with different lengths and orientations. These parameters introduce spatial variability into the system. For 3DCP, the gravitational load is applied to the system as nodal forces. As mentioned before, the continuum system is discretized by a series of lattice nodes connected by the Timoshenko beams using the Delaunay triangulation. Given that the Voronoi tessellation is used to partition the domain, the volume of each Voronoi cell is used to assess the self-weight load for the corresponding nodal force (Pan et al., 2017) as shown in Figure 7. This volume can then be used to determine the equivalent gravitational force acting on each node in the mesh as = × × ( 2) where f represents the nodal force; V and ρ are the nodal volume and unit weight of the material, respectively; g is gravity acceleration (g ≈ 9.8 m/s 2 ).

Element failure criterion
As discussed before, the printing object is discretized by a series of lattice beams. Herein, the removal of lattice beams is utilized to account for localized damage that occurs during the additive manufacturing process. Again, in 3DCP, two competing parameters determine whether an object can be successfully printed or not. Specifically, the increasing strength and stiffness need to keep up with the gradually increasing load as more layers are deposited from the nozzle. The former can be incorporated into the lattice model through the development of mechanical properties in time, and the latter can be reflected via the incremental load in each step. Subjected to non-uniform gravitational load and boundary conditions, linear elastic calculations are performed to determine the stress distribution within the lattice mesh. The element failure criterion, consisting of normal force and bending moment determines the critical element and is given as Here, F is the normal force in the lattice beam element; M i and M j stand for the bending moments in the local coordinate system; and σ yield stands for the material yield stress. Specifically, a lattice element is assumed to break in compression or tension in accordance with the maximal stress theory, where the limiting stresses are represented by the uniaxial compressive and tensile strengths (Suiker, 2018;Suiker et al., 2020). A is the element cross-sectional area, and W is the section modulus (W = πD 3 /32, where D is the effective diameter of the lattice element). The coefficients α N and α M are the normal force factor and the bending influence factor, which control the degrees to which normal force and bending moment govern the failure mode. The value for α N is commonly adopted as 1.0 (Schlangen et al., 1997;D. -W. Zhang et al., 2016;H. Zhang et al., 2017), and the α M value is set to 0.05 in agreement with the literature Qian, 2012;Qian et al., 2017;Schlangen, 1993;Vassaux et al., 2014;D. -W. Zhang et al., 2016).

Structural failure criterion
In previous 3D printing models (Suiker et al., 2020;Wolfs et al., 2018), the structure is considered to have failed when any point of the printed object has reached the material yielding. This failure criterion may be reasonable for a homogeneous system under uniform loading. Under non-uniform loading, however, the point reaching the yield strength generally causes only localized damage instead of system failure. The printed structure may still be able to retain the shape and allow for the printing process to continue since some stress redistribution could occur. Therefore, this kind of failure criterion may underestimate the buildability performance under the complex printing geometry or the non-uniform gravitational loading.
To accurately model the printing process and quantify the buildability, a different structural failure criterion should be established. In the 3DCP experiments, the structures usually fail because of the system instability or material yielding. This means that a relatively large offset occurs between the new printing layer and the designed position. As a consequence, a new printing layer may fail to be placed into the system at a stable stage. Even if the object somehow does not collapse, such excessive deformation may be considered a "failure" when the aim is to print a predefined geometry, since the printing process is usually pre-programmed and does not correct for such deviations. In this numerical model, failure occurs when the lateral offset of a layer surpasses the width of an individual layer (i.e., when the next printing layer fails to find a position to be placed into the printed system) as F I G U R E 8 Structural failure for 3DCP due to large offset of the top layer F I G U R E 9 Plastic collapse dominant failure mode for 3D printing. The plotted grayscale indicates the displacement magnitude (unit: mm) indicated in Figure 8. This printing process stops since the structure has reached the maximum printing height.
Regarding a point in the structure reaching the limiting stress in tension or compression, the corresponding element will be removed from the printed object, representative of the localized damage during the printing process. Such local damage is typically accompanied by stress redistribution, which allows for the printing process to be continued. In other words, by modeling failure as an incremental process rather than as a single event, the critical printing height is less likely to be underestimated.

2.9
Failure mode Due to the absence of molds, the object may collapse during the printing process. A model must be able to capture the typical failure modes to quantify buildability performance. Herein, using the lattice model, plastic collapse results from the accumulation of local damage in the lower layers as depicted in Figure 9. The details related to the calibration of input parameters and experimental validation are given in the next section.

Model applicability
In this model, the solution procedure is based on the lattice beam model with the assumption of the linear elastic F I G U R E 1 0 Dimensions of the printing objects: (a) cylinder structure, (b) square structure. The rounding radius for the corners of the square structure is 50 mm (unit: mm) theory and small deformation. In practical applications of concrete 3D printing, the onset of elastic buckling or plastic collapse is the only thing that matters in a 3D printing experiment. Once the system fails, excessive deformation will occur. In that sense, large deformation is not relevant for predicting the critical printing height. Some structures, such as the wall geometry, are sensitive to elastic buckling. For these structures, the geometric non-linearity plays a significant role in the buildability quantification since even a small offset can introduce significant stress in the structure. Other geometries usually fail due to material yielding. This study focuses on the plastic collapse dominant failure mode with a relatively short printing time.
The second-order effects are not considered, and structural instability due to the loss of geometry (e.g., elastic buckling) is outside of the scope of the current study.
In addition to the mechanical analysis, in practice, the process of 3DCP is influenced by hydrothermal phenomena through a combination of cement hydration (and associated heat development and cooling), shrinkage (autogenous, plastic and drying) and viscoelastic deformation (creep). All of these factors may, to a certain extent, influence the buildability quantification, especially for cases with long printing time. The influence of most of these factors has not yet been experimentally studied and is poorly understood. The present model can therefore be seen as an initial processing step for the buildability performance quantification. Extensive investigations on the aforementioned factors are beyond the scope of this research and will be studied in the future.

NUMERICAL ANALYSES
The experimental validation of the lattice model for buildability quantification is performed by simulating two printing geometries (see Figure 10) using two types of materials, namely, (a) a cylindrical geometry with material model A and (b) a square structure with material model B.
The 3DCP numerical analyses are performed in three steps using the lattice model. First, the computational uniaxial compression tests are carried out for calibration of the time-dependent stiffness and strength. Then, two printing parameters including the non-uniform gravitational load and randomly distributed material properties are quantitively analyzed to study their impact on the printing characteristics, namely, the critical printing height and failure mode. On the basis of these parametric analyses and calibrated time-dependent material properties, two 3D printing numerical campaigns are carried out for model validation, in which the derived printing characteristics are compared with experimental results from the literature.

Computational uniaxial compression test
In this section, the unconfined uniaxial compression tests (UUCT) are simulated to calibrate the time-dependent stiffness and strength of lattice elements. As previously noted, two types of printable materials are considered. Material model A is regarded as the case study to demonstrate the detailed steps for parameters calibration.
During the 3D printing process, cementitious materials undergo a transition from a fluid-to a solid-like stage at an early age. The development of yield stress (static yield stress to be more precise) can be regarded as the transition between these two stages. When it comes to the buildability performance, the material is actually in the transition process, in which the cement hydration continues, impacting the physical and rheological behavior. To be specific, when the material stress reaches the yield stress, the hydrate bridges between percolated cement particles break and they flow. As a consequence, the printable materials display a viscous behavior and the shear rate is proportional to dynamic yield stress based on the plastic viscosity (Buswell et al., 2020;Roussel, 2018;Roussel et al., 2012Roussel et al., , 2020. Once the structure reaches the critical weight load, it will collapse. In such a "load-controlled" situation, further stress redistribution within the whole system is not possible as explained by Suiker (2018). It can be concluded that the yield stress plays a dominant role in the final printing characteristics. Therefore, in the initial stage of this model developed in this research, the local post-peak behavior is not considered: Only the material yield stress and elastic modulus are considered for the buildability performance quantification.
Regarding the determination of material properties at different curing times, five computational uniaxial compressive tests from the literature (Wolfs et al., 2018) are simulated to determine the properties of lattice elements to reflect the time-dependent material properties. In the experiment, the printable material is cast into steel F I G U R E 1 1 Computational compressive test. Loading is applied in the z-direction cylindrical molds and undergoes a compaction process three times on a 30 Hz vibration table for a homogeneous sample realization. Then, cylindrical samples with a diameter of 70 mm and a height of 140 mm are loaded through an INSTRON setup in which a Teflon sheet is placed on both sides of the specimen to minimize the friction at the supports. The corresponding mechanical properties for fresh concrete at t = 0, 15, 30, 60, and 90 min are tested.
Because the printable material subjected to uniaxial compression test undergoes a compaction process, a homogeneous sample with increased strength and stiffness and decreased compressibility is achieved as explained by Wolfs et al. (2018). For the computational uniaxial compression test, a numerical model of a cylinder structure with the same size as the experimental specimen is created. Considering the computational time and accuracy, a cell size of 5 mm is selected in the numerical campaign. The model in UUCT consists of 4317 nodes and 30,614 Timoshenko beam elements. Considering the low friction support condition in the experiments, radial deformation of the top and bottom boundary conditions are free. Nodal displacements on one side are applied and the opposite side is fixed as shown in Figure 11. In the lattice model, there are four input parameters that determine the material properties: Elastic modulus, shear modulus, compressive strength, and tensile strength. The compressive strength is assumed to be 10 times higher than tensile strength to consider the fact that the cementitious materials have a higher resistance to compressive loading than tensile loading (H. Zhang et al., 2018a;. Concerning the calibration process, Young's modulus can be computed dependent on the initial slope of the curve, and the compressive strength is calibrated based on the peak load. The detailed calibration process is based on the previously published research on hardened cementitious materials and may need to be reassessed in the future . After calibration, the time-dependent stiffness and strength used in the lattice model are listed in Table 2. Note that, in general, the local mechanical properties of the individual lattice element differ from the global ones (Delaplace et al., 2007). For example, the lattice element mechanical properties with the compressive strength equal to 10.96 kPa and elastic modulus equal to 77.9 kPa are calibrated by the computational uniaxial compression test in 0 min. The derived Young's modulus and compressive strength are compared to the experimental data (see Figure 12). In the literature, two types of "curing" functions have been proposed for describing the development of early age mechanical properties of 3DCP: Linear curing and exponentially decaying curing functions (Perrot et al., 2016;Suiker, 2018). Here, the linear function is adopted for the material properties calibration based on the previous research (Wolfs et al., 2018). The development of the time-dependent strength and stiffness within the first 90 min can be determined as indicated by the black solid lines shown in Figure 12 and the R2 values are calculated through the linear regression function and the average experimental results on different concrete stages.
The time-dependent stiffness and strength for lattice elements are described using the following linear relation:  (Wolfs et al., 2018) where E represents Young's modulus (kPa), which affects the elastic deformation; f c stands for the material compressive strength (kPa), and t is curing time within 90 min; the superscripts A and B stand for the material type.

Parametric analyses
As mentioned in Section 3.1, the cylinder specimen for UUCT underwent a compaction process (Wolfs et al., 2018). This vibration process reduces the local porosity, breaks the material internal structure, and promotes the re-flocculation and structuration process, thereby probably resulting in a (relatively) homogeneous sample with higher material strength, compared to the one in the 3D printing process. The material in the 3D printing process goes through a pump and is extruded through a hose and a nozzle, which may also result in higher local porosity . Consequently, the material subjected to a 3D printing process will probably have lower mechanical properties (i.e., strength, stiffness) and larger variations in randomly distributed material properties. Thus, considering the variability of material properties may therefore be crucial for realistic buildability quantification of 3DCP. Besides, as discussed in Sections 2.1 and 2.4, the continuous printing process in the experiment leads to the non-uniform gravitational load to the structure. This type of loading condition leads to localized damage, which may eventually result in structural failure. In the following sections, the influence of non-uniform gravitational load and randomly distributed material properties are investigated for the buildability quantification.

Non-uniform gravitational load
In this section, the impact of non-uniform gravitational load on the printing characteristics and geometry will be quantitatively analyzed. Two kinds of geometries under various non-uniform gravitational load scenarios (see Table 3) are analyzed. For the numerical analysis, each printing layer is divided into several printing pieces, which stands for the continuous printing process in the experiment and reflects the influence of non-uniform gravitational load as explained in Section 2.4. Here, seven numerical models with a different number of printing pieces and geometries are established to study their impact on the printing characteristics, and the corresponding input parameters can be found in Table 4. The model with a number equal to 1 means that one layer is added into the system in each step, while the model with a piece number equal to 10 will divide one layer into 10 sequential pieces to print. Note that the high friction at the bottom is utilized as a boundary condition. Figure 13 provides the critical printing height for two kinds of printing geometries. It can be illustrated from Figure 13(a) that the non-uniform gravitational load decreases the number of printing layers at failure for the cylinder structure. The critical printing height converges to a stable value by increasing the number of printing pieces. On the other hand, the non-uniform gravitational load plays a negligible influence on the square wall layout as shown in F I G U R E 1 3 Relation between the printing divisions per layer and printing characteristics: (a) cylinder geometry, (b) square wall layout Figure 13(b). Regarding the failure mechanism, the nonuniform gravitational load may reduce the critical printing height due to the uneven displacement for each layer. Specifically, the non-uniform gravitational load introduces uneven stress in each layer, leading to localized damage to the structure. On the other hand, the model with more printing pieces updates the system mechanical property in a finer time step, and as a consequence, stronger loading capacity for each layer is accompanied. These two opposing factors codetermine the critical printing height. Which one plays a more significant role depends on the printing geometry and the printing time. Based on this numerical analysis, it is observed that the non-uniform gravitational load influences the critical height for the model with relatively low printing speed or the numerical campaign with relatively large size, while the non-uniform loading hardly affects the printing critical height for the objects with fast printing speed and small printing size. Figure 14 shows the ultimate global failure related to the influence of non-uniform gravitational load, clearly showing that non-uniform gravitational load influences the final failure mode. For the model under uniform gravitational load, namely, models 1 and 5, the system deforms in an axial-symmetric way. In contrast, the failure zone in all models subjected to non-uniform loading arise from the non-uniform stress distribution to the structure, resulting in the uneven deformation per layer. Consequently, the radial deformation in the top area increases, such that the next printing layer may fail to be placed on the system within the printing accuracy.

3.2.2
Randomly distributed material properties As mentioned before, two methods are used to introduce spatial variability within the printed sytstem: Random placement of the lattice nodes (within each sub-cell) and random assignment of the material properties using TA B L E 4 Input parameters for the numerical analyses on non-uniform gravitational load Gaussian distributions. Due to the nature of the printing process, variability of the material properties might be high, yet this phenomenon remains scarcely investigated. Although investigations regarding the effects of vibration on early age concrete have been performed (Dunham et al., 2007;Hulshizer et al., 1984), they mainly focused on the 7 or 28 d material strength influenced by the vibration. Some studies do provide information about the strength measured 24 h after casting, but this is far beyond the timeframe of a concrete 3D printing process (C. Zhang et al., 2005). Herein, randomly distributed material properties are introduced into the 3D printing model. Specifically, the beam elements have compressive/tensile strength values that are randomly determined from a Gaussian distribution (with an average value of 10.96 kPa for the case of compressive strength at t = 0 s). This is an effective means for assessing the potential influences of material heterogeneity, in lieu of more advanced representations of material heterogeneity. Although the value of Young's modulus should correlate with strength, it is taken as constant herein at a given printing time for single variable comparison. Two standard deviation values, that is, 2 and 3 kPa, are chosen for qualitatively evaluating the effect of randomly distributed material properties on the buildability performance of 3D printing. During the printing process, the compressive strength and Young's modulus increase with the printing time following the growth rate given by Equation (4). For the material strength distribution, compressive strength is assumed to be 10 times greater than tensile strength for each element as discussed before. Note that, concerning the Gaussian distribution, a limited number of elements should be assigned with negative strength values, which is not physically possible. Therefore, those elements are removed from the mesh prior to the analysis.
Two numerical analyses of the hollow cylinder are performed with the same loading and boundary conditions as described in Section 3.2.2. Figure 15 shows the maximum radial displacement and the ultimate failure-deformation mode for these two examples before failure. It can be concluded that allowing for disorder in the lattice system, in terms of strength variation, reduces the number of layers that can be printed. To be specific, heterogeneous material behavior leads to uneven displacement for the same layer under self-weight. Some points of the target object reach their strength limits producing localized damage. The accumulation of localized damage leads to extensive deformation, progressive collapse and eventually system failure. Increasing the standard deviation of strength reduces the achievable printing height. Introducing strength variations also produces a more asymmetric failure mode relative to that of the homogeneous model. It can be stated that, in practice, the printing process should aim to ensure a homogeneous material deposition (e.g., in terms of density and porosity) to maximize buildability. F I G U R E 1 6 Failure-deformation mode during the printing process (unit: mm)

Model validation
In this section, two 3D printing experimental campaigns, namely, a cylinder and a square structure, are modeled employing the lattice model for buildability quantification as indicated in Figure 10.

Cylinder geometry
A cylinder structure with a diameter of 500 mm, a thickness of 40 mm, and a layer height of 10 mm is modeled to quantify buildability performance. The details about cylinder geometry can be found in Figure 10 and Table 3. The numerical model is in line with the specimen in the trial printing (Wolfs et al., 2018) as indicated in Figure 16.
In the numerical program, each layer is subdivided into 60 pieces for sequential printing reflecting the non-uniform gravitational load during the printing process. The interface between the different pieces is modeled using tie constraint as described in Section 2.2. The computational printing process advances until structural failure. More specifically, as soon as the linear elastic calculation fails to converge, the next printing piece in the numerical analysis cannot be placed due to the excessive deformation as described in Section 2.7. The time interval between two layers is set to 0.31 min, determined by the printing speed in the experimental campaign. The development of material stiffness and strength depends on the age of the printing piece. For example, after the placement of 5 layers, the initial printing piece in the first layer is 1.55 min old, and the corresponding mechanical properties are assigned to this piece, while the current printing piece is assigned the initial mechanical properties (i.e., those corresponding to t = 0 min). In the UUCT (Wolfs et al., 2018), the printable materials show a standard deviation value equal to 1.07 kPa for UUCT. Based on this global standard deviation, a standard deviation of 1.96 kPa for the lattice elements can be derived and this value is adopted in this numerical analysis.
In the 3D printing model, the cell size is set to 5 mm, the same as the computational uniaxial compression test. The same cell size can eliminate the influence of element size on the material properties. Each layer consists of 5000 nodes connected by around 360,000 beam elements. Concerning the boundary condition, the bottom layer is fixed to model the high frictional resistance of the printing bed, corresponding to the conditions of the experiment (Wolfs et al., 2018). The relevant gravitational load with the average material density of 2070 kg/m 3 is calculated by employing Equation (2). Figure 16 shows the maximum structural deformation affiliated with the printing process and ultimate global failure. Specifically, the final failure-deformation is affected by the lateral restraint in the radial direction and increasing non-uniform gravitational load. The stepwise addition of printing pieces causes eccentric loading, introducing bending moment to the system. The element stresses increase under loading, such that the randomly assigned strength values are exceeded at some locations, causing localized damage. This localized damage negatively affects the structural stability and significantly governs the critical printing height. Eventually, the system fails at the 41st printing layer, with the radial deformation equal to 14.42 mm, which occurs at the corresponding z position of 115.36 mm. Figure 17 shows the relation between the radial deformation and corresponding position in the z-direction, clearly showing a gradual increase of radial deformation as more layers are printed. The model with a total of 41 filament layers is successfully printed with a total building height of around 410 mm within 12.71 min. Concerning the printing process, excessive deformation occurs when placing the 41st layer, and the system eventually collapses after the placement of the 42nd layer. The excessive radial deformation of the top area at the layer 41 stage is noticeable, demonstrating that the current printing layer cannot be placed in the predetermined position. Thus, the object fails to ensure the printing quality and meet printing accuracy due to the excessive deformation in the top zone; therefore, TA B L E 5 The overview of printing characterizations between the current lattice model and previous simulations and experiments the 41st layer is regarded as the critical printing height. The maximum radial deformation and the corresponding z position are derived from the stage before the system collapse (i.e., layer 41). The final failure mode is dominated by the accumulation of damage leading to plastic collapse. Table 5 gives a comparison of the final printing characterizations between the newly developed lattice model, and the FEM and experimental result from the literature (Wolfs et al., 2018). The average number of 41 layers derived from the lattice model overestimates the experimental findings for critical printing height (29 layers) by 41.38%. The numerical average radial deformation before failure is equal to 14.42, which yields a relative difference below 1% with the experimental values. The corresponding z position is 115.36 mm, which differs from the experimental result by less than 6%.
There is good general agreement between the results presented herein and those of the FEM employed by Wolfs et al. (2018). In particular, both sets of results overestimate the critical printing height but provide accurate estimates of maximum radial displacement and the height at which it occurs. This agreement is remarkable considering the large fundamental differences between the modeling approaches.
The fact that both differing modeling approaches overestimate the critical printing height suggests that there are aspects of the printing process that are underappreciated. The remaining discrepancy between the lattice model results and experimental data might be due to disregarding the influence of manufacturing imperfections in the printing process. The failure of cylindrical structures is notoriously sensitive to imperfections in geometry (Amazigo et al., 1972;Koiter, 1967;Leipholz et al., 1975). Similarly, the initial imperfections generated in 3DCP may decrease the critical printing height due to similar forms of structural instability. Furthermore, experimental techniques are needed for determining the in situ properties of the printed materials. The use herein of UUCT data, F I G U R E 1 8 Final failure mode for square wall layout (unit: mm): (a) final failure mode, (b) localized damage obtained from specimens that were consolidated through vibration, may be another source of error.

3.3.2
Square structure In this section, a relatively small square wall structure (in Figure 10(b)), with the layer width equal to 250 mm and a corner radius of 50 mm, is simulated for model validation and buildability quantification. Additional information about this model can be found in Figure 10(b) and Table 3.
In the numerical campaign, the bottom layer is fixed to correspond with the high friction on the print bed (Suiker et al., 2020). The time to place one layer is set to 9.6 s, determined by the printing speed in the experimental campaign. Each layer consists of around 4300 nodes, connected by around 31,000 lattice beams; each layer is divided into four printing pieces, which enables the representation of the non-uniform gravitational load. The time-dependent material stiffness and strength are described using Equation (5): Material model B. The nodal force is calculated based on the average material density of 2100 kg/m 3 employing Equation (2). Figure 18 shows the failure-deformation mode and the localized damage in plane A-A. A comparison with experimental data is given in Table 6. As for the analyses of the cylindral structure, under the gradually increasing gravity load, the lattice elements eventually exceed the material strength, causing localized damage in the printed object. This kind of damage negatively affects the structural stability and governs the critical number of printing layers. Eventually, the system fails at the printing of the 24th layer because of the accumulation of damage in the bottom layer. The lattice model is verified via comparison with the experimental results of Suiker et al. (2020). With respect to critical printing height, there is a relative difference of around 10% (Table 6), which is smaller than the relative difference observed for the cylinder structure. The square layout might be less sensitive to imperfections generated by the printing process, but this possibility needs further study.

CONCLUSION
In this study, the lattice model is adopted to quantify the buildability for 3DCP. This model, incorporating an element birth technique, the time-dependent nature of the mechanical properties, printing velocity, non-uniform gravitational load, localized damage and spatial variation of the printed object, is able to reproduce the plastic collapse failure modes reported in the literature. Using this model, computational uniaxial compression tests are first conducted to calibrate the time-dependent material stiffness and strength in the range of 0 to 90 min. Thereafter, parametric analyses regarding the non-uniform gravitational load and randomly distributed material properties are conducted to evaluate their impact on printing characteristics, including the failure-deformation response and critical printing height. Based on the calibrated material properties and parametric analyses results, the model is finally validated by comparison with two well-documented 3D printing experiments from the literature. The main conclusions are summarized below.

A new failure criterion for buildability quantification is
proposed and applied in the numerical model. Failure of the system is assumed to occur when the next printing layer fails to be placed on the current printed system. This failure criterion allows for the occurrence and accumulation of localized damage. Relative to a single event, stress-based failure criteria, the proposed criterion provides a more realistic measure of critical printing height. 2. In the 3DCP numerical analysis, the spatial variability of the printed object is approximated in two ways. First, the random positioning of nodes within the sub-cells produces an irregular lattice. This leads to non-uniform stresses in the elements even under uniform loading. Second, the strength properties of the lattice elements are randomly assigned, according to a Gaussian distribution, to mimic the material heterogeneity. A higher standard deviation value causes a lower critical printing height based on the parametric analysis results. 3. By means of the element birth technique, the continuous printing process can be simulated more realistically. The quantitative influence of non-uniform gravitation load can be reflected through the number of segments for each layer. The non-uniform gravitational load influenced the critical printing height of the cylindrical model, which was relatively large in size and constructed with relatively low printing velocity. The effect of non-uniform loading less affected the critical printing height of the square structure.
Lattice modelings of the square structure and cylinder structure produce the correct failure-deformation patterns and quantitative agreement with experimental data with respect to critical printing height and maximum lateral displacement of the printed structure. Furthermore, the lattice model results are comparable to those of the FEM employed by Wolfs et al. The tendency for both models to overestimate the critical printing height suggests there are aspects of the printing material behavior, and of the printing process, that need further study. In particular, the consequences of imperfections generated by the printing process could be investigated through further numerical studies.
Through this research work, a new numerical method is proposed to quantify the buildability for 3DCP. The present lattice model can be seen as an intermediate step toward accurate buildability quantification. We anticipate extending the lattice model to study the effects of geometric non-linearity, early-stage drying shrinkage and creep on 3DCP buildability.