Experimental and numerical evaluation of multi-directional compressive and flexure behavior of three-dimensional printed concrete

Abstract

Three-dimensional concrete printing (3DCP) can proliferate the industrialization of the construction sector, which is notoriously conservative and indolent toward changes. However, the mechanical behavior of 3DCP should be characterized and modeled considering the interfaces when its performance is thoroughly compared to that of the existing concrete construction methods. This study presents an experimental and numerical investigation of uniaxial compression and three-point bending (TPB) tests on extruded 3DCP beams in different loading directions. The orientation of translational and depositional interfaces with respect to the direction of loading influenced the strength. Both the elastic and post-damage behavior of the 3DCP specimens were compared with those of the conventionally cast specimen under quasi-static loading conditions. Despite the higher compressive strength of the casted specimen, the flexural strength of the 3DCP specimens was higher. This study employed the finite element and cohesive zone models of the appropriate calibrated traction-separation law to model fracture in the notched TPB specimens. Furthermore, the real-time acoustic emission test revealed the nature of failure phenomenon of three-dimensional-printed specimens under flexion, and accordingly, the cohesive law was chosen. The predicted load−displacement responses are in good agreement with the experimental results. Finally, the effects of cohesive thickness and notch shape on the performance under bending were explored through parametric studies.

Graphical abstract

Keywords

three-dimensional printing / anisotropy / flexure / compression / cohesive zone model / finite element model

1 Introduction

Extrusion-based three-dimensional concrete printing (3DCP) is an additive manufacturing process that employs the computer-controlled placement of extruded cementitious mortar to realize physical objects in a layer-by-layer fashion [1]. This automated method of constructing concrete has many advantages over the conventional construction method [2–6]. However, the layer-wise nature of concrete construction results in heterogeneous nature with anisotropic behavior [7,8]. The formation of random voids in the layer and in the interface of several layers results in weak structural capability ultimately. The interface between the layers likely influences the hardened properties of concrete components. Buswell et al. [2] investigated the relationship of freshly made and hardened paste, mortar, and concrete material properties on the final geometry of the fabricated part. Findings are categorized according to the construction application to produce a conglomerate of problems that pinpoint the range of upcoming research inquiry in this developing topic. Panda et al. [9] examined the three-dimensional (3D)-printed geopolymer mortar’s tensile bond strength considering the printing time gap, nozzle speed, and nozzle standoff. For the same batch of materials, it was discovered that a longer interval between layers weakens the material; however, a slower printing speed and closer nozzle standoff distance have a stronger effect.

Meurer and Classen [10] developed a simple experimental approach to evaluate the flexural, tensile, shear, compressive, and fracture energy of the 3D-printed concrete and mortars. Their experimental results demonstrated an anisotropic material behavior. In addition, they verified that interval duration has a significant effect on the attributes of the material under investigation, severely reducing strength and fracture energy for longer interval times. Liu et al. [11] investigated the 3DCP specimens with and without reinforcement considering varied pitch angles. When steel fibers were included in the specimens during the printing process, it was discovered that the effect of pitch angle on the impact resistance was more pronounced than in the specimens without fibers. The fiber-reinforced specimens exhibited a mixed-mode complex cracking process that included micro-cracks, branching, twisting, deflection, and twisting.

In a full-scale bench structure printed by Le et al. [12], 30% lower compressive strength was reported compared to that of the casted concrete structure. Panda et al. [13] explored the directional properties of 3DCP geo-polymer mixed with different fiber sizes and concluded that the addition of short glass fibers improved the concrete flexural strength. However, the flexural and tensile strength indicated a material anisotropic behavior. Nerella et al. [14] studied the effect of layer interface on the mechanical properties of the concrete element manufactured using the extrusion-based 3D-printing method. They concluded that the uniaxial compression (UC) test could not accurately evaluate the interface strength; however, it presented as an indicator of specimen anisotropic behavior when conducted in all three principal material directions. Xiao et al. [15] performed the finite element (FE) analysis under flexural and compression tests and showed that the printing parameters dominate the anisotropic mechanical characteristics of the 3D-printed concrete. Owing to the layer-wise extrusion approach, small voids are likely to occur between the extruded layers, resulting in an anisotropic behavior [16,17]. Wolfs et al. [18] studied the influence of various process parameters on the mechanical performance of 3D-printed concrete and stated that the interface strength was essential for the structural integrity of the printed specimen. In general, a specimen failure occurs in the form of interfacial shear slip, interlayer cracking owing to the shear and compression, layer delamination, crushing, or interfacial tensile cracking [19]. In addition, the failure modes for 3D-printed specimens were confirmed experimentally.

With the advent of sophisticated computational resources, the prediction of the failure response of materials using the finite element method (FEM) has become simpler [20]. Nguyen-Van et al. [21] identified and covered some of the most significant issues for modeling 3DCP properties. Xiao et al. [15] used the concrete damage plasticity (CDP) model along with the traction-separation law (TSL) in the commercial ABAQUS^® software package to study the effect of interface bond properties in the flexural and compression tests. In addition, the FE model of 3D-printed concrete accounted for the interface bond characteristics using TSL, and the results were validated using the experimental data. Jasim et al. [22] validated the experimental findings obtained by the other authors using the CDP model in the FE tool (ABAQUS), and the findings satisfactorily agreed with the deep beams. van den Heever et al. [23] used the anisotropic Rankine-Hill continuum multi-surface plasticity FE model to simulate the four-point bending vs crack mouth opening displacement of the 3DCP specimen; the response satisfactorily agreed with the modeling results. Yang et al. [24] found that the FE model of the reinforced concrete (RC) arch based on the CDP model could reasonably mimic the cracking pattern and load−deflection curves within an error of 10%. Peng et al. [25] investigated the behavior of RC structures under different loading conditions. They verified the applicability of seven interaction methods in LS-DYNA and relevant parameters by simulating the pull-out and three-point flexural tests based on the calibrated concrete model. The validations of the concrete model and corresponding parameters were verified through a comparison with the empirical formulae, related specifications, and test results. Liu et al. [26] used CDP to simulate the performance of some cementitious composites. Experiments with cementitious composites mixed on single-material 3DCP cubes were performed. The stress−strain data from the experiment was exported to the FE model to evaluate the numerical performance of multi-layered concrete cubes. Bedon et al. [27] used a cohesive zone model (CZM)-based FE numerical model to investigate the pull-out performance of a timber structure and validate previous experiments. Anas et al. [28] used a number of numerical damage models to simulate the nonlinear response on various slabs and experimentally validated their findings. Nguyen-Van et al. [29] simulated the experimental results using micromechanical modeling, and both results satisfactorily agreed in terms of the stress−strain curve and fracture patterns. The performance of 3D-printed beams that have a unidirectional printing pattern or a modest pitch angle is superior to those with bigger pitch angles when subjected to bending. Nguyen-Van et al. [30] employed a nonlinear dynamic FE model. The simulation and experimental impact forces agreed, notwithstanding a slight difference in the peak values. In addition, the damage pattern in the simulation demonstrated an overall consistency between the experimental and numerical results. To the best of our knowledge, the FE and cohesive law models appropriate for the 3DCP parts have not been combined to simulate the load vs. displacement behavior of TPB specimens under different directions of loading. Moreover, the effect of cohesive layer thickness and shape of the starter crack on the predicted mechanical response is investigated within the FE framework.

The existing studies have proved that the anisotropic behavior is obvious for 3DCP in the hardened state owing to the presence of multiple interfaces and internal porosity. Although mechanical properties of 3DCP have been studied, the numerical modeling of printed concrete behavior considering interlayer properties has not been completely explored. Therefore, in this study, an FE model of 3DCP considering the interfacial bond properties is proposed and experimentally verified. Moreover, the UC and three-point bending (TPB) tests were carried out to capture the load−displacement and fracture behavior of 3D-printed concrete along three loading directions. In addition, an acoustic emission (AE) sensor was used for the first time and based on the AE data, and an appropriate TSL was chosen and calibrated to predict the concrete response under flexural loading and gain an insight into the salient damage and fracture mechanism.

2 Materials and experimental methods

2.1 Three-dimensional printable concrete mix

A 3D printable concrete mixture was prepared by blending an equal proportion of fine aggregate with the total binder. In this study, the binder contained 43-grade ordinary Portland cement complying with IS 8112: 2013, Class F Fly ash complying with IS 3812 Part 1: 2013, and densified silica fume that had a specific gravity of 3.15, 2.08, and 2.17, respectively. The specific proportion (wt.%) of different binder contents is listed in Tab.1. The fine aggregate used in this study was a regular river sand (Zone-II) with a maximum particle size of 2.36 mm, specific gravity of 2.65, and fineness modulus of 2.75. To tweak the mixture for smooth extrusion and printing, 0.1 wt.% of Methyl cellulose-based viscosity modifying admixture (VMA) was added to the binder. All these dry materials were mixed in a laboratory-scale Hobart mixture and rotated at 107 revolutions per minute for 1 min to achieve a homogeneous mixture. Subsequently, 75% of the total water was added to the dry mixture at 198 revolutions per minute for 2 min, followed by stirring for 30 s. At the final stage of mixing, the rest of water (25%) was added and mixed at 361 revolutions per minute for 1 min to gain a consistent mix appropriate for printing.

Tab.1 Mixture proportion used for the printing

material	ratio
sand/binder	1.0
cement/binder	0.7
fly ash/binder	0.2
silica fume/binder	0.1
water/binder	0.36
VMA/binder	0.01

2.2 Three-dimensional concrete printer and sample preparation

The concrete printing was performed using an extrusion-based 3D printing machine with a build volume of 1 m × 1 m × 1 m in the x-, y-, and z-direction, respectively, as shown in Fig.1(a). The printing was performed with a circular nozzle (diameter = 20 mm) with a stand-off distance of 10 mm from the print head to have a continuous rectangular filament. The schematic representation of the 3D-printed slab and 3DCP process flow is depicted in Fig.1(b) and Fig.2, respectively. As shown in Fig.1(b), translational and depositional interfaces exist. The two interfaces should be considered in the analysis as the translational interface is generated before the depositional interface. Xiao et al. [15] studied the effect of these two types of interfaces on the mechanical behavior of 3DCP structures. The mechanical properties of both interfaces are not the same owing to differences in the curing times.

Fig.1 (a) Laboratory-scale 3-axis concrete printer; (b) schematic of printing directions.

Full size|PPT slide

Fig.2 Process flow of 3D concrete printing.

Full size|PPT slide

2.3 Uniaxial compression test

The specimen for the UC test was prepared according to ASTM C109 [31]. The 3DCP specimen for the UC test is shown along with its geometry specifications and layer orientation in Fig.3. Three loading directions (x-, y-, and z-direction) are considered in this study owing to anisotropy.

Fig.3 Compressive test specimen with dimension: (a) x-direction; (b) y-direction; (c) z-direction.

Full size|PPT slide

The desired geometric specifications of the specimen for the UC test were extracted from a 3DCP slab with a dimension of 700 mm × 100 mm × 70 mm using a concrete cutter. For each direction, five specimens were prepared to guarantee the repeatability of the mechanical behavior. Before testing, it was confirmed that the specimen and load cell surfaces were perfectly intact without any gap, and tests were performed at a constant loading rate of 100 kN/min.

2.4 Three-point bending test

The specimen for the TPB test was prepared according to ASTM C880 [32], where the dimension of the beam specimen was 160 mm × 40 mm × 40 mm. In this study, a notch was introduced directly along the line of action of loading at the bottom of the specimen, as presented in Ref. [10]. The notch had a length and width of 10 and 3 mm, respectively; it was prepared at the mid-distance of the prepared specimen. The schematic of the 3DCP specimen for the TPB test, along with the specimen geometry and layer orientation is shown in Fig.4. In Fig.4(b), the notch is shown at the middle of the layer. In the simulation, the location and dimensions of the notch correspond to the experimentally observed position and geometrical dimensions. The notch prepared for the TPB specimen loaded in the y-direction did not coincide along the interface. In the simulation, following the physical TPB specimen, the notch was made slightly offset from the interface. The TPB test was performed to explore the flexural and fracture behavior. For each direction, six specimens were prepared to guarantee the repeatability. The span between the two lower supports was fixed at 120 mm and the test was performed at a loading rate of 2600 N/min according to the ASTM standard. The mechanical characterization studies including UC test and TPB were conducted on 250 kN BISS universal testing machine and 100 kN SHIMADZU universal testing machine, respectively, after a typical curing period of 28 d, as shown in Fig. S1 in the Electronic Supplementary Material.

Fig.4 Schematic of TPB specimen with dimensions: (a) x-direction; (b) y-direction; (c) z-direction.

Full size|PPT slide

The 3DCP specimens for the TPB tests were extracted using a concrete cutter from the 3DCP slab with dimensions of 700 mm × 300 mm × 70 mm. As shown in Fig.5, many studies have shown that when specimens marked y and $y^{\prime }$ were loaded along the y-direction exhibited a negligible difference in the mechanical behavior [15]. Therefore, only the y specimen was considered in this study.

Fig.5 Four directions for the flexure specimen in the 3DCP slab.

Full size|PPT slide

2.5 Determination of elastic modulus and Poisson’s ratio

The elastic modulus of the 3DCP concrete and conventionally cast specimens are determined using an empirical formula provided by the American Concrete Institute (ACI). Most of the empirical formulae provided by various codes are a function of the compressive strength. According to ACI 318-14 Subsection 19.2.2, the modulus of elasticity for normal-weight concrete [33] is evaluated using

(1)$E_{\mathrm{c}}=4700\sqrt{{\sigma }_{\mathrm{c}}},$

where ${\sigma }_{\mathrm{c}}$ represents the compressive strength. To evaluate Young’s modulus, the compressive test in all three orthogonal directions of the 3DCP and conventionally cast specimens was performed. Poisson’s ratio was considered $0.2$ [34] in the numerical analysis.

2.6 Acoustic emission sensor test

A commercial piezoelectric sensor system, AEwin (physical acoustics) provided by MISTRAS Private Limited was used to achieve AE during the test in terms of voltage and time. The AE sensor had a DAQ frequency response in the range of 1 kHz–1 MHz whereas the resonance frequency was 300 kHz. The sensor had a peak sensitivity of 61 dB and a temperature range from 60 to 170 °C. Moreover, details of the sensor were presented in Ref. [35]. The pencil-lead break and wave velocity tests were performed for calibration and placement of the sensor. The sensor was placed at the bottom of the specimen approximately10 mm away from the notch.

3 Results and discussion

3.1 Uniaxial compression test

Fig.6 shows the mean compressive strength of the conventionally cast and 3DCP specimens loaded in three directions. As expected, the compressive strength of the conventionally cast specimen is comparatively higher than that of the 3DCP specimens in all directions. The 3DCP specimens loaded in the x-direction has the lowest compressive strength, whereas the specimen loaded in the y-direction has the highest compressive strength. Thus, the material strength dependence on the direction of the loading demonstrates the anisotropy.

Fig.6 Mean compressive strength for conventionally cast and 3DCP beams (loaded in the x-, y-, and z-direction) under the UC test.

Full size|PPT slide

The fractured UC test specimens were inspected to gain insight into the failure mechanism. The yellow color shade in Fig.7 shows the interface region, which indicates the onset and propagation of crack leading to failure in the 3DCP specimens. The presence of a material interface clearly results in a lower strength. The conventionally cast specimen failed in the column manner, as shown in Fig.7(a)–Fig.7(c). Although the failure pattern is column-like in all the tested cast specimens, the onset of failure cannot be pinpointed to a specific location. This is because the pores are randomly distributed in the conventionally cast specimen, and the failure nucleus is generally associated with the pore experiencing the highest stress concentration.

Fig.7 Failure images of (a) isometric view of conventionally cast specimen; (b) Face-1 of conventionally cast specimen; (c) Face-2 of conventionally cast specimen; (d) isometric view of 3DCP specimen loaded in x-direction; (e) Face-1 of 3DCP specimen loaded in x-direction; (f) Face-2 of 3DCP specimen loaded in x-direction; (g) isometric view of 3DCP specimen loaded in y-direction; (h) Face-1 of 3DCP specimen loaded in y-direction; (i) Face-2 of 3DCP specimen loaded in y-direction; (j) isometric view of 3DCP specimen loaded in z-direction; (k) Face-1 of 3DCP specimen loaded in z-direction; (l) Face-2 of 3DCP specimen loaded in z-direction.

Full size|PPT slide

In the case of 3DCP specimens, the y-direction loading has the maximum mean compressive strength, followed by the z- and x-direction loading. The compressive strength of the 3DCP specimens in the x-direction is significantly lower than that of the conventionally cast specimen. The failures are majorly confined to the interfaces in all the 3DCP specimens subject to the UC test. For specimens loaded in the x-direction, multiple depositional interfaces exist along the loading direction, which are parallel to each other, as shown in Fig.7(d). The arrangement of interfaces is shown in the bold black line. The crack propagation is in a straight line along the interface in the specimen because of the presence of the only depositional interface on the plane perpendicular to the loading direction, as shown in Fig.7(d) and Fig.7(e). The individual compressive strength and coefficient of variation (CoV) for the specimens loaded in the x-direction are shown in Table S1 in the Electronic Supplementary Material. The specimens loaded in the x-direction exhibit a minimum and maximum compressive strength of 23.2 and 28.2 MPa, respectively.

Upon compression in the x-direction, the material stretches in the y-z plane. As all depositional interfaces are parallel to the z plane, the failure originates in one of the interfaces close to the point of attaining the maximum compressive strength. Dominant cracks are present in the middle interface layers, as shown in Fig.7(e). The interfaces are generally weak compared to the parent material owing to the presence of voids, as discussed in Ref. [36].

3.1.1 Compressive behavior in the y-direction

The specimens in the y-direction have a minimum and maximum strength of 33.3 and 36.2 MPa, respectively, as presented in Table S2 in the Electronic Supplementary Material. The 3DCP specimen loaded in the y-direction have a higher compressive strength. Although multiple depositional interfaces exist along the z-direction similar to that of the specimen loaded in the x-direction, a translational interface perpendicular to the depositional interfaces is also present. Notably, crack branching is observed in the specimen after failure. This may be attributed to the presence of a translational interface between the depositional interfaces, which can arrest the free crack propagation along the loading direction, as shown in Fig.7(h). The UC tests conducted in different directions containing a combination of all types of interfaces (only translational (z-direction), only depositional (x-direction) and both (y-direction)) revealed that the cross-section perpendicular to the loading direction contained interfaces not parallel to each other to increase the failure load.

The crack lines are visible in the parent material, as shown in Fig.7(g), and along the depositional interfaces, as shown in Fig.7(h). However, no crack line along the translational interface appears. Indeed, translation interfaces of a layer are generated before depositional interfaces, and among depositional and translational interfaces, where the bonding is weaker in the latter.

3.1.2 Compressive behavior in the z-direction

In the z-direction, the minimum and maximum strengths are near 32 and 35 MPa, respectively, as presented in Table S3 in the Electronic Supplementary Material. Along the z-direction, a single translational interface exists. Moreover, no crack lines in the parent material appears, as shown in Fig.7(j). However, a dominant straight crack line, similar to the x-direction specimen, along the translational interface exist, as shown in Fig.7(k), which can be attributed to the failure of material when loaded along the z-direction.

3.2 Three-point bending behavior

In the TPB test, the notched 3D-printed specimens exhibited the highest strength compared to the conventionally cast beam, as reported by Zhang et al. [37]. The highest strength is owing to the careful selection of the appropriate material mix design and printing parameters resulting in robust layer bonding. Reinforcement techniques can further improve the strength with the same material design. In the TPB test, the same crack propagation patterns, that is, Mode-1, are observed for the casted and 3DCP beams.

Fig.8(a)–8(d) depict the fracture pattern propagating along the direction of the notch in casted and 3DCP beams. The casted and 3DCP specimens loaded in the x-, y-, and z-direction exhibit mean bending loads of 0.92, 1.06, 0.99, and 1.27 kN, respectively. The quality of the connection between layers resulting from printing parameters and curing processes is demonstrated by a greater strength in the printing specimens than the casted specimen, despite pores.

Fig.8 Crack propagation pattern in 3DCP beam: (a) conventionally cast; (b) x-direction; (c) y-direction; (d) z-direction.

Full size|PPT slide

In all cases, the crack line does not coincide with the depositional or translational interface, as shown in Fig.8. The crack propagation line of 3DCP specimens loaded in the x-, y-, and z-direction is evidence of the Mode-1 failure. Fig.9 illustrates the mean bending load concerning loading direction in the conventionally cast and 3DCP specimens.

Fig.9 Mean flexural load vs. loading direction under TPB.

Full size|PPT slide

3.2.1 Bending behavior in the x- and y-direction

The 3DCP beams in the TPB test loaded in the x- and y-direction exhibit approximately similar behavior in terms of the load-bearing capacity. The minimum and maximum load in the x-direction is 0.98 and 1.03 kN, respectively, whereas in the y-direction, they are 0.84 and 1.14 kN, respectively. Tables S4 and S5 in the Electronic Supplementary Material represent the individual load of beams under TPB and CoV, respectively. The flexural strength for the specimens loaded in x- and y-direction is lower than that in the z-direction because of the presence of interfaces along the longitudinal axis of the beam. The x-direction loading has a smaller strength that that of the y-direction loading as the cross-sectional planes of the beam are separated by the interface.

3.2.2 Bending behavior in the z-direction

The 3DCP beams subjected to bending load in the z-direction exhibit a minimum and maximum bending load of 1.23 and 1.33 kN, respectively. The beams loaded in the z-direction exhibit the maximum load compared to that of the x- and y-direction and casted beams. The maximum load in the z-direction is owing to the absence of any interfaces between the cross-sectional planes of bending. Table S6 in the Electronic Supplementary Material represents the individual load of beams under TPB and CoV, respectively. Additionally, the AE data during the TPB test was analyzed to gain insight into the salient damage and fracture mechanism.

3.3 Acoustic emission response

The AE sensor response was studied to explore the failure mechanism by observing the energy released during the bond breakage process during the TPB test. The sensor detects the maximum absolute energy released by the specimen at the time of failure in terms of voltage and time. For each specimen, the sensor detects several hits, and each hit contains many features like energy, absolute energy, frequency, and amplitude. Time is a typical parameter between the outputs of the AE sensor and universal testing machine. Based on the time set, the relationship between the load and absolute energy with respect to displacement is built, as shown in Fig.10.

Fig.10 Relationship between the mean load and absolute energy obtained from the TPB test and AE sensor, respectively, for (a) conventionally cast specimen; (b) 3DCP specimen loaded in the x-direction; (c) 3DCP specimen loaded in the y-direction; (d) 3DCP specimen loaded in the z-direction.

Full size|PPT slide

The maximum energy dissipated was recorded for the 3DCP specimen loaded in the y-direction. This was expected as multiple interfaces between the cross-sectional planes in bending exist. Each of the interfaces may act as a source of energy dissipation at the time of crack initiation. The absolute energy of the specimen loaded in the z-direction at the crack initiation is lower owing to the absence of any interfaces. Therefore, the fracture toughness is lower in the y-direction and higher in the z-direction.

4 Numerical modeling of the notched three-dimensional concrete printing specimens under three-point bending

4.1 Finite element simulation

A nonlinear numerical model was generated in ABAQUS (Simulia, Dassault Systèmes) FE package with 4-noded bilinear two-dimensional (2D) warping element (WARP2D4) to predict the mechanical behavior of notched 3DCP specimens under TPB. A 2D computer-aided design (CAD) model of the concrete beam was first imported into the software as a continuous part. A schematic for the 2D FE model is presented in Fig.11, in which an extremely thin cohesive layer is placed above the initial notch to complete the crack propagation. The model is identical to that of the experimental specimen. The out-of-plane thickness is 40 mm as the dimension of the TPB specimen is 160 mm × 40 mm × 40 mm. The bottom left support was constrained in the vertical direction only, while the bottom right support was constrained in the vertical and horizontal directions. To complete the TPB test, a load was required at the top mid-point to create tension at the notch faces for crack propagation with the combination of appropriate boundary conditions. In addition, a loading in the form of displacement of 2 mm was prescribed at the mid-point of the geometry in the negative y-direction during the simulation.

Fig.11 Schematic of the TPB specimen.

Full size|PPT slide

The directional behavior of the 3DCP specimens was simulated using both the linear and nonlinear elastic constitutive laws. The former pertains to the bulk material, whereas the latter is for the cohesive zone along the interfaces. The material properties of conventionally cast and 3DCP specimens are listed in Tab.2. The mechanical response of the 3DCP specimens depends on the loading direction.

Tab.2 Material properties

property	Young’s modulus, E (MPa)	Poisson’s ratio, $\upsilon$
conventionally cast	28254.37	0.2
3DCP specimen loaded in the x-direction	23879.44	0.2
3DCP specimen loaded in the y-direction	27873.77	0.2
3DCP specimen loaded in the z-direction	27128.13	0.2

In the next step, a mesh convergence study was performed in the x-, y-, and z-direction for the conventionally cast and 3DCP specimens at various global mesh sizes, such as 0.5, 1, 1.5, and 2 mm. One of the considerable aspects of FEM is the dependency of its result accuracy on the mesh size as the FEM with a coarse mesh is not representative of the continuous model and leads to deviations from exact results. Thus, a mesh sensitivity analysis was performed by constantly decreasing the mesh size or increasing the number of elements to reach mesh size-independent results. It is observed that the load−displacement results are similar for global mesh sizes less than 2 mm, as shown in Fig.12. Further mesh coarsening away from the crack tip can be considered, as with the adopted mesh size the simulation was completed in approximately 30–40 s. Thus, the reduction in computational time with further mesh coarsening can be neglected as the computational time is low for the existing mesh size. Therefore, a global mesh size of 2 mm was adopted to conduct the FE analysis. However, a fine mesh was used near the cohesive zone layer.

Fig.12 Mesh convergence study for: (a) casted specimen; (b) printed specimen in the x-direction; (c) printed specimen in the y-direction; (d) printed specimen in the z-direction.

Full size|PPT slide

4.2 Cohesive zone model

CZM is a popular computational tool, where the entire fracture process is grouped into the crack line and is characterized by a cohesive law that connects traction and separation over crack surfaces. Contrary to linear elastic fracture mechanics-based techniques, CZM can be used to accurately analyze crack initiation and growth. It has been effectively used to simulate fracture in materials such as concrete, metals, polymers, and functionally-graded materials [38–41]. In CZM, a thin layer, known as the cohesive layer, is considered ahead of the crack tip, as shown in Fig.13. The constitutive response of the cohesive element is defined by TSL. The crack propagates when the separation at the start of the cohesive zone exceeds the failure separation ${\delta }_{\mathrm{f}}$ or the dissipated energy is equal to the fracture energy $G_{\mathrm{F}}$ of the cohesive element. The mathematical expression for the fracture energy is expressed as

Fig.13 Schematic representation of a CZM.

Full size|PPT slide

(2)$G_{\mathrm{F}}={\int }_0^{{\delta }_{\mathrm{f}}}\sigma \left(\delta \right)\mathrm{d}\delta ,$

where $\sigma$ and $\delta$ represent the cohesive traction and separation at the crack surfaces, respectively.

The traction separation law is generally defined as [42]

(3)$\sigma ={\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}f\left(\frac{\delta }{{\delta }_{\mathrm{f}}}\right),$

where ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ represents the maximum cohesive traction, and $f$ is a dimensionless function used to define the profile of the cohesive law.

A bilinear cohesive law can efficiently simulate the fracture behavior of quasi-brittle materials, such as rock and concrete, because it highly agrees with the maximum failure load at a low-yielding factor [43]. The traction in a bilinear cohesive law linearly varies with the separation until the damage occurs. After the damage initiates, the cohesive traction decreases, as shown in Fig.14.

Fig.14 Schematic of bilinear CZM.

Full size|PPT slide

Until the initiation of the damage, the constitutive relation can be described as follows:

(4)$\sigma =K\delta ,$

where K represents the stiffness, which is a nonphysical quantity of the cohesive zone.

Large cohesive stiffness values are not recommended as they cause oscillations in cohesive tractions. Moreover, a small value creates fictitious compliance in the model [44]. Once the damage begins, the cohesive stiffness $K$ decreases, which is expressed by $\left(1-d\right)K,$ where $d$ is a damage variable and is considered to be zero till the damage initiation, that is, at ${\delta }_{\mathrm{c}}$ and 1 when it reaches the final separation ${\delta }_{\mathrm{f}}$. During the progression of damage, the instantaneous cohesive stress can be expressed as

(5)$\sigma =\left(1-d\right)K\delta ,$

(6)$d=\frac{{\delta }_{\mathrm{f}}(\delta -{\delta }_{\mathrm{c}})}{\delta ({\delta }_{\mathrm{f}}-{\delta }_{\mathrm{c}})},d\in \left[0,1\right].$

4.2.1 Determination of cohesive parameters

A triangular TSL would be ideal as the AE responses exhibit a sudden spike in the energy at the time of crack initiation, that is, at the instant of the material attaining its flexural strength. The AE sensor data are used to estimate the CZM parameters. The AE data indicate a spike in the voltage signal for a quasi-brittle material as the damage propagates. This will be helpful in deciding the form of TSLs while dealing with crack propagation modeling in different materials.

The inverse approach was used to determine the cohesive parameters, in which the cohesive parameters were iteratively chosen. Many studies have used a similar approach to obtain the cohesive parameters [44,45]. In this study, the initial value of the cohesive strength was considered as the experimental flexural strength. The cohesive stiffness is iteratively determined using a numerical trial and error procedure such that the load in terms of displacement response of the FE analysis and experiments satisfactorily agree with each other. The average iterations are less than 5 in each material direction. The critical displacement ${\delta }_{\mathrm{c}}$ is determined using Eq. (4). As the concrete specimens are brittle, the fracture energy was calculated using the critical displacement and maximum cohesive traction. These parameters were supplied to the CZM to conduct the FE analysis.

Fig.15 shows the obtained TSL parameters. The values of the cohesive stiffness, cohesive strength, and fracture energy for different loading conditions are listed in Tab.3. As expected, the fracture energy for the 3DCP specimen in the z-direction is higher. This is in line with the observations from the AE responses for the 3DCP specimens. The AE responses are comparable for the same type of specimen. A direct comparison between the AE responses of casted and 3DCP specimens is beyond the scope of this study.

Tab.3 Properties of CZM

type of specimens	Knn (MPa/mm)	nominal stress normal-only mode (MPa)	fracture energy (N/mm)
conventionally cast	20	2.10	0.21
x-direction	31	2.55	0.21
y-direction	34	2.40	0.18
z-direction	26	3.10	0.27

Fig.15 Cohesive traction in terms of separation.

Full size|PPT slide

The specimen considered in the simulation had a notch of 10 mm × 3 mm. The final separation ${\delta }_{\mathrm{f}}$ for a particular cohesive fracture energy $G_{\mathrm{F}}$ can be computed using a bilinear cohesive law as

(7)${\delta }_{\mathrm{f}}=\frac{2G_{\mathrm{F}}}{{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}}.$

4.2.2 Effect of cohesive layer thickness

In the CZM, a thin cohesive layer is considered to conduct the FE analysis. In the physical experiments, no separate cohesive layer exists. Therefore, the effect of the thickness of cohesive layer on the flexural mechanical response was studied. Four cohesive layer thicknesses of 0.01, 0.05, 0.1, and 0.5 mm were considered while all other parameters were constant to examine the influence of cohesive layer thickness. In the conventionally cast specimen, the thickness of the cohesive layer has no noticeable effect on the peak load and critical displacement nor in the elastic/post-damage region, as shown in Fig.16(a). However, the 3DCP specimens loaded in the x-, y-, and z-direction exhibit an extremely small variance in the post-damage region depending on the thickness of the cohesive layer, as shown in Fig.16(b)–Fig.16(d). The responses obtained from varied cohesive thicknesses showed that cohesive thicknesses in the range of 0.01–0.5 mm generated similar responses.

Fig.16 Effect of cohesive layer thickness on the load in terms of displacement behavior under the TPB test for different directions: (a) conventionally cast; (b) 3DCP loaded in the x-direction; (c) 3DCP loaded in the y-direction; (d) 3DCP loaded in the z-direction.

Full size|PPT slide

4.3 Finite element model validation

Fig.17(a)–Fig.17(d) show the experimental and numerical load in terms of the mid-point displacement in the TPB test for the conventionally cast and 3DCP beams. The numerical simulation using the 2D plane-strain model was used to conduct the FE analysis. The experimental results revealed that the crack propagation leading to the failure was a Mode-1 crack. In Fig.17(a)–Fig.17(d), the black lines and cloud plot show the mean bending load and standard deviation, whereas the red line shows the load in terms of the displacement response obtained from the simulation. The mean values of experimental data and numerical simulation results are comparable. The significant difference can be explained by the fact that in physical analysis, the CZ layer thickness is not there, whereas in the FE analysis, a thin layer is assumed to predict the post-damage behavior of the beam under the TPB test. The cohesive layer was placed just above the initial notch and not among the horizontal and vertical layers’ interfaces, because no shear failure among the layers existed according to the experimental results; indeed, it only fails in Mode-1. Therefore, placing the cohesive layers among all interfaces increases the computational time.

Fig.17 Load in terms of displacement for: (a) conventionally casted specimen; (b) printed specimen in the x-direction; (c) printed specimen in the y-direction; (d) printed specimen in the z-direction.

Full size|PPT slide

For conventionally cast and 3DCP beams, softening behavior is witnessed after the initiation of the damage. The mean loads and their respective mean critical displacements obtained from the experiment and numerical simulations are shown in Fig.17, whereas the deformed experimental and simulation pictures are shown in Fig.18. The experimental and numerical simulation results demonstrate a satisfactory agreement in the elastic and post-damage regions. Variations in the experimental result can be owing to flaws in the specimen preparation process, such as air cavities, which ultimately results in pores. A comparison between the experimental and numerical responses shows that the CZM within the FE framework is an excellent tool to predict the elastic and post-damage response in the 3DCP specimens.

Fig.18 Experimental and simulation cracked pictures of the beam loaded in: (a) conventionally cast; (b) x-direction; (c) y-direction; (d) z-direction specimens.

Full size|PPT slide

The effect of the direction of 3DCP beams on the performance of the cement-based beam was calculated by comparing the mean load under TPB. Experimental and numerically simulated bending loads satisfactorily agree, as shown in Fig.19(a). When the number of interfaces changes, the maximum bending load experimentally varies according to the loading direction. The bending strength of the 3DCP specimen in the y-direction is the lowest; however, it is approximately comparable to that of the specimen loaded in the x-direction. However, the specimens in the x- and z-direction reflect slightly more bending load-bearing capacity as compared to the conventionally cast specimen. The numerical simulation results for all specimens in terms of the maximum bending load sit between the highest and lowest bending strength. The bending load obtained from the simulation matched 99.88%, 97.52%, 99.48%, and 95.23% with experimentally observed mean bending load for conventionally cast specimen and 3DCP specimens in the x-, y-, and z-direction, respectively. In addition to the mean peak bending load, the elastic and post-damage behavior was obtained satisfactorily. A slight difference between the experimental and simulation results can be attributed to the fact that the beam used in the experiment could have minor flaws, such as entrapped air, which could be caused by insufficient material mixing and level of material mix in changes, resulting in pressure differences and a lack of material flow, which can lead to the formation of voids.

Fig.19 Comparison of the mean values obtained from experiments and simulations: (a) mean load; (b) mean mid-span displacement.

Full size|PPT slide

In addition, the mean mid-span deflection of conventionally cast and 3DCP beams were investigated. The experimental and numerical simulations revealed that the conventionally cast specimen reflected a higher capacity to deform compared to that of the 3DCP beams, except the beam loaded in z-direction under TPB. Fig.19(b) shows the experimental and simulation mean mid-span deflection of the beam under TPB.

The simulation results obtained for the mid-span deflection are considerably overestimated except for the 3DCP specimen loaded in z-direction, that is, 107.14%, 105.69%, 104.34%, and 95.23% with experimentally observed mean mid-span deflection for conventionally cast specimen and 3DCP specimens in the x-, y-, and z-direction, respectively.

4.4 Effects of the shape of notch

The effects of the notch shape was numerically examined using the CZM for the conventionally cast and 3DCP specimens to examine the load−displacement behavior. Fig.20 displays the notch shape varieties. Each notch had a dimension of 10 mm × 3 mm. The element developed during the simulations had a different shape as a result of adjusting the shape of notch.

Fig.20 Notch shape: (a) rectangular notch; (b) V-shape notch; (c) semi-ellipse notch.

Full size|PPT slide

Changes in the number of elements and shape of the element generated were responsible for a slight variation in the load−displacement responses for various notch shapes. The notch shape mostly determines the element’s shape and size for the entire model in the conventionally cast specimen. As for the 3DPC specimens, the effect of notch shape in the element generation is limited to the specific layer, which allows achieve equivalent results for various notch shapes. From a manufacturing perspective, a specialized cutting tool is needed to prepare the V-shape and semi-ellipse notch, whereas a regular cutting tool, such as a multi-cut machine or even a hacksaw, can be used to readily create a notch of rectangular shape. Therefore, in this study, owing to the ease of manufacturing, a rectangular shape notch was made in the TPB specimens. Nevertheless, the FE analysis revealed that the notch shape has a slight effect on the load vs. displacement response, as shown in Fig.21.

Fig.21 Load−displacement response of different notch shapes for: (a) conventionally cast specimen; (b) 3DCP specimen in the x-direction; (c) 3DCP specimen in the y-direction; (d) 3DCP specimen in the z-direction.

Full size|PPT slide

5 Conclusions

In this study, the directional properties of 3DCP specimens were investigated using both the numerical simulations and experiments. The AE analyses were carried out to explore the printed concrete fracture behavior in different loading directions under the TPB test. The FE model was able to reasonably predict the deformation patterns and failure observed in the bending experiments. The FEM results, corroborated by experiments, elucidated the role of the interface. Based on the experimental and simulation results, the following conclusions are drawn.

1) The linear elastic material responses were compared with those of the conventionally casted specimen. The compressive strength of the 3DCP specimen in the x-direction was lower compared to that of the other two principal material directions and that of the conventionally cast specimen. However, the conventionally cast and 3DCP cube specimens loaded in the y- and z-direction exhibited similar mechanical behavior in terms of the strength. The lowest strength in the x-direction during the UC test was attributed to the presence of depositional layers and porosity. Interestingly, owing to the presence of translational and depositional interfaces on the cross-section, the compressive strength in the y-direction increased significantly. The crack pattern revealed that crack initiation and propagation were confined to translational or depositional interfaces.

2) The load in terms of displacement behavior of the 3DCP specimens under the TPB test principally depended on the structure of the cross-sectional planes in bending. In addition, depositional interfaces were weaker compared to translational interfaces. Consequently, the flexural strength of the 3DCP specimen in the z-direction was higher than that of the conventionally cast specimen.

3) The AE responses from the TPB studies demonstrated a higher energy loss for cross-sectional bending planes separated by depositional interfaces, indicating damage initiation at multiple interface locations. Thus, an architectural design, including the positioning of depositional and translational layers, for a certain loading scenario would offer superior mechanical properties than those of the conventionally cast specimen. Based on the AE responses, an appropriate CZM was used within the FE framework to predict the load in terms of the displacement response. The experimental and numerical results satisfactorily agreed with each other even beyond the elastic regime.

4) The fracture toughness of the cohesive zone was inversely related to the absolute energies of the AE responses under flexion. The proposed FE model combined with the CZM could be used to simulate the behavior of 3DCP specimens satisfactorily. Parametric analysis results confirmed the insignificant effect of notch shape and cohesive layer thickness on the load-bearing capacity.

References

[1] Labonnote N, Rønnquist A, Manum B, Rüther P. Additive construction: State-of-the-art, challenges and opportunities. Automation in Construction, 2016, 72: 347–366

[2] Buswell R A, Leal de Silva W R, Jones S Z, Dirrenberger J. 3D printing using concrete extrusion: A roadmap for research. Cement and Concrete Research, 2018, 112: 37–49

[3] Menna C, Mata-Falcón J, Bos F P, Vantyghem G, Ferrara L, Asprone D, Salet T, Kaufmann W. Opportunities and challenges for structural engineering of digitally fabricated concrete. Cement and Concrete Research, 2020, 133: 106079

[4] Wu P, Wang J, Wang X. A critical review of the use of 3-D printing in the construction industry. Automation in Construction, 2016, 68: 21–31

[5] Bos F, Wolfs R, Ahmed Z, Salet T. Additive manufacturing of concrete in construction: Potentials and challenges of 3D concrete printing. Virtual and Physical Prototyping, 2016, 11(3): 209–225

[6] de Schutter G, Lesage K, Mechtcherine V, Nerella V N, Habert G, Agusti-Juan I. Vision of 3D printing with concrete—Technical, economic and environmental potentials. Cement and Concrete Research, 2018, 112: 25–36

[7] Feng P, Meng X, Chen J F, Ye L. Mechanical properties of structures 3D printed with cementitious powders. Construction & Building Materials, 2015, 93: 486–497

[8] Ma G, Li Z, Wang L, Wang F, Sanjayan J. Mechanical anisotropy of aligned fiber reinforced composite for extrusion-based 3D printing. Construction & Building Materials, 2019, 202: 770–783

[9] Panda B, Paul S C, Mohamed N A N, Tay Y W D, Tan M J. Measurement of tensile bond strength of 3D printed geopolymer mortar. Measurement, 2018, 113: 108–116

[10] Meurer M, Classen M. Mechanical properties of hardened 3D printed concretes and mortars-development of a consistent experimental characterization strategy. Materials, 2021, 14(4): 1–23

[11] Liu J, Li S, Fox K, Tran P. 3D concrete printing of bioinspired Bouligand structure: A study on impact resistance. Additive Manufacturing, 2022, 50: 102544

[12] Le T T, Austin S A, Lim S, Buswell R A, Law R, Gibb A G F, Thorpe T. Hardened properties of high-performance printing concrete. Cement and Concrete Research, 2012, 42(3): 558–566

[13] Panda B, Chandra Paul S, Jen Tan M. Anisotropic mechanical performance of 3D printed fiber reinforced sustainable construction material. Materials Letters, 2017, 209: 146–149

[14] Nerella V N, Hempel S, Mechtcherine V. Effects of layer-interface properties on mechanical performance of concrete elements produced by extrusion-based 3D-printing. Construction & Building Materials, 2019, 205: 586–601

[15] Xiao J, Liu H, Ding T. Finite element analysis on the anisotropic behavior of 3D printed concrete under compression and flexure. Additive Manufacturing, 2021, 39: 101712

[16] Paul S C, Tay Y W D, Panda B N, Tan M J. Fresh and hardened properties of 3D printable cementitious materials for building and construction. Archives of Civil and Mechanical Engineering, 2018, 18(1): 311–319

[17] Panda B, Paul S C, Hui L J, Tay Y W D, Tan M J. Additive manufacturing of geopolymer for sustainable built environment. Journal of Cleaner Production, 2017, 167: 281–288

[18] Wolfs R J M, Bos F P, Salet T A M. Hardened properties of 3D printed concrete: The influence of process parameters on interlayer adhesion. Cement and Concrete Research, 2019, 119: 132–140

[19] Lourenço P B, Rots J G, Blaauwendraad J. Continuum model for masonry: Parameter estimation and validation. Journal of Structural Engineering, 1998, 124(6): 642–652

[20] Mobasher B. M&S Highlight: Hillerborg (1985). The theoretical basis of a method to determine the fracture energy GF of concrete. Materials and Structures, 2022, 55(2): 56–57

[21] Nguyen-Van V, Li S, Liu J, Nguyen K, Tran P. Modeling of 3D concrete printing process: A perspective on material and structural simulations. Additive Manufacturing, 2023, 61: 103333

[22] Jasim W A, Tahnat Y B A, Halahla A M. Behavior of reinforced concrete deep beam with web openings strengthened with (CFRP) sheet. Structures, 2020, 26: 785–800

[23] van den Heever M, Bester F, Kruger J, van Zijl G. Mechanical characterisation for numerical simulation of extrusion-based 3D concrete printing. Journal of Building Engineering, 2021, 44: 102944

[24] Yang J, Xia J, Zhang Z, Zou Y, Wang Z, Zhou J. Experimental and numerical investigations on the mechanical behavior of reinforced concrete arches strengthened with UHPC subjected to asymmetric load. Structures, 2022, 39: 1158–1175

[25] Peng Q, Wu H, Jia P C, Ma L L, Fang Q. Numerical studies on rebar-concrete interactions of RC members under impact and explosion. Structures, 2023, 47: 63–80

[26] Liu H, Egbe K J I, Wang H, Matin Nazar A, Jiao P, Zhu R. A numerical study on 3D printed cementitious composites mixes subjected to axial compression. Materials, 2021, 14(22): 6882

[27] Bedon C, Rajcic V, Barbalic J, Perkovic N. CZM-based FE numerical study on pull-out performance of adhesive bonded-in-rod (BiR) joints for timber structures. Structures, 2022, 46: 471–491

[28] Anas S M, Alam M, Umair M. Experimental and numerical investigations on performance of reinforced concrete slabs under explosive-induced air-blast loading: A state-of-the-art review. Structures, 2021, 31: 428–461

[29] Nguyen-Van V, Liu J, Li S, Zhang G, Nguyen-Xuan H, Tran P. Modeling of 3D-printed bio-inspired Bouligand cementitious structures reinforced with steel fibers. Engineering Structures, 2023, 274: 115123

[30] Nguyen-Van V, Liu J, Peng C, Zhang G, Nguyen-Xuan H, Tran P. Dynamic responses of bioinspired plastic-reinforced cementitious beams. Cement and Concrete Composites, 2022, 133: 104682

[31] ASTMC109. Standard Test Method for Compressive Strength of Hydraulic Cement Mortars (Using 50mm Cube Specimens). West Conshohocken, PA: ASTM, 2010

[32] ASTMC880. Standard Test Method for Flexural Strength of Dimension Stone. West Conshohocken, PA: ASTM, 1999, 98–100

[33] MahmoudK S. Design of concrete Structures According to ACI 318M-14 Third Stage, 2016

[34] Gonçalves R, Giacon Júnior M, Lopes I M. Determining the concrete stiffness matrix through ultrasonic testing. Engenharia Agrícola, 2011, 31(3): 427–437

[35] MISTRASGroup Inc.—Products & Systems Division. USB-AE Node & AEwin for USB Software User’s Manual. Princeton Junction: MISTRAS Group Inc., 2010

[36] Sanjayan J G, Nematollahi B, Xia M, Marchment T. Effect of surface moisture on inter-layer strength of 3D printed concrete. Construction & Building Materials, 2018, 172: 468–475

[37] Zhang Y, Zhang Y, Yang L, Liu G, Chen Y, Yu S, Du H. Hardened properties and durability of large-scale 3D printed cement-based materials. Materials and Structures, 2021, 54(1): 45

[38] Bachir Bouiadjra B, Achour T, Berrahou M, Ouinas D, Feaugas X. Numerical estimation of the mass gain between double symmetric and single bonded composite repairs in aircraft structures. Materials & Design, 2010, 31(6): 3073–3077

[39] Odi R A, Friend C M. An improved 2D model for bonded composite joints. International Journal of Adhesion and Adhesives, 2004, 24(5): 389–405

[40] Feih S, Shercliff H R. Adhesive and composite failure prediction of single-L joint structures under tensile loading. International Journal of Adhesion and Adhesives, 2005, 25(1): 47–59

[41] Magalhães A G, de Moura M F S F, Gonçalves J P M. Evaluation of stress concentration effects in single-lap bonded joints of laminate composite materials. International Journal of Adhesion and Adhesives, 2005, 25(4): 313–319

[42] Jin Z H, Sun C T. A comparison of cohesive zone modeling and classical fracture mechanics based on near tip stress field. International Journal of Solids and Structures, 2006, 43(5): 1047–1060

[43] Schellekens J C J, de Borst R. A non-linear finite element approach for the analysis of mode-I free edge delamination in composites. International Journal of Solids and Structures, 1993, 30(9): 1239–1253

[44] de Moura M F S F, Gonçalves J P M, Magalhães A G. A straightforward method to obtain the cohesive laws of bonded joints under mode I loading. International Journal of Adhesion and Adhesives, 2012, 39: 54–59

[45] deMoura M F S FCampilhoR D S GGonçalvesJ P M. Crack equivalent concept applied to the fracture characterization of bonded joints under pure mode I loading. Composites Science and Technology, 2008, 68(10−11): 2224−2230

Acknowledgements

This work was supported by the Science and Engineering Research Board, India, under the scheme “Early Career Research Award (No. ECR/2018/001638)”, DST-SERB and VSSC, and ISRO of the project titled “Functionality Enhancement Through Design and Development of Advanced Finite Element Algorithms for STR Tools” under the IMPRINT.IIC (IMP/2019/000276) scheme. The authors acknowledge the Smart Materials and Structure Laboratory for facilitating the AE sensor test and DST-FIST funded facility, SHIMADZU AGS-X 100-kN universal testing machine, at the Department of Mechanical Engineering, Indian Institute of Technology Guwahati. Biranchi Panda would like to thank the Science and Engineering Research Board (SERB), India, for the start-up grant (No. SRG/2021/000052).

Electronic Supplementary Material

Supplementary material is available in the online version of this article at https://doi.org/10.1007/s11709-023-0004-z and is accessible for authorized users.

Conflict of Interest

The authors declare that they have no conflict of interest.

RIGHTS & PERMISSIONS

2023 Higher Education Press