3D 콘크리트 프린팅에서 철근 통합에 대한 전산 유체 역학 모델링 및 실험적 분석

Md Tusher Mollah, Raphaël Comminal, Wilson Ricardo Leal da Silva, Berin Šeta, Jon Spangenberg

Abstract

A challenge for 3D Concrete Printing is to incorporate reinforcement bars without compromising the concrete-rebar bonding. In this paper, a Computational Fluid Dynamics (CFD) model is used to analyze the deposition of concrete around pre-installed rebars. The concrete is modelled with a yield-stress dependent elasto-viscoplastic constitutive model. The simulated cross-sections of the deposited layers are compared with experiments under different configurations and rebar sizes, and found capable of capturing the air void formation with high accuracy. This proves model robustness and provides a tool for running digital experiments prior to full-scale tests. Additionally, the model is employed to conduct a parametric study under three different rebar-configurations: i) no-rebar; ii) horizontal rebar; and iii) cross-shaped (horizontal and vertical) rebars. The results illustrate that air voids can be eliminated in all investigated cases by changing the toolpath, process parameters, and rebar joint geometry, which emphasizes the great potential of the digital model.

Keywords

3D Concrete Printing (3DCP); Reinforcement bars (rebars); Computational Fluid Dynamics (CFD); Multilayer deposition; Air voids

1. Introduction

3D Concrete Printing (3DCP) [1] is an extrusion-based automated construction process that belongs to Digital Fabrication with Concrete (DFC) [2,3]. The 3DCP offers high-quality built-structures with customizable structural design in a cost- and time-efficient manner [[4], [5], [6], [7]]. Structures in 3DCP are fabricated in a layer-by-layer approach, where a concrete extrusion nozzle is controlled by a robotic arm, cylindrical robot, gantry system, or delta system [[8], [9], [10], [11]]. Despite the enormous potential of 3DCP, one of its crucial limitations is the integration of reinforcement for the production of load-bearing structures.
Most structural applications require the use of reinforcement to withstand tensile forces and introduce structural ductility [[12], [13], [14], [15]]. However, the introduction of reinforcement with 3DCP has never been an easy task, and difficulties were recognized at early stages of the technology [4] and various design solutions have been tested in practice to either circumvent the need for reinforcement or integrate reinforcement after the concrete is printed [[16], [17], [18], [19], [20]]. As a result, several reinforcement techniques have been proposed, such as bar reinforcement [21], micro-cable reinforcement [22,23], fiber reinforcement into the cementitious material [[24], [25], [26]], steel reinforcement using robotic arc welding [27,28], and in-process mesh reinforcement [29]. For comprehensive details on the reinforcement strategies, refer to [30]. Nevertheless, these reinforcement strategies are still rudimentary in many instances.
This study focuses on bar reinforcement methods, where rebars are integrated with freshly deposited cementitious material. A few approaches can be found in the literature, for example, penetration of vertical bars through multiple printed layers [31,32], placement of horizontal bars into a printed layer along the printing direction and then covered by the next layer on top [33,34], and depositing around pre-installed bi-directional rebars [35]. However, in most approaches, the bonding between the rebar and concrete was compromised by the air void around the rebar [21,36]. To overcome this constraint, a large amount of trial and error is required, which is costly and time-consuming.
An approach to mitigate extensive experimental campaigns is to apply numerical models. In the context of 3D printing technologies, like Fused Filament Fabrication (FFF), Robocasting, and 3DCP, CFD modelling has been found to be very beneficial [[37], [38], [39], [40], [41], [42], [43], [44]]. The morphology of the deposited strands in FFF was studied by Comminal et al. [45]. Furthermore, Serdeczny et al. [46] addressed how to reduce the porosities and enhance the bonding between subsequent layers. Mollah et al. [[47], [48], [49]] studied ways to minimize the deformation and thereby stabilize layers printed by Robocasting, while for 3DCP, the geometrical shapes of the single- and multiple-deposited layers have been investigated in detail in [[50], [51], [52]].
This paper uses the CFD model and extends the preliminary results recently published in [53]. The model uses elasto-viscoplastic constitutive equations to approximate the rheology of the concrete. The CFD model is validated by comparison with a number of experiments, and the model is subsequently exploited to make an in-depth analysis of air void formation between rebars and concrete using the cross-sections of the deposited part and the calculated volume fraction of air voids. Different material properties, such as yield stress and plastic viscosity, and processing parameters, like the rebar diameter, nozzle-rebar distance, a geometric ratio (i.e., the distance from nozzle to the substrate divided by the nozzle diameter), as well as a speed ratio (i.e., the printing speed divided by the extrusion speed) are varied. Section 2 describes the methodology of the study, along with the experimental and numerical details. Next, Section 3 presents and discusses the results. Finally, Section 4 summarizes the results with the conclusion.

2. Methodology

2.1. Materials’ properties and 3DCP experiments

A fresh cement-based mortar was used to perform the 3DCP experiment around the rebars. The mortar includes a binder system with white cement CEM I 52.5 R-SR 5 (EA), limestone filler with sand of maximum particle size 0.5 mm, admixtures, and water. The binder was prepared with a 75 L Eirich Intensive Mixer Type Ro8W. The water to cement ratio was 0.39. The admixtures dosage (by weight of cement) was set at 0.1 % high-range water-reducing agent, 0.1 % viscosity-modifying agent, and 0.5 % hydration retarder.

The rheological characterization of the mortar was done using an Anton Paar rheometer MCR 502, as used in [50,54]. The rotational and oscillatory tests were performed with a vane-in-cup measuring device. The obtained flow curve of the mortar from the rotational rheometric tests, with a ramp-down controlled shear rate (CSR), was fitted by a linear regression to determine the yield stress τ0= 630 Pa and plastic viscosity ηP= 7.5 Pa·s. The oscillatory test showed that the constitutive behavior of the unyielded mortar had a factorized relationship between the storage modulus G′ and loss modulus G′′ within the linear viscoelastic (LVE) region, where G′= 200 kPa was captured. Therefore, the mortar’s rheology was modelled as a yield stress limited elasto-viscoplastic material, where the storage modulus is used as the linear elastic shear modulus of the unyielded mortar. Furthermore, the rheological characterization showed that the mortar exhibited time-independent rheological characteristics within the actual printing process, see [50] for more details.

The setup for 3DCP experiments around ribbed rebars is presented in Fig. 1. It comprised a 6-axis industrial robot (Fanuc R-2000iC/165F) with a custom-designed nozzle ∅20 mm (i.e., nozzle diameter, Dn= 20 mm) made by fused filament fabrication of ABS thermoplastic, cf. Fig. 1-a. The robot also included a progressive cavity pump (NETZSCH) equipped with a hopper and a long steel-wire rubber hose (cf. Refs. [50, 52] for details). A 25 mm thick plywood plate was used as the built substrate as seen in Fig. 1-b. The 1000 mm long rebars of diameter Dr= 8 and 12 mm were placed horizontally on top of the substrate at a distance Hr= 14 mm. The horizontal rebars were held in place by two vertical rebars with a height of 37 mm. The setup was used to print a structure of four successive layers of parallel strands around the rebars. Details on the printing toolpath around the rebars are illustrated in the subsections below.

Fig. 1. 3DCP experiment around rebars: (a) 6-axis robotic arm [50]; (b) plywood built platform with integrated rebars; (c) example of printing (picture is taken during printing of the third layer).

The extrusion nozzle was placed above the substrate with a nozzle height Dn/2 for the first layer, whereas for subsequent number of layers (Nl), the nozzle height was set at Nl∗Dn/2. Thus, the nominal height of a layer was h=Dn/2. The print was done with a material extrusion rate 0.91 dm3/min and nozzle speed 35 mm/s. An example of a physical print is presented in Fig. 1-c. After the prints hardened, cross-sections were collected to investigate the rebar-concrete bonding. The cross-section slices were taken at specific positions to analyze the print around the horizontal rebar and cross-shaped rebar (i.e., horizontal and vertical rebars). To avoid destroying the specimens while cutting them, the printed part were impregnated with epoxy resin in a vacuum chamber.

2.2. Computational models and governing equations

Three different CFD models are built. The first model only simulated the mortar flow to understand the void formation pattern without rebars. The last two models simulated the 3DCP experiment around rebars: one model simulated the mortar flow around the horizontal rebar, while the other considered the cross-shaped rebar. This subdivision enabled the CFD models to consider a smaller computational domain than if the two scenarios were combined.

The CFD models comprised a cylindrical nozzle, a solid-substrate, and an artificial solid component (at the top) within the computational domain of size 8.5D×6D×2D+4h as shown in Fig. 2 (top), where Model 1 excluded rebars (left), Model 2 included the horizontal rebar (middle), and Model 3 considered the cross-shaped rebar (right). The printing toolpath of the models are illustrated in Fig. 2. The toolpath for Models 1 and 2 are presented in 3D (left bottom figure), where the only difference was the presence of the horizontal rebar. The toolpath in 2D presented at the bottom right is for Model 3. For all the models, the toolpath of the extrusion nozzle kept a distance of Dnr from the axis of the nearest rebar. The lengths of the horizontal and vertical rebars were 50 and 40 mm, respectively. The other printing parameters were similar to the ones used in the experiment, cf. Section 2.1. Finally, the models were used to simulate four successive layers with a length of 125 mm. Note that the rebars are modelled as cylindrical solid objects (i.e., smooth rebar).

Fig. 2. Model geometry with the extrusion nozzle, substrate, integrated rebars, and computational domain (top) and toolpath (bottom).

The computational domain was meshed by a uniform Cartesian grid. A mesh sensitivity test was performed for different meshes with cell sizes 0.9, 1.0, and 1.1 μm. Even if the change in absolute size of the cells were small, the total number of cells within the domain was 1.1, 1.5, and 2.0 million, respectively. A cell size of 1.0 μm was chosen as that was found to be time-efficient and had a negligible effect on the accuracy of the results. The top plane of the domain was an inlet boundary, where the artificial solid component was defined in order to prevent material flow outside the nozzle orifice, cf. Fig. 2 (top). On the bottom plane, a wall boundary was applied to represent the solid substrate. The other planes were assigned continuative boundary conditions, but had no influence on the results. Furthermore, no-slip boundary conditions were applied between fluids and solids.

Table 1 lists the printing parameters and their values for each of the investigated cases. All the models and cases are simulated for 4 successive layers.

Table 1. Description of case IDs with printing parameters and accompanying values. The reference values (corresponding to the experimental print) are written in bold, while the parameter change for each case is highlighted by underlining the value.

The cementitious mortar flow was assumed transient and isothermal. Thus, the flow dynamics of the mortar are governed by the mass and momentum conservation equations of incompressible fluid:

where u is the velocity vector, ρ is the density, g=00−g is the gravitational acceleration vector, t is the time, p is the pressure, and σ is the deviatoric stress tensor.

The rheological behavior of the mortar was modelled by the following elasto-viscoplastic constitutive equation that represents σ as the sum of the deviatoric part of the viscous stress σV and elastic stress σE tensors; i.e.:

The deviatoric viscous stress tensor was predicted as:

is the deformation rate tensor, and T represents the transpose notation.

The deviatoric elastic stress tensor was modelled by the Hookean assumption of a small strain rate tensor E between each small time steps Δt=t−t0, to represent the elastic response of unyielded materials 

Ewhere G is the shear modulus and Et=Et0+ΔtDT is the incremental strain rate tensor approximated by integrating the deformation rate tensor over Δt.

The incremental representation of Eq. (5) can be written as:

 is the vorticity tensor. The first term of the left-hand side of Eq. (6) represents the change in stress at a fixed location in space. The change in stress due to advection and rotation of material particle is approximated by the second and third terms, respectively. The right-hand side takes into account the change in stress due to shearing.

The elastic stress tensor of the yielded material was approximated by imposing the yield stress τ0 limit as follows:

where σvM is the von Mises stress predicted as:

where IIσE∗=trσE∗2 is the second invariant of σE∗. The material was yielded when σvM exceeded the yield stress. The properties of the material used in the different models and cases are presented in Table 2. Note that the CFD model does not include the solidifications of the printed layers.

Table 2. Material properties.

2.3. Numerical method

The computational model was developed in the commercial CFD tool FLOW-3D® (V12.0; Flow Science Inc.) [55]. It uses the FAVOR technique (Fractional Area/Volume Obstacle Representation) to embed solid objects (i.e., the nozzle, rebars, substrate, etc.) in the computational domain. The computational domain was meshed with a Cartesian grid and discretized with the Finite Volume Method.
The governing equations of the mortar flow were solved by the implicit pressure-velocity solver GMRES (Generalized Minimum Residual) [[56], [57], [58]]. The predictions of pressure and forces near solid objects were modelled by the immersed boundary method [59]. The yield stress limited elasto-viscoplastic criterion was built in the software, and the elastic stress was calculated explicitly. An implicit technique, successive under-relaxation, was used to solve the viscous stress of the momentum equation (Eq. (2)). The free surface of the mortar was captured by the Volume-of-Fluid technique, see details in Ref. [60, 61]. The momentum advection was calculated explicitly by an upwind-difference technique and ensured first-order accuracy. The time step size was controlled dynamically with a stability limit in order to avoid numerical instabilities [55]. All the simulations were run with 20 cores on a high-performance computing cluster. The study was carried out with a first-order accuracy in both space and time in order to reduce the computational time, which was extensive due to the elasto-viscoplastic material model that computes both viscous and elastic stresses (e.g., the computational time was six days for model 1 case 1). In this regard, one should note that the model can simulate a multitude of scenarios simultaneous.

2.4. Results post-processing

The simulated results were processed in two steps. The first step was to show the cross-sectional shapes which were done in the post-processing tool FLOW-3D®POST, and the second step was to calculate the volume fraction of air voids inside the printed structure. The cross-sectional shapes were used to investigate the interior of the structure and the rebar-concrete bonding. The cross-sectional shapes were obtained in the plane parallel to the yz-plane at the middle of the layer’s length, as shown in Fig. 3-d. Fig. 3-c sketches the nominal positions of the air void creation in a four-layered structure. The positions of the air voids were defined as outer-bottom, bottom, mid, top, and outer-top in this study. The presence of the air void was quantified by calculating the volume fraction of air voids around the middle of the layer’s length. The calculation enabled to capture the presence of air voids for Model 3, i.e., around the vertical rebar, as seen in the dashed-box of Fig. 3-b.

Fig. 3. Post-processing of results; (a) introduction of volume sampling cuboid to calculate the volume fraction of air voids; (b) presence of air voids around the vertical rebar in the experiment; (c) schematic of air void creation; (d) cross-sectional shapes and void sampling area for the different models.

In order to calculate the volume fraction of air voids, a cuboid of size 20×25×3h mm3 was introduced to the CFD models as a volume sampling object, as seen in Fig. 3-a. Note that the size and position of the object were the same for all the models. The volume sampling was a three dimensional data collection tool built in the software that enabled calculating the amount of material as well as air void inside of it. Finally, the volume fraction of air voids VV was calculated as below:

3. Results and discussions

his section compares the simulated and experimental results of 3DCP around the rebars. Furthermore, it discusses the influence of different parameters on the air void formation in the cross-sectional shapes of printed parts and the volume fraction of air voids inside the structure. The parametric study includes the material properties (i.e., yield stress and plastic viscosity) and the printing properties (i.e., rebar diameter, rebar-nozzle distance, geometric ratio, and speed ratio).

3.1. Experiments and validation of the CFD model

The CFD models (Models 2 and 3) are compared and validated with the experiments. The results are presented in Fig. 4. Two rebar diameters, i.e., Dr= 8 and 12 mm, are taken into account, where the nozzle-rebar distances are 20 and 24 mm, respectively. The other printing and material parameters are kept constant in the experiments with different rebar diameters, see Table 1 (cases 1 and 3). In the case of the simulations, all parametric details are the same as implemented in the experiments except for the elastic shear modulus. The choice of shear modulus is subject to the analysis presented in Appendix A. The shear modulus is reduced to 100 kPa from the measured value (i.e., 200 kPa) to compare the simulated results with experiments. Furthermore, a shear modulus of 20 kPa is chosen for the parametric study in the later sections. These assumptions seem reasonable to avoid extensive computational time consumption since the differences found in void formation are limited (see Fig. A.1).

Fig. 4. 3DCP experiments (left column), simulations (middle column), and comparison (right column). (a) Horizontal rebar with Dr= 8 mm and Dnr= 20 mm; (b) horizontal rebar with Dr= 12 mm and Dnr= 24 mm; (c) cross-shaped rebar with Dr= 8 mm and Dnr= 20 mm; and (d) cross-shaped rebar with Dr= 12 mm and Dnr= 24 mm. The blue part in the experiments is epoxy resin. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The cross-sectional shapes of the 3DCP experiments around the horizontal rebar of
8 mm in Fig. 4-a illustrate the presence of a top air void as well as mid and bottom air voids positioned respectively above and below the horizontal rebar. The mid and bottom air voids are found to be smaller than the top one. This is due to the presence of the rebar that occupies the mortar’s flowable space as well as the deformation of the previously printed layers. For a detailed analysis of the deformation pattern, refer to [53]. The air voids at the top and bottom are significant when the rebar diameter is increased to 12 mm (Fig. 4-b). This is due to the fact that the nozzle-rebar distance was increased, which enhances the flowable space between the strands. In addition, the larger size of the rebar creates larger channels below and above itself, where the mortar of the second and third layers is forced to be squeezed into. The mid-air void is found to be absent as its area is fully occupied by the larger rebar. In the case of the cross-shaped rebars, the existence of the vertical rebar seems to restrict the merging of parallel strands, and therefore, the presence of air void content increases, cf. Fig. 4-c, d. This limitation is found to be pronounced for the larger rebar diameter with the larger nozzle-rebar distance.
The cross-sectional shapes of the simulations (middle column in Fig. 4) illustrate high accuracy predictions of the position and size of the air voids when compared with the experiments. Particularly, the models capture small details around the vertical rebar for both diameters. This can clearly be seen in the comparison of experiments and simulation, cf. right column in Fig. 4. A discrepancy is found in the strand’s width of the bottom layer as well as the shape of the printed part, specifically in the shape of strands of all the layers next to the vertical rebar and the height of the part for the smaller rebar diameter. These could be due to a combined effect of the idealized rheological model, as well as slight differences in the processing parameters, e.g. nozzle height above the printing surface, nozzle-rebar distance, as well as printing- and extrusion-speed. Note that the height of the vertical rebar in the experiments is a bit shorter than the one in the simulations (40 mm), although it does not influence the results.

3.2. Influence of materials properties

The influence of the material properties, yield stress and plastic viscosity, on the air void formation is presented in Fig. 5, Fig. 6, Fig. 7. The process parameters of Case 1 (cf. Table 1) are utilized.

Fig. 5. Air void formation in the cross-sections of the printed parts for different yield stress.

Fig. 6. Air void formation in the cross-sections of the printed parts for different plastic viscosity.

Fig. 7. Volume fraction of air voids for different models as a function of (a) yield stress and (b) plastic viscosity.

Yield stress

Fig. 5 presents cross-sectional shapes for different yield stress, 400, 630, and 800 Pa. It can be seen that Models 1 and 2 predict a top air void, while for Model 3 the two topmost strands to the left are not in contact with the vertical rebar. The cross-sections illustrate that an increased yield stress causes less deformation of the printed layers and create stands with less round shape. This is due to the reduced effective gap (i.e., the distance between the nozzle and previous printed layer), which results in a reduced air void content for most models as seen quantitatively in Fig. 7-a. This behavior is converse to conventional concrete castings where a more fluid material (e.g. self-compacting concrete) can lead to a lower void content. An exception to the observed behavior that a higher yield stress leads to less voids formation is seen in case of Model 2 with
800 Pa, where the top air void is slightly larger as compared to the one for
630 Pa. Another exception is that an outer bottom air void appears for Models 1 and 2 when increasing the yield stress to 800 Pa. Both exceptions are a consequence of the yield stress now restricting the flow in confined spaces, which illustrates that it is a non-trivial task to fully eradicate air voids only by increasing the yield stress.

Plastic viscosity

As the plastic viscosity is varied, cross-sections for Model 1 show a slight change in air voids, see Fig. 6. A mid-air void is produced when the plastic viscosity is small, while the two largest plastic viscosities only produce the top air void. This is due to the increase in extrusion pressure that leads to larger deformation of the printed layers when the plastic viscosity is increased, cf. details in ref. [47]. When integrating a horizontal rebar (see Model 2), the air void formation increases at higher plastic viscosities. This could be due to the fact that the sideway flow of the depositing material (i.e., y-velocity) is limited by the flow resistance that comes from both the larger plastic viscosity and the presence of the solid rebar. No noticeable change can be seen in the cross-sections of Model 3 for different plastic viscosities. The same findings are quantitatively highlighted in Fig. 7-b, which illustrates that the volume fraction of air voids is not influenced much by the plastic viscosity except for Model 2.

3.3. Influence of processing conditions

The influence of processing conditions such as rebar diameter, nozzle-rebar distance, geometric ratio, and speed ratio on air void formation is presented in terms of cross-sectional shapes and volume fraction of air void, see Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12. Models 1 to 3 are simulated with the reference material properties, cf. Table 2.

Fig. 8. Air void formation in the cross-sections of the printed parts for different rebar diameters.

Fig. 9. Air void formation in the cross-sections of the printed parts for different nozzle-rebar distances.

Fig. 10. Air void formation in the cross-sections of the printed parts for different geometric ratios.

Fig. 11. Air void formation in the cross-sections of the printed parts for different speed ratios.

Fig. 12. Volume fraction of air voids for different models as a function of (a) rebar diameter, (b) nozzle-rebar distance, (c) geometric ratio, and (d) speed ratio.

Rebar diameter

The influence of the rebar diameter on the air void formation is presented in Fig. 8, Fig. 12-a. Models 2 and 3 are simulated (Model 1 does not contain a rebar) for Cases 1, 2, and 3, cf. Table 1. Fig. 8 illustrates that the top air void appears almost constant for Model 2, whereas the air void below the rebar increases with an enlarged rebar diameter. Two phenomena with opposite effects on the void formation play a role in this regard. On the one hand, increasing the rebar size reduces the space that the strands need to occupy to fully merge and thereby eliminate voids. On the other hand, the resistance towards flow and merging of the parallel strands next to the rebar increases proportionally with the size of the reinforcement. The latter effect is dominating in this case. For Model 3, the air void formation also increases when increasing the reinforcement. In addition to the previously mentioned argument, this is due to the left strands having to flow longer to reach the vertical rebar (i.e., Dnr+Dr/2). Conversely, additional air voids take place on the right-hand side of the vertical rebar for the smallest rebar diameter, because the nozzle-rebar distance Dnr= 20 mm is kept constant. Fig. 12-a underlines quantitatively that the volume fraction of air voids reduces when the rebar diameter is small. The trend is more pronounced for the cross-shaped rebar (Model 3), but in absolute values the air voids are substantially less for Model 2.

Nozzle-rebar distance

Fig. 9, Fig. 12-b show the effect of different nozzle-rebar distances on the formation of air voids. All the models are simulated for Cases 1, 4, and 5, cf. Table 1. Fig. 9 shows that the presence of air voids is reduced when the nozzle-rebar distance reduces. This is because the flowable space around the rebars shrinks. Interestingly, no significant air voids are present in Model 1 and 2 when Dnr= 18 mm, see Fig. 12-b. Fig. 12-b also depicts that the trend is more pronounced for the cross-shaped rebar (Model 3). One should be careful though about decreasing the nozzle-rebar distance too much, as a ridge is forming on the top layer since material from the left strand starts to flow on top of the right strand (Fig. 9), which potentially could affect the final shape of the structure. In case of the cross-shaped rebar model, air voids are formed for all investigated Dnr. One could potentially with benefit reduce the Dnr further, but not more than the sum of half of the nozzle diameter (10 mm), the nozzle wall thickness (2.5 mm), and half of the rebar diameter (4 mm), i.e., 16.5 mm, otherwise the nozzle will collide with the rebar.

Geometric ratio

The effect of the geometric ratio on the formation of air voids is presented in Fig. 10, Fig. 12-c. The considered simulations are Cases 1, 6, and 7 cf. Table 1. Fig. 10 illustrates that decreasing the geometric ratio can reduce the presence of air voids in the cross-sections. This is because a smaller geometric ratio results in wider strands, which then occupy more of the flowable space around the rebars. Note that when Gr= 0.50, 0.45, and 0.40 the layer height is 10, 9, and 8 mm, respectively. No air voids are formed for Model 1 and 2 when Gr= 0.45, and 0.40. However, for the smallest ratio ridges are obtained on either side of the strands as clearly seen for the top layer. These ridges can as previously mentioned have a negative effect on the final shape of the structure. Consequently, Gr= 0.45 is preferable for these two models. In the case of Model 3, air voids are still present next to the vertical rebar, even for the smallest investigated geometric ratio. The volume fraction of air voids is approximately 1.5 %, see Fig. 12-c. The geometric ratio could be reduced further in order to decrease the air voids even more, but the ridges already form at Gr= 0.40. Consequently, it is not possible to fully eliminate air voids while at the same time avoiding ridges when only varying the geometric ratio for Model 3.

Speed ratio

Fig. 11, Fig. 12-d illustrate the formation of air voids for different speed ratios. The considered simulations are Case 1, Case 8, and Case 9, cf. Table 1. The three speed ratios are obtained by applying an extrusion speed of 48.4 mm/s, 51.5 mm/s, and 53.8 mm/s. Fig. 11 show that less air voids are formed when decreasing the speed ratio (i.e., higher extrusion speed). Reducing the speed ratio increases the cross-sectional area of the strands proportionally, thereby decreasing the air voids. Model 1 obtains no air voids for the two smallest ratios, and the same is almost the case for Model 2; only a very small air void is formed when Sr= 0.68, see Fig. 12-d. Model 3 forms air voids for all speed ratios. For the lowest speed ratio, the third strand to the left is in contact with the vertical rebar, but air voids are still formed around the horizontal reinforcement, which underlines the fact that it is difficult to fully eliminate air voids for the cross-shaped rebars.

3.4. Cross-shaped reinforcement

Based on the above analysis, it is clear that the air voids around the horizontal rebar can be eliminated by changing some of the processing conditions, such as the nozzle-rebar distance, geometric ratio, and speed ratio. However, it remains unsolved to fully omit the presence of air voids around the cross-shaped rebar, although the processing conditions can reduce the volume fraction of air voids. A parametric study was conducted by varying some combinations of processing conditions; however, the same conclusion was achieved that the air voids could not be fully eliminated. In order to solve this predicament, a new stepped toolpath is investigated (see Fig. 13) along with three different rebar geometries: 1) cylindrical rebars as in the previous analysis, see Fig. 14-a; 2) a squared horizontal rebar, cf. Fig. 14-b; and 3) cylindrical rebars with a smooth transition between them, see Fig. 14-c. In all scenarios, the speed ratio is 0.665, the size (i.e., diameter or side of square) of the rebars are 6 mm, and the horizontal rebar is placed at a height of 8 mm from the substrate. The other processing parameters are the same as for Case 2 and reference material properties are applied. For scenarios one and two small air voids are formed, but for scenario three air voids are eliminated, see Fig. 14. This numerical analysis illustrates that although it is difficult to get rid of air voids for the cross-shaped rebars, one can do it when carefully selecting the material- and processing-parameters and remembering to have a smooth transition between the rebars.

Fig. 13. New toolpath planning around the cross-shaped rebar.

Fig. 14. Simulated structure with new toolpath and different rebar geometries. (a) Cylindrical horizontal- and vertical-rebar, (b) square horizontal rebar and cylindrical vertical rebar, (c) cylindrical horizontal- and vertical-rebar with smooth transition. Note that Dr= 6 mm, Sr= 0.665, and Hr= 8 mm.

4. Conclusions

A CFD model was used to predict the morphology of strands and the formation of air voids around reinforcement bars when integrated with 3DCP. The model used an elasto-viscoplastic constitutive model to mimic the cementitious mortar flow. The CFD model was compared with experiments that constituted a horizontal and a cross-shaped rebar configuration. The results illustrated that the model with high-accuracy could predict the air void formation in the structures. The simulations had though slightly less wide bottom strands as compared to the experimental counterpart, which was attributed to small differences in material behavior and processing parameters.
The CFD model was exploited to investigate the effect of material properties on the air void formation. The results illustrated that by increasing the yield stress less air voids were formed due to the reduced effective gap. However, the air voids could not be eliminated as the increased yield stress also restricted the flow in confined spaces. In contrast to the effect of the yield stress, the void formation decreased somewhat when decreasing the plastic viscosity (although not enough to omit air voids fully).
The process parameters were found to have a substantial effect on the air void formation. The air void formation increased when increasing the rebar diameter, because the resistance towards flow around the reinforcement and thereby merging of the strands increased proportional with the size of the rebars. The air voids could be reduced and in some of the horizontal cases fully avoided by reducing the nozzle-rebar distance, but it could come with the expense of ridges (which could affect the final geometry of the structure), since material from one strand would flow on top of a previously deposited stand. Similarly, decreasing the geometric ratio was found to reduce the presence of air voids, because a smaller geometric ratio resulted in wider stands that occupied more of the space around the rebars. However, the smallest ratios also resulted in ridges. It was also found that less air voids were formed when decreasing the speed ratio, since the cross-sectional area of the strands increased proportionally, thereby occupying more space around the rebars.
By decreasing the nozzle-rebar distance, geometric ratio, and speed ratio, voids were omitted around the horizontal rebar, but air voids would still be introduced for the cross-shaped rebar. Those air voids could be eliminated by changing the toolpath and some processing parameters, as well as altering the geometry of the reinforcement joint to a smooth transition between the horizontal and vertical rebar. The results highlight that it is non-trivial to avoid air voids when integrating rebars in 3DCP, but that the CFD model is a very strong digital tool when it comes to securing a good bonding between the reinforcement and concrete.
A limitation of the CFD model is that with the current computational power it is not possible to simulate a full 3DCP structure. Nevertheless, the CFD model is powerful when it comes to understanding and optimizing printing strategies for individual reinforcement details. In future research, it would be interesting to exploit the CFD model to systematically investigate various reinforcement setups and based on the model results generate response surfaces or lookup tables, which can be coupled with the 3DCP toolpath software. This approach would have several benefits: 1) the computational heavy CFD model could run all scenarios in parallel, thereby minimizing the physical time spent on calculations; 2) one could avoid repetition of individual reinforcement studies; and 3) the response surfaces or lookup tables could in a computational light manner come up with printing strategies for all reinforcement details in a full 3DCP structure.

CRediT authorship contribution statement

Md Tusher Mollah: Conceptualization, Methodology, Investigation, Formal analysis, Visualization, Writing – original draft, Writing – review & editing. Raphaël Comminal: Conceptualization, Formal analysis, Writing – review & editing, Supervision. Wilson Ricardo Leal da Silva: Conceptualization, Formal analysis, Writing – review & editing, Resources. Berin Šeta: Investigation, Formal analysis, Writing – review & editing. Jon Spangenberg: Conceptualization, Investigation, Formal analysis, Writing – review & editing, Supervision, Resources, Project administration.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors would like to acknowledge the support of the Danish Council for Independent Research (DFF) | Technology and Production Sciences (FTP) (Contract No. 8022-00042B). Also, the authors would like to acknowledge the support of the Innovation Fund Denmark (Grant No. 8055-00030B: Next Generation of 3D-printed Concrete Structures and Grant no. 0223-00084B: ThermoForm – Robotic ThermoSetting Printing of Large-Scale Construction Formwork), Moreover, the support of FLOW-3D® regarding licenses is acknowledged.

Appendix A.

This analysis varies the shear modulus (i.e., 20, 50, and 100 kPa) in the case Dr= 8 mm, as seen in Fig. A.1, which presents the cross-sectional shapes (top) and the volume fraction of air voids (bottom). It can be seen that increasing the shear modulus slightly reduces the air void formation. This is because the larger shear modulus enhances the ability of the material to act against the shear deformation. However, an increase in shear modulus also extensively increases the computational time of solving the non-linear elastic response of the elasto-viscoplastic material. For example, the computational time of Model 3 is about 6, 12, and 18 days for a shear modulus of 20, 50, and 100 kPa, respectively. Therefore, the shear modulus is reduced to 100 kPa from the measured value (i.e., 200 kPa) to compare the simulated results with experiments. Furthermore, the shear modulus 20 kPa is chosen for the parametric study in 3.2 Influence of materials properties, 3.3 Influence of processing conditions, 3.4 Cross-shaped reinforcement. These assumptions seem reasonable to avoid extensive computational time consumption since the differences found in Fig. A.1 are not substantial.

Fig. A.1. Air void formation in the cross-sections of the printed parts (top) and the volume fraction of air voids (bottom) for different shear modulus.

References

[1] B. Khoshnevis, R. Russel, H. Kwon, S. Bukkapatnam
Contour crafting–a layered fabrication technique
Spec. Issue IEEE Robot. Autom. Mag., 8 (2001), pp. 33-42
[2] R.A. Buswell, W.R.L. da Silva, F.P. Bos, H.R. Schipper, D. Lowke, N. Hack, H. Kloft, V. Mechtcherine, T. Wangler, N. Roussel
A process classification framework for defining and describing Digital Fabrication with Concrete
Cem. Concr. Res., 134 (2020), Article 106068, 10.1016/j.cemconres.2020.106068
[3] D. Asprone, C. Menna, F.P. Bos, T.A.M. Salet, J. Mata-Falcón, W. Kaufmann
Rethinking reinforcement for digital fabrication with concrete
Cem. Concr. Res., 112 (2018), pp. 111-121, 10.1016/j.cemconres.2018.05.020
[4] B. Khoshnevis
Automated construction by contour crafting—related robotics and information technologies
Autom. Constr., 13 (2004), pp. 5-19, 10.1016/j.autcon.2003.08.012
[5] R. Hague, I. Campbell, P. Dickens
Implications on design of rapid manufacturing
Proc. Inst. Mech. Eng. C J. Mech. Eng. Sci., 217 (2003), pp. 25-30, 10.1243/095440603762554587
[6] V. Mechtcherine, V.N. Nerella, F. Will, M. Näther, J. Otto, M. Krause
Large-scale digital concrete construction – CONPrint3D concept for on-site, monolithic 3D-printing
Autom. Constr., 107 (2019), Article 102933, 10.1016/j.autcon.2019.102933
[7] S. Lim, R.A. Buswell, T.T. Le, S.A. Austin, A.G.F. Gibb, T. Thorpe
Developments in construction-scale additive manufacturing processes
Autom. Constr., 21 (2012), pp. 262-268, 10.1016/j.autcon.2011.06.010
[8] Meet ICON’s next generation Vulcan construction system
ICON
(n.d.)
https://www.iconbuild.com/vulcan, Accessed 10th Jul 2022
[9] FABRICANT D’IMPRIMANTES 3D POUR LA CONSTRUCTION AUTOMATISÉE DE BÂTIMENTS, Constr.-3D
(n.d.)
https://www.constructions-3d.com/, Accessed 9th Jul 2022
[10] 3dWasp, 3dWasp
(n.d.)
https://www.3dwasp.com/, Accessed 12th Aug 2022
[11] G. Ma, R. Buswell, W.R. Leal da Silva, L. Wang, J. Xu, S.Z. Jones
Technology readiness: a global snapshot of 3D concrete printing and the frontiers for development
Cem. Concr. Res., 156 (2022), Article 106774, 10.1016/j.cemconres.2022.106774
[12] L. Wang, G. Ma, T. Liu, R. Buswell, Z. Li
Interlayer reinforcement of 3D printed concrete by the in-process deposition of U-nails
Cem. Concr. Res., 148 (2021), Article 106535, 10.1016/j.cemconres.2021.106535
[13] S. Kristombu Baduge, S. Navaratnam, Y. Abu-Zidan, T. McCormack, K. Nguyen, P. Mendis, G. Zhang, L. Aye
Improving performance of additive manufactured (3D printed) concrete: a review on material mix design, processing, interlayer bonding, and reinforcing methods
Structures., 29 (2021), pp. 1597-1609, 10.1016/j.istruc.2020.12.061
[14] L. Gebhard, J. Mata-Falcón, A. Anton, B. Dillenburger, W. Kaufmann
Structural behaviour of 3D printed concrete beams with various reinforcement strategies
Eng. Struct., 240 (2021), Article 112380, 10.1016/j.engstruct.2021.112380
[15] P. Feng, X. Meng, H. Zhang
Mechanical behavior of FRP sheets reinforced 3D elements printed with cementitious materials
Compos. Struct., 134 (2015), pp. 331-342, 10.1016/j.compstruct.2015.08.079
[16] F.P. Bos, C. Menna, M. Pradena, E. Kreiger, W.R.L. da Silva, A.U. Rehman, D. Weger, R.J.M. Wolfs, Y. Zhang, L. Ferrara, V. Mechtcherine
The realities of additively manufactured concrete structures in practice
Cem. Concr. Res., 156 (2022), Article 106746, 10.1016/j.cemconres.2022.106746
[17] H. Kloft, M. Empelmann, N. Hack, E. Herrmann, D. Lowke
Reinforcement strategies for 3D-concrete-printing
Civ. Eng. Des., 2 (2020), pp. 131-139, 10.1002/cend.202000022
[18] F.P. Bos, Z.Y. Ahmed, R.J.M. Wolfs, T.A.M. Salet
3D printing concrete with reinforcement
D.A. Hordijk, M. Luković (Eds.), High Tech Concr. Technol. Eng. Meet, Springer International Publishing, Cham (2018), pp. 2484-2493, 10.1007/978-3-319-59471-2_283
[19] J. Xiao, G. Ji, Y. Zhang, G. Ma, V. Mechtcherine, J. Pan, L. Wang, T. Ding, Z. Duan, S. Du
Large-scale 3D printing concrete technology: current status and future opportunities
Cem. Concr. Compos., 122 (2021), Article 104115, 10.1016/j.cemconcomp.2021.104115
[20] Z. Wu, A.M. Memari, J.P. Duarte
State of the art review of reinforcement strategies and technologies for 3D printing of concrete
Energies., 15 (2022), p. 360, 10.3390/en15010360
[21] B. Baz, G. Aouad, P. Leblond, O. Al-Mansouri, M. D’hondt, S. Remond
Mechanical assessment of concrete – steel bonding in 3D printed elements
Constr. Build. Mater., 256 (2020), Article 119457, 10.1016/j.conbuildmat.2020.119457
[22] G. Ma, Z. Li, L. Wang, G. Bai
Micro-cable reinforced geopolymer composite for extrusion-based 3D printing
Mater. Lett., 235 (2019), pp. 144-147, 10.1016/j.matlet.2018.09.159
[23] Z. Li, L. Wang, G. Ma
Mechanical improvement of continuous steel microcable reinforced geopolymer composites for 3D printing subjected to different loading conditions
Compos. Part B, 187 (2020), Article 107796, 10.1016/j.compositesb.2020.107796
[24] A.R. Arunothayan, B. Nematollahi, R. Ranade, S.H. Bong, J. Sanjayan
Development of 3D-printable ultra-high performance fiber-reinforced concrete for digital construction
Constr. Build. Mater., 257 (2020), Article 119546, 10.1016/j.conbuildmat.2020.119546
[25] B. Nematollahi, P. Vijay, J. Sanjayan, A. Nazari, M. Xia, V. Naidu Nerella, V. Mechtcherine
Effect of polypropylene fibre addition on properties of geopolymers made by 3D printing for digital construction
Materials., 11 (2018), p. 2352, 10.3390/ma11122352
[26] Y. Zhang, F. Aslani
Development of fibre reinforced engineered cementitious composite using polyvinyl alcohol fibre and activated carbon powder for 3D concrete printing
Constr. Build. Mater., 303 (2021), Article 124453, 10.1016/j.conbuildmat.2021.124453
[27] V. Mechtcherine, J. Grafe, V.N. Nerella, E. Spaniol, M. Hertel, U. Füssel
3D-printed steel reinforcement for digital concrete construction – manufacture, mechanical properties and bond behaviour
Constr. Build. Mater., 179 (2018), pp. 125-137, 10.1016/j.conbuildmat.2018.05.202
[28] J. Müller, M. Grabowski, C. Müller, J. Hensel, J. Unglaub, K. Thiele, H. Kloft, K. Dilger
Design and parameter identification of wire and arc additively manufactured (WAAM) steel bars for use in construction
Metals., 9 (2019), p. 725, 10.3390/met9070725
[29] T. Marchment, J. Sanjayan
Mesh reinforcing method for 3D concrete printing
Autom. Constr., 109 (2020), Article 102992, 10.1016/j.autcon.2019.102992
[30] V. Mechtcherine, R. Buswell, H. Kloft, F.P. Bos, N. Hack, R. Wolfs, J. Sanjayan, B. Nematollahi, E. Ivaniuk, T. Neef
Integrating reinforcement in digital fabrication with concrete: a review and classification framework
Cem. Concr. Compos., 119 (2021), Article 103964, 10.1016/j.cemconcomp.2021.103964
[31] T. Marchment, J. Sanjayan
Bond properties of reinforcing bar penetrations in 3D concrete printing
Autom. Constr., 120 (2020), Article 103394, 10.1016/j.autcon.2020.103394
[32] T. Marchment, J. Sanjayan
Lap joint reinforcement for 3D concrete printing
J. Struct. Eng., 148 (2022), Article 04022063, 10.1061/(ASCE)ST.1943-541X.0003361
[33] V2 Vesta Beton-3D-Drucker baut kleines Haus
(n.d.)
https://3druck.com/drucker-und-produkte/v2-vesta-beton-3d-drucker-baut-kleines-haus-2846225/, Accessed 27th Oct 2022
[34] Traditional reinforcement in 3D concrete printed structures, Eindh
Univ. Technol. Res. Portal
(n.d.)
https://research.tue.nl/en/studentTheses/traditional-reinforcement-in-3d-concrete-printed-structures, Accessed 31st May 2023
[35] M. Classen, J. Ungermann, R. Sharma
Additive manufacturing of reinforced concrete—development of a 3D printing technology for cementitious composites with metallic reinforcement
Appl. Sci., 10 (2020), p. 3791, 10.3390/app10113791
[36] L. Gebhard, L. Esposito, C. Menna, J. Mata-Falcón
Inter-laboratory study on the influence of 3D concrete printing set-ups on the bond behaviour of various reinforcements
Cem. Concr. Compos., 133 (2022), Article 104660, 10.1016/j.cemconcomp.2022.104660
[37] R. Comminal, J.H. Hattel, J. Spangenberg
Numerical simulations of planar extrusion and fused filament fabrication of non-Newtonian fluids
Annu. Trans. Nord. Rheol. Soc, 25 (2017)
[38] M.P. Serdeczny, R. Comminal, D.B. Pedersen, J. Spangenberg
Experimental validation of a numerical model for the strand shape in material extrusion additive manufacturing
Addit. Manuf., 24 (2018), pp. 145-153, 10.1016/j.addma.2018.09.022
[39] H. Xia, J. Lu, G. Tryggvason
A numerical study of the effect of viscoelastic stresses in fused filament fabrication
Comput. Methods Appl. Mech. Eng., 346 (2019), pp. 242-259, 10.1016/j.cma.2018.11.031
[40] M.P. Serdeczny, R. Comminal, Md.T. Mollah, D.B. Pedersen, J. Spangenberg
Numerical modeling of the polymer flow through the hot-end in filament-based material extrusion additive manufacturing
Addit. Manuf., 36 (2020), Article 101454, 10.1016/j.addma.2020.101454
[41] M.T. Mollah, M.P. Serdeczny, R. Comminal, B. Šeta, M. Brander, J. Spangenberg
A numerical investigation of the inter-layer bond and surface roughness during the yield stress buildup in wet-on-wet material extrusion additive manufacturing
2022 Summer Top. Meet. Adv. Precis. Addit. Manuf, American Society for Precision Engineering (2022), pp. 131-134
[42] M.P. Serdeczny, R. Comminal, M.T. Mollah, D.B. Pedersen, J. Spangenberg
Viscoelastic simulation and optimisation of the polymer flow through the hot-end during filament-based material extrusion additive manufacturing
Virtual Phys. Prototyp., 17 (2022), pp. 205-219, 10.1080/17452759.2022.2028522
[43] B. Šeta, M.T. Mollah, V. Kumar, D.K. Pokkalla, S. Kim, A.A. Hassen, J. Spangenberg
Modelling fiber orientation during additive manufacturing-compression molding processes
Proc. 33rd Annu. Int. Solid Free. Fabr. Symp, The University of Texas at Austin (2022), pp. 906-919
[44] M.T. Mollah, A. Moetazedian, A. Gleadall, J. Yan, W.E. Alphonso, R. Comminal, B. Šeta, T. Lock, J. Spangenberg
Investigation on corner precision at different corner angles in material extrusion additive manufacturing: an experimental and computational fluid dynamics analysis
Solid Free. Fabr. Symp. 2022 33rd Annu. Meet, The University of Texas at Austin (2022), pp. 872-881
[45] R. Comminal, M.P. Serdeczny, D.B. Pedersen, J. Spangenberg
Numerical modeling of the strand deposition flow in extrusion-based additive manufacturing
Addit. Manuf., 20 (2018), pp. 68-76, 10.1016/j.addma.2017.12.013
[46] M.P. Serdeczny, R. Comminal, D.B. Pedersen, J. Spangenberg
Numerical simulations of the mesostructure formation in material extrusion additive manufacturing
Addit. Manuf., 28 (2019), pp. 419-429, 10.1016/j.addma.2019.05.024
[47] M.T. Mollah, R. Comminal, M.P. Serdeczny, D.B. Pedersen, J. Spangenberg
Stability and deformations of deposited layers in material extrusion additive manufacturing
Addit. Manuf. (2021), Article 102193, 10.1016/j.addma.2021.102193
[48] M.T. Mollah, R. Comminal, M.P. Serdeczny, D.B. Pedersen, J. Spangenberg
Numerical predictions of bottom layer stability in material extrusion additive manufacturing
JOM. (2022), 10.1007/s11837-021-05035-9
[49] M.T. Mollah, R. Comminal, M.P. Serdeczny, B. Šeta, J. Spangenberg
Computational analysis of yield stress buildup and stability of deposited layers in material extrusion additive manufacturing
Addit. Manuf., 71 (2023), Article 103605, 10.1016/j.addma.2023.103605
[50] R. Comminal, W.R. Leal da Silva, T.J. Andersen, H. Stang, J. Spangenberg
Modelling of 3D concrete printing based on computational fluid dynamics
Cem. Concr. Res., 138 (2020), Article 106256, 10.1016/j.cemconres.2020.106256
[51] R. Comminal, W.R.L. Da Silva, T.J. Andersen, H. Stang, J. Spangenberg
Influence of processing parameters on the layer geometry in 3D concrete printing: experiments and modelling
F.P. Bos, S.S. Lucas, R.J.M. Wolfs, T.A.M. Salet (Eds.), Second RILEM Int. Conf. Concr. Digit. Fabr, Springer International Publishing, Cham (2020), pp. 852-862, 10.1007/978-3-030-49916-7_83
[52] J. Spangenberg, W.R. Leal da Silva, R. Comminal, M.T. Mollah, T.J. Andersen, H. Stang
Numerical simulation of multi-layer 3D concrete printing
RILEM Tech. Lett., 6 (2021), pp. 119-123, 10.21809/rilemtechlett.2021.142
[53] J. Spangenberg, W.R. Leal da Silva, M.T. Mollah, R. Comminal, T. Juul Andersen, H. Stang
Integrating reinforcement with 3D concrete printing: experiments and numerical modelling
R. Buswell, A. Blanco, S. Cavalaro, P. Kinnell (Eds.), Third RILEM Int. Conf. Concr. Digit. Fabr., Springer International Publishing, Cham (2022), pp. 379-384, 10.1007/978-3-031-06116-5_56
[54] N. Ranjbar, M. Mehrali, C. Kuenzel, C. Gundlach, D.B. Pedersen, A. Dolatshahi-Pirouz, J. Spangenberg
Rheological characterization of 3D printable geopolymers
Cem. Concr. Res., 147 (2021), Article 106498, 10.1016/j.cemconres.2021.106498
[55] FLOW-3D® Version 12.0 [Computer software]
https://www.flow3d.com (2019), Accessed 10th May 2022
[56] S.F. Ashby, T.A. Manteuffel, P.E. Saylor
A taxonomy for conjugate gradient methods
SIAM J. Numer. Anal., 27 (1990), pp. 1542-1568, 10.1137/0727091
[57] R. Barrett, M. Berry, T.F. Chan, J. Demmel, J.M. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H.V.D. Vorst
Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods
(1994)
[58] Y. Saad
Iterative Methods for Sparse Linear Systems
The PWS Series in Computer Science (1996)
(http://books.google.com/books?id=jLtiQgAACAAJ)
[59] Modeling Capabilities- Immersed Boundary Method | FLOW-3D | Solving The World’s Toughest CFD Problems, FLOW Sci
(n.d.)
https://www.flow3d.com/modeling-capabilities/immersed-boundary-method/, Accessed 21st Aug 2021
[60] C.W. Hirt, B.D. Nichols
Volume of fluid (VOF) method for the dynamics of free boundaries
J. Comput. Phys., 39 (1981), pp. 201-225, 10.1016/0021-9991(81)90145-5
[61] R. Comminal, J. Spangenberg, J.H. Hattel
Cellwise conservative unsplit advection for the volume of fluid method
J. Comput. Phys., 283 (2015), pp. 582-608, 10.1016/j.jcp.2014.12.003