Rheological Property Criteria for Buildable 3D Printing Concrete

Fresh concrete used in 3D printing should ensure adequate yield stress, otherwise the printed concrete layer may suffer intolerable deformation or collapse during the printing process. In response to this issue, an analytical study was carried out to derive the initial yield stress and hardening coefficient of fresh concrete suitable for 3D printing. The maximum shear stress distribution of fresh concrete was calculated using a stress transformation equation derived from the equilibrium condition of forces. In addition, the elapsed time experienced by fresh concrete during the printing processes was estimated and was then substituted into the elapsed time-yield stress function to calculate the yield stress distribution. Based on these results, an algorithm capable of deriving both the initial yield stress and the hardening coefficient required for printing fresh concrete up to the target height was proposed and computational fluid dynamics (CFD) analyses were performed to verify the accuracy of the proposed model.


Introduction
Reinforced concrete (RC) structures have been widely utilized because they are extremely easy and economical to build. However, the disadvantages of the construction method include the facts that it causes an inefficient consumption of materials, problems with the labour force and issues with the construction period in the process of manufacturing the formwork, processing/placing the reinforcement and pouring the concrete [1]. Furthermore, there has been a growing demand for freeform structures that require complicated and heavy formwork [2]. 3D concrete printing technology can be an effective alternative, which can remove the formwork and reduce cost, labour force and the construction period [3]. Therefore, numerous studies have been conducted on the concrete 3D printing technology [4][5][6][7][8][9][10][11][12].
Khoshnevis [4] developed a contour crafting (CC) that uses computer control to exploit the excellent surface-forming capability of trowelling and it is known as the first example of 3D concrete printing. Dini developed a D-shape process that produces a concrete element by spraying a binder onto a layer of powder in a manner similar to powder-bed printing or binder jetting [5]. This printing method is advantageous in that it can produce complicated elements more easily than fused deposition modelling (FDM). Le et al. [6,7] established the key concept of fresh concrete properties, such as buildability, extrudability, open time and the workability of fresh concrete, suitable for concrete 3D printing. In addition, they closely examined the hardened properties, such as the void and dry shrinkage, as well as the mechanical properties of printing concrete. Gosselin et al. [8] used the 3D printing technology to print out a complicated 3D form of concrete members without the use of formworks. Using this particular method, they also developed concrete members with improved performance in heat insulation and sound insulation. Hambach and Volkmer [9] implemented a 3D printing technology able to artificially control fibre orientation in fibre-reinforced Portland cement paste (hereafter referred to as FRC). Their research results confirmed that the flexural strength of the 3D printed FRC varied depending on the printing path and fibre type and it exhibited the highest flexural strength when the volume fraction of carbon fibre was 1% and the 3D printing process progressed along the optimized printing path [9]. Bos et al. [10] summarized the history of the study on additive manufacturing of concrete and addressed the characteristics of 3D concrete printing (3DCP) product geometry and structure. Buswell et al. [11] outlined the technical issues on 3DCP and presented a guide and vision for future research. They also described the effects of the fresh and hardened paste, mortar and concrete properties on geometry of the created object in detail. Labonnote et al. [12] performed mapping study by classifying additive construction into four main categories; material science, engineering, architecture and market analysis.
The constitutive equations of fresh concrete are commonly represented by the Bingham model or the Herschel-Bulkley model [13]. In the 3D printing process, fresh concrete may suffer deformation and collapse when the shear stress generated in concrete exceeds the yield stress [14]. Therefore, in previous studies [15][16][17][18], the maximum shear stress of fresh concrete was calculated based on the equilibrium of forces and compared with the yield stress in order to predict the slump of fresh concrete. Unlike in the slump test, in the 3D concrete printing process, the increase in the yield stress of the fresh concrete should be taken into consideration. The increase in the yield stress according to the elapsed time of fresh concrete is represented by a linear or exponential function [19][20][21]. Perrot et al. [22] predicted the failure of fresh concrete during the printing process by reflecting the increase of the yield stress of fresh concrete over time and they proposed a model for calculating the optimal build-up velocity for concrete printing. They also verified the proposed model through a failure test on fresh concrete. However, in a test conducted by Perrot et al. [22], the fresh concrete showed many deformations prior to its failure and the test results showed that, since these deformations can cause the outputs of other modified forms different from those of the input 3D modelling data, the mixture design of printing concrete should be done to avoid deformation and collapse, rather than failure. Meanwhile, Wangler et al. [23] suggested the printing speed required to prevent collapse or deformation of fresh concrete during layering process. Suiker [24] suggested a mechanical model that could consider buckling collapse and Roussel [25] also proposed a model that reflects the stability failure. Wolfs et al. [26] conducted research on early age mechanical behaviour of 3D printed concrete through numerical simulations and experiments. From their research findings, it was confirmed that the numerical simulations well predict the early age behaviour of 3D printed concrete. Although many researches have been conducted to investigate the behaviour of fresh concrete and to prevent the collapse of layered concrete, it is still necessary to develop analytical models that can evaluate the behaviour of 3D printed concrete in a simple or efficient way without complex numerical simulations.
The purpose of this study is to quickly obtain the initial yield stress and hardening coefficient required to prevent the deformation and collapse of fresh concrete during printing. To this end, the maximum shear stress distribution generated in fresh concrete during the built-up process was calculated using the force equilibrium condition based on the slump prediction model of Murata [15] and the concrete yield stress distribution considering the overtime was calculated using the elapsed time-yield stress equation. In the proposed model, it was assumed that the printed concrete collapses when the maximum shear stress exceeds the yield stress. In other words, the collapse of printed concrete was assumed to be dominated by material failure, excluding stability failure. In addition, ANSYS CFX (ANSYS Inc., Canonsburg, PA, USA), a commercial software program for numerical analysis, was used to verify the reasonability of the proposed model by checking the deformation of layered printed concrete. Figure 1 shows that even in printing work carried out for the same plan, different printing paths are possible depending on the printing start point of each layer (hereafter referred to as the printing start point) can have and the printing start point has a significant effect on the printing time gap. According to the test results conducted by Le et al. [7], as the printing time gap increases, the interfacial bond performance of the hardened 3D printed concrete decreases. This phenomenon was also reported by Feng et al. [27] and Nerella et al. [28]. Therefore, the printing start point needs to be carefully considered before the printing process begins, since it affects the mechanical performance of the hardened 3D printed concrete. This printing start point can be set in the slicing software and is classified into optimized printing start points for the fastest printing velocity (hereafter referred to as OP), in which the closest point after the output of one layer is designated as the printing start point of the next layer in order to minimize the printing time as shown in Figure 1a and the selected printing start points closest to specific location (SP), where the output of one layer, which is returned to the original point after its output, progresses again with the printing start point fixed to a specific location, as shown in Figure 1b. As Figure 1a shows, for the OP, the filaments of each layer exhibit different printing time gaps in the length direction of the printing path. On the other hand, in terms of the SP, the filaments of each layer show the same printing time gap at any point in the length direction of the printing path, as shown in Figure 1b. In this regard, this study assumed the case of SP, which can have the same printing time gap and ensure uniform interfacial bond performance and calculated the yield stress distribution and elapsed time in the printing of the concrete.   The yield stress of fresh concrete increases with the increasing elapsed time and thus the time at which fresh concrete is supplied affects buildability. Figure 2a shows an ideal continuous supply method, in which a continuous mixer supplies a fixed amount of fresh concrete per second required for printing. In this case, fresh concrete remains in the hopper for a very short time. On the other hand, the example shown in Figure 2b is of a discontinuous supply method in which fresh concrete has a certain amount of elapsed time before printing as it is supplied only when the hopper is completely emptied. In this case, fresh concrete is kept in the hopper for a relatively long time compared to a continuous supply method. A typical concrete mixer adopts the same method as that of the pan mixer shown in Figure 2b. Therefore, this study sought to derive the initial yield stress and hardening coefficient required for the build-up layers of fresh concrete under the assumption of the discontinuous supply shown in Figure 2b. In addition, the following assumptions were made to simplify the calculation of the maximum shear stress in fresh concrete.

1.
The accelerator is injected from the nozzle, as shown in Figure 3 and a fixed quantity is injected in proportion to the unit discharge amount.

2.
The cross section of fresh concrete is a perfect rectangle, as shown in Figure 4.

3.
Fresh concrete is built up vertically and the width (w), height (h) and length (l) of each layer of printed concrete are the same.

Maximum Shear Stress Distribution (Demand Curve)
In Figure 4, n represents the total number of built up layers and i means the order of the layers built up from the bottom surface. Normal compressive stress (σ) generated from any height z to an xy plane can be calculated as shown below [15][16][17][18].
where ρ is the density of fresh concrete and g is the gravity acceleration. The maximum shear stress (τ max ) that occurs in the layered fresh concrete can be calculated from the equilibrium condition of forces, as shown below [15][16][17][18]. Figure 4a shows the maximum shear stress distribution according to Equation (2). On the other hand, Figure 4b shows the maximum shear stress distribution when fresh concrete is not fully formed but partially formed within the same layer. In that case, the normal compressive stress and maximum shear stress generated from the position (y , z ) of a build-up layer point to the xy plane can be calculated as shown below.
In addition to the self-weight of fresh concrete, considering the extra shear stress (τ ex ) caused by discharge pressure at the point where the printing is carried out, if any, the maximum shear stress (τ max ) in the layered fresh concrete can be re-expressed, as follows.

Elapsed Time
The yield stress distribution of fresh concrete varies from hour to hour, even during the printing process, because the yield stress of fresh concrete changes with the elapsed time. The elapsed time can be calculated by dividing the length of path by the velocity of nozzle movement.  As shown in Figure 5a, if fresh concrete is vertically built up with the same width, length and height in each layer, the length of the printing path (l pr,i ) from the i th layer to arbitrary time (t) is classified into the length of the extrusion path (l ex,i ), the length of the horizontal returning path (l rh,i ) and the length of vertical returning path (l rv,i ), as shown in Figure 5b and this can be represented by the equation below.
Therefore, the total length (l ex ,l rh ,l rv ) considering the total printing process and the length of the total printing path (l pr ) can be calculated by summing up the printing paths at each layer, as shown below.
If the velocity of nozzle movement along each path is v ex ,v rh , v rv , the time(t ex,i ,t rh,i t rv,i ) required to move along each path at the i th layer is shown below.
Therefore, the time (t pr,i ) required for printing the i th layer is the sum of t ex,i , t rh,i , t rv,i and this can be represented by the equation below.
In the same way, the total time taken to move along each path (t ex , t rh , t rv ) and the total time spent on the total printing process (t pr ) can be calculated as follows.
On the other hand, since fresh concrete is formed only in the extrusion path during the printing process, the relationship between l ex and t pr is represented by a discontinuous function, as shown in Figure 5c and this can be represented using the floor function, as follows. Figure 6 shows the volume of fresh concrete supplied to the hopper and pipe. In this study, it was assumed that the amount of fresh concrete as much as the sum of the volume of the empty space (V e ) from the start point of the pipe to the end point of the nozzle and the volume of the hopper (V h ) is first supplied and then fresh concrete is supplied as much as V h . The length of the filament (l s ) that can be printed by fresh concrete as much as V h can be calculated as shown below.
Where w and h are the width and height of the filament, respectively. If the number of times for concrete supply is m, fresh concrete has to be supplied at a point in time l ex = (m − 1)l s . In this case, the mth fresh concrete supply point t s,m can be calculated by substituting (m − 1)l s into Equation (11), as shown below. Figure 7 shows the extrusion phase of fresh concrete. At the beginning of printing, the 1st-supplied fresh concrete contained in the pipe is extruded as shown in Figure 7a. After that, the 1st-supplied fresh concrete and the 2nd-supplied fresh concrete begin to be mixed in the pipe as shown in Figure 7b and then the mixed fresh concrete is extruded as shown in Figure 7c. After this step, the 2nd-supplied fresh concrete is extruded as shown in Figure 7d. Consequently, it is very difficult to predict the yield stress of the mixed fresh concrete in this multiphase flow process. Therefore, to calculate the yield stress on safe side and make it simple, this study set the elapsed time of fresh concrete to be zero when the 2nd fresh concrete is supplied into the hopper, although the remaining fresh concrete in the pipe is the 1st-supplied fresh concrete or mixed fresh concrete. The elapsed time (t e ) experienced by fresh concrete is the time elapsed from the mixture of concrete to the moment of printing. Thus, it can be calculated by the difference between t pr and t s,m , as shown below.
In the study conducted by Le et al. [6,7], an accelerator was added to rapidly increase the yield stress of the fresh concrete. Since the accelerator is injected in the vicinity of the nozzle, the velocity of concrete hardening increases from the time when concrete is extruded. Therefore, the relationship between elapsed time and the yield stress of fresh concrete changes from the point in time of the extrusion and the total elapsed time (t e ) experienced by fresh concrete can be classified into the elapsed time experienced before extrusion(t b ) and the elapsed time (t a ) experienced after extrusion.
If the length of the extrusion path for concrete present at an arbitrary location to be printed is l o and the time when concrete located at l o is extruded is t o , the t o can be calculated by substituting the l o into t pr (l ex ) function, as shown below.
When the fresh concrete located at l o is extruded, the number of times concrete must be supplied can be calculated using the floor function, as follows.
The t s,m of concrete located at l o can be estimated by substituting m calculated from Equation (17) into Equation (13). Since t b is the elapsed time experienced from the mixture of concrete to the extrusion of concrete, it is the difference between t o and t s,m . Likewise, as t a is the elapsed time experienced from the extrusion of concrete to the moment of printing, it is the difference between t pr and t o . Therefore, t b and t a can be represented by the following equations, respectively. Figure 8 shows the elapsed time t b , t a calculation procedure of concrete located at points A, B and C in the printing completion stage during the printing process of building up four layers. For example, since concrete at point A is located in l o,A , the time when printed is t o,A and concrete located at A is supplied from t s,1 , the equation becomes t b,A = t o,A − t s,1 . In addition, since the printing end point is

Relationship between Elapsed Time and Yield Stress
To define failure criteria, Wolfs et al. [26] used Mohr-Coulomb criterion, while Murata [15], Christensen [16], Pashias et al. [17] and Saak et al. [18] used Tresca criterion, in which the yield stress of fresh concrete was obtained from rheometer measurements. They reported that both Tresca and Mohr-Coulomb criteria well predicted the failure of fresh concrete and thus this study adopted Tresca criterion that is simpler than Mohr-Coulomb criterion.
Previous research [19][20][21] found that a liner or exponential function can be used to approximate the increase in the yield stress of fresh concrete according to the elapsed time. That is, the elapsed time-yield stress equations can be represented as shown below. τ y (t e ) = τ 0 e αt e (19a) τ y (t e ) = αt e + τ 0 (19b) Where α and τ 0 are the hardening coefficient and the initial yield stress when the accelerator is not used. Equation (19a) or Equation (19b) can be selected and applied based on the yield stress measurement results of fresh concrete. However, the yield stress derived from the elapsed time-yield stress equation should be always smaller than the measured value in order to obtain safe results. It is noted that yield stress of fresh concrete should be measured after discharge because the yield stress can be changed due to the energy introduced into the concrete during the pumping and transport processes.
If the accelerator is injected into fresh concrete, the yield stress of the fresh concrete increases more rapidly from the elapsed time t b of the accelerator injection. To take this phenomenon into consideration, a hardening coefficient β was additionally introduced in this study, where β is a constant that can be obtained by measuring the yield stress of fresh concrete according to the elapsed time after the accelerator is injected into the fresh concrete. The elapsed time-yield stress equation can be represented as shown below.

Yield Stress Distribution (Capacity Curve)
If the accelerator is not injected into the fresh concrete, the yield stress according to the position of the fresh concrete (l o ) and the position where the fresh concrete is extruded (l ex ) can be calculated using Equations (11), (13), (14), (17) The yield stress calculated using Equation (22a,b) are a function with respect to l o , l ex and the yield stresses of Equation (22a,b) should be represented as a function of (y,z) for comparison with the maximum shear stress (τ max ) represented by Equation (5). In addition, since the maximum shear stress is largest at the bottom of each layer, the yield stresses for each layer should be compared at that same location. The following transformation equations can be used to represent l o and l ex as a function of (y,z) and a function of (y ,z ) at the bottom of each layer.
That is, the yield stress τ y (y, z, y , z ) at all positions on the y − z coordinate plane can be calculated using Equation (24a-d).

Demand τ o , α and β
Deformation occurs when the maximum shear stress generated in layered fresh concrete exceeds the yield stress [14][15][16][17][18]. To prevent this, the yield stress of the fresh concrete should be greater than the maximum shear stress to satisfy the following relation.
Therefore, the yield stress (τ 0 ) required to build up one layer of fresh concrete at the beginning of printing without deformations should be satisfied in accordance with the following equation.
If the accelerator is not injected into fresh concrete, the α required to prevent deformation that may occur in the build-up process of fresh concrete can be calculated using Equations (5) and (24a), (24b) and (25).
If the accelerator is injected into fresh concrete, the β required to prevent deformation that may occur in the build-up process of fresh concrete can be derived using Equations (5) and (24c), (24d) and (26).
Based on Equations (27a) and (27b), the α(y, z, y , z ) required to prevent collapse and deformation can be calculated by changing (y , z ) from the start point of extrusion to the coordinates (y f , z f ) of the end point with respect to all the (y, z) within the area in which fresh concrete is layered. In addition, the α(y, z, y , z ) should have a larger than maximum value (α max ) among all the α values calculated in the printing process based on Equation (27a,b), since the mixture proportions of fresh concrete used in printing are all the same. If the accelerator is injected into fresh concrete, the α should be assumed to be as shown in Equation (28a,b) in order to calculate the β(y, z, y , z ). Even if fresh concrete with the same mixture proportion is used in the build-up process, the β(y, z, y , z ) may vary depending on the accelerator injection time (t b ). Therefore, the maximum value (β max (y, z)) among the β(y, z, y , z ) values calculated from the arbitrary (y, z) is calculated based on the β(y, z, y , z ) values calculated using Equation (28a,b), and represented on the t b − β max plane and the mixture design should be done so that β of fresh concrete can always have β max or more. Figures 9 and 10 show algorithms for calculating τ 0 and α max when the accelerator is not injected into the fresh concrete and τ 0 and β max when the accelerator is injected, respectively.   Table 1 shows the summary of variables for two example cases that are selected for the validation purpose of the proposed model. The size of the 3D printed structure and the width and height of each layer were selected in the consideration of the time required for the CFD analysis. The main variable is the hardening coefficient (β) that changes with the accelerator input and the other properties in Cases 1 and 2 were assumed to be the same. From Equation (26), the τ 0 value, which is the minimum yield stress to prevent collapse of a layer while the fresh concrete of upper layer is being discharged, is calculated as 122.5 Pa. According to Petit [19], the hardening coefficient (α) of self-compacting concrete and mortar depends on the proportions of the components of the mixture or temperature of the fresh concrete, ranging from 0.00013 to 0.00025/s. In this example, the α value was assumed to be 0.0003/s based on the fact that the hardening rate of printable concrete is higher than that of the self-compacting concrete. In addition, τ 0 was set to 500 Pa that is higher than 122.5 Pa, which was determined to induce the occurrence of deformation immediately after a predetermined number of layers are built up.  Figure 11 shows the comparison between the maximum shear stress distribution and the yield stress distribution analysed using variables corresponding to Case 1 of Table 1, during which the coordinates of the printing end point are (y f , z f ). In Figure 11, the red line indicates the maximum shear stress (i.e., shear demand) due to layering of the extruded concrete and the blue line indicates the yield stress (i.e., shear capacity). The demand curve and the capacity curve do not meet each other until the time at which the fifth layer is built up so that deformation does not occur in the layered concrete. When the sixth layer is built up, the deformation of the concrete in the bottom surface occurs when the maximum shear stress exceeds the yield stress.    Figure 12 shows the calculation result of β max , which does not cause the collapse of fresh concrete during the build-up process under the same printing conditions as found in Case 1, according to the elapsed time before extrusion (t b ) with the use of the algorithm shown in Figure 10. The analytical results obtained through the proposed model show that the collapse of the fresh concrete does not occur when the β max is greater than 0.0069/s and the value β was set to 0.0069/s in Case 2. In the following section, CFD analysis was performed for Cases 1 and 2 and the reasonability was verified through a comparison between the proposed model and analysis results.

Validation Using Computational Fluid Dynamic Analysis
Computational fluid dynamics (CFD) is one of the analytical techniques used to predict the behaviour of fluids. It has been widely used to simulate the behaviour of fresh concrete analysis [29]. In this study, ANSYS CFX, a commercial software program for numerical analysis, was used to perform a CFD analysis in order to verify the reasonability of the proposed algorithm. The analysis parameters used in the verification are Cases 1 and 2, shown in Table 1. Figure 13 shows the modelling for the CFD analysis. As shown in Figure 13a, the layered printed concrete consists of the air domain and the fresh concrete domain, while a rectangular mesh composed of 41,924 nodes and 30,698 elements was used as shown in Figure 13b. As for the boundary conditions, an opening was applied to the upper part of the model and a wall (i.e., no slip condition) to the other sides. The modelling also ensured that the force of gravity acts on the whole domain and the Bingham model was used in the CFD analysis. The yield stress value obtained when the analysis was performed under the conditions of Cases 1 and 2, shown in Table 1, was applied with respect to the yield stress (τ 0 ).   Table 1. The left graphs of Figure 14a,b show the yield stress distribution and the maximum shear stress distribution in the case of y = y f = 0.075 m in the proposed model. It should be noted that in Cases 1 and 2, the other properties are the same except for the hardening coefficient (β), due to the injection of accelerator. The β of Case 2 is the hardening coefficient, which is calculated through the algorithm proposed in this study required to prevent the collapse of fresh concrete. In Case 1, the proposed model well predicted the CFD simulation results, in which the fresh concrete collapsed at a height of approximately 0.018 m. On the other hand, in Case 2, the fresh concrete did not collapse in both results by the proposed model and the CFD analysis.

Conclusion
This study aimed to develop an algorithm to quickly derive the rheological properties needed to prevent the collapse of fresh concrete during the printing process. In the proposed model, the fresh concrete was regarded as Herschel-Bulkley fluid and it was thus assumed that no deformation occurs before the maximum shear stress of layered concrete exceeds its yield stress. The maximum shear stress in the layered fresh concrete was calculated from the force equilibrium condition and the yield stress of the fresh concrete at any location during the printing process was estimated considering the elapsed time that depends on the printing path. Based on this research, the following conclusions were drawn: 1.
The proposed model could reflect the increase in the yield stress of the fresh concrete during the 3D concrete printing process by substituting the elapsed time into the yield stress-elapsed time equation.

2.
Based on the CFD analysis results on the cases examined in this study, it was confirmed that the proposed model can provide very accurate estimations on the occurrence and location of collapse. 3.
The proposed model would be very helpful to obtain the rheological properties of fresh concrete required for 3D printing, such as hardening coefficient and initial yield stress, without any complex numerical simulations.

4.
Further researches are still required to check the stability problems of the layered concrete for vertically complicated structures.