Computational modeling of fiber orientation during 3D-concrete-printing

During 3D-printing of fiber-reinforced concrete, fibers tend to align with the printing direction due to strong shearing deformation of the material, allowing for the controlled production of components with desired fiber orientation states. The accurate prediction of the fiber orientation state in printed components poses a major challenge due to the large number of processing and material parameters involved and due to the complex mechanisms of flow and fiber reorientation during printing. This contribution presents a novel incorporation of the Folgar–Tucker fiber orientation model within a fluid dynamics framework based on the Particle Finite Element Method for simulations of the fiber orientation evolution during 3D-concrete-printing. The fiber orientation state is represented using a second-order orientation tensor, which is coupled with a new anisotropic Bingham constitutive model used for the viscous fiber-concrete mixture to account for the effect of fiber orientation on the velocity field. Further, the orientation distribution function is reconstructed from the second-order orientation tensor, following the maximum entropy method for a more convenient interpretation of the results. The model is validated by comparing the simulated orientation numbers of a 3D-printed concrete layer for different extrusion nozzle diameters with experimental values from the literature. Several parametric studies are performed to examine the flow and fiber reorientation mechanisms and the influence of process parameters on the fiber orientation state in printed components. Stronger fiber alignment in the printing direction is obtained for higher printing speeds or smaller extrusion nozzles, associated with higher shear stresses developing in the extrusion nozzle.


Introduction
Automated construction techniques, such as additive manufacturing with fresh concrete, have become a rapidly evolving field of research and application in recent years. In particular, extrusion-based approaches, like 3D-concrete-printing, are very popular among technologies based on the additive manufacturing using fresh concrete. 3D-concrete-printing is characterized by an automated layer-wise deposition of a "stiff" fresh concrete material through an extrusion nozzle to produce structures without formwork. Some state-of-theart reports about additive manufacturing using fresh concrete are given in [1][2][3][4]. To make 3D-concrete-printing a more reli-B Günther Meschke guenther.meschke@rub.de 1 Institute for Structural Mechanics, Ruhr-University Bochum, Bochum, Germany able and efficient production technology, several challenges, as discussed in the state-of-the-art reports above, must be overcome, among which the incorporation of reinforcement is one of the biggest to be solved [5]. The utilization of fiberreinforced concrete as the printing material is one strategy that addresses the lack of reinforcement in printed components. Several experimental studies revealed that fibers tend to align in printing direction during 3D-concrete-printing [5][6][7][8], which is a substantial advantage compared to traditionally cast fiber-reinforced concrete components. The studies also indicated that the degree of fiber orientation could be controlled by adjusting process parameters, i.e., flow rate or printing speed, during printing [8]. However, the prediction of the fiber orientation in printed components must be well understood to allow for a reliable printing process. Due to the multitude of process and material parameters involved during extrusion and deposition, in combination with the complex fiber reorientation and flow mechanisms, such predictions are highly complex and require extensive experimental and numerical investigations. In this context, by introducing a numerical approach based on the Particle Finite Element Method (PFEM) [9,10] and a fiber orientation model according to Folgar and Tucker [11], the extrusion process of fiber-reinforced 3D-concrete-printing is analyzed in this paper.
Numerical modeling strategies for fiber-reinforced fresh concrete are based on either discrete approaches or on models built on Jeffery's orbit [12], like the probabilistic approach according to Folgar and Tucker [11]. The Discrete Element Method (DEM) is one candidate for numerical modeling of fiber-reinforced fresh concrete using the discrete approach, as proposed in [13], in the framework of simulations of the slump flow test. Smoothed Particle Hydrodynamics (SPH) was also used for simulations of the slump flow test using fiber-reinforced fresh concrete based on the discrete approach in [14] and for simulations of the L-box test in [15]. A well-calibrated discrete fiber orientation model should generally provide the most accurate fiber orientation predictions due to the possibility to account for various interactions between fibers, aggregates, and boundaries. However, the high computational cost and the question of the scale on which the heterogeneity of fresh concrete suspension is modeled restrict applications of this approach to small scale simulations, e.g. small laboratory tests. On the other hand, the fiber orientation state can be modeled based on a probabilistic approach, representing the fiber orientation by second-order tensors, according to Folgar and Tucker [11]. These models have the advantage of being computationally more efficient than the discrete approach at the cost of accuracy, i.e., direct contacts between fibers, aggregates, and boundaries are not accounted for directly. In addition, when using an orientation tensor to represent the fiber orientation state, information of the actual fiber distribution is lost and a reconstruction of the orientation distribution is required for a better comparison with experiments. An SPH approach for fiber-reinforced fresh concrete based on the Folgar-Tucker model was developed in [16]. This model revealed that the simulated fiber orientation states predicted with the model qualitatively agreed with the experiments of an L-box test. Lattice-Boltzmann-based formulations for flow simulations of fiber-reinforced fresh concrete were given in [17,18], where the casting processes of self-compacting concrete into formwork were discussed. Computational Fluid Dynamics simulations were used to predict the flow-induced fiber orientation of ultra high performance fiber reinforced concrete elements in [19]. Apart from fiber-reinforced fresh concrete, numerical modeling using the Folgar-Tucker model mainly focuses on the application using polymers [20,21]. For example, in respect to additive manufacturing, numerical modeling strategies for polymers in Fused Deposition Modeling (FDM) were developed by using a finite element model in [22] and an SPH approach in [23] based on a Folgar-Tucker model. These simulations were conducted in 2D assuming a planar flow, not considering the additional complexity due to the flow in lateral direction relevant in 3D.
Despite the numerical studies on related topics listed above, the extrusion process of fiber-reinforced 3D-concreteprinting and the correlations between process parameters and the resulting fiber orientation state in the printed filaments, to the best knowledge of the authors, have not been analyzed numerically so far. While in general any of the numerical methods discussed above can be used for such simulations, this work builds upon a previously developed PFEM model [24][25][26], which has been successfully applied to extrusion simulations during 3D-concrete-printing. To circumvent the high computational cost associated with a discrete fiber modeling approach, in the proposed model, the orientation state of fibers is calculated using a Folgar-Tucker fiber orientation model that also builds upon a previous implementation in the context of fresh concrete flow using SPH [16]. In this model, the fiber orientation state is represented by a secondorder orientation tensor according to Jeffery [12] and Folgar and Tucker [11,27]. The fiber orientation model is summarized in Sect. 2. Section 3 gives an overview of the adopted anisotropic Bingham model and the upscaling formulas that relate the rheological parameters with the fiber orientation state. An overview of the PFEM and an efficient solution algorithm for the fiber orientation model within a PFEM framework is covered in Sect. 4. This section also introduces a methodology for reconstructing the orientation distribution function based on the maximum entropy method [28], which is used as a post-processing tool in the numerical analyses to quantify fiber orientation states, allowing for a more convenient interpretation of the results. The current work improves the previous SPH model given in [16] in terms of the new anisotropic Bingham constitutive model and the post-processing tool based on the maximum entropy method. In addition, this work presents an extension of the previously developed PFEM model (see [24][25][26]) in view of the implementation of the fiber orientation model, which has also not been implemented in a PFEM framework in general. Finally, numerical analyses of the extrusion process during 3D-printing of fiber-reinforced concrete are performed in Sect. 5. The numerical results are compared with experimental results from the literature, and correlations between process parameters and the resulting fiber orientation states in printed layers are analyzed in parametric studies. The presented results may help tuning the process parameters during the 3D-printing of fiber-reinforced concrete to obtain desired and optimized properties of printed structures concerning the fiber orientation state.

Folgar-Tucker fiber orientation model
Fiber suspensions can be categorized into different concentration regimes based on the fiber volume fraction c and the aspect ratio r of a single fiber [29]. The fiber aspect ratio r is defined as the ratio of the fiber length L f and the fiber diameter D f as r = L f /D f ; see Fig. 1. In the dilute regime (c < 1/r 2 ), the motion of fibers is solely influenced by drag forces from the liquid matrix. As the number of fibers increases, the semi-concentrated regime (1/r 2 < c < 1/r ) is reached, and hydrodynamic interactions between fibers are no longer negligible. In the highly concentrated regime (c > 1/r ), the motion of fibers is dominated by direct contacts between neighboring fibers. In this regime, fiber clumps and contact networks between fibers may be formed. Based on these definitions, much work has been invested to develop theoretical frameworks for fiber suspension flow models, mainly in the context of polymers, assuming a perfectly homogeneous Newtonian liquid matrix (see, for instance, [20,21]). Only few publications are focusing on developing theoretical frameworks for fiber suspension flow models assuming a viscoplastic [30] or fresh concrete matrix [31][32][33]. In the case of fresh concrete, the fiber volume fraction is mostly sufficient to allow for hydrodynamic interactions and direct contacts between the fibers [34], which are the semi-concentrated and highly concentrated regimes. The cases studied in this work are in the range of the semi-concentrated regime.
This Section introduces the evolution equation for the fiber orientation state by using a well-known approach by Folgar and Tucker [11,27]. According to this approach, a secondorder orientation tensor represents the fiber orientation state of a group of fibers, and its evolution equation is solved using a closure approximation of the fourth-order orientation tensor.

Evolution equation for the fiber orientation state
The orientation state of a single fiber in a Cartesian frame is determined by the unit vector p pointing in the fiber direction, which is defined using spherical coordinates as with the azimuth angle φ and the polar angle θ , as illustrated in Fig. 1. To define arbitrary orientation states of a larger population of fibers, the orientation distribution function ψ(θ, φ) is introduced as a continuous function over the surface of a unit sphere [11]. The orientation distribution function must be periodic (ψ(θ, φ) = ψ(π − θ, φ + π)) and must satisfy the normalization condition, given as the integral of the orientation distribution function over the surface of a unit sphere S as While the evolution of a single fiber suspended in a Newtonian fluid was given by Jeffery [12], Folgar and Tucker combined Jefferey's model with a diffusivity term for fiber interactions and recast it into a probabilistic representation [11]. The proposed continuity equation of the orientation distribution function ψ(θ, φ) by Folgar and Tucker is given as [27] with the gradient operator in orientation space (surface of a unit sphere) ∇ S , the strain rate tensor , which is close to unity for the aspect ratios studied in this work, and the equivalent strain rate |γ | = 2D i j D i j . The hydrodynamic contributionṗ h in Eq. (3) is the wellknown evolution law for the motion of a single spheroidal particle suspended in a Newtonian fluid given by Jeffery [12]. According to [11], q h is a phenomenological diffusivity term, which accounts for the interactions between fibers. Since these interactions are random, their influence on the fiber orientation is considered diffusive. The existence of the diffusivity term results in a steady state fiber orientation state obtained during flow, which was also observed in case of fresh concrete [31,34]. In this work, the isotropic rotary diffusivity factor C I is obtained from the relation C I = 0.03(1.0 − e −0.224cr ) suggested in [35]. As seen in this relation, the diffusivity factor increases with the fiber volume fraction c due to the higher probability of interactions between the fibers. Note that in case of fresh concrete, aggregates can also contribute to the diffusivity term due to additional interactions. However, to the best knowledge of the authors, no experimentally validated diffusivity model to account for the additional interactions due to aggregates exists so far.

Fiber orientation modeling using orientation tensors
The direct solution of the continuity equation of the orientation distribution function Eq. (3) is computationally expensive. Instead, it was found that the evaluation of moments of the distribution function is computationally far more efficient and sufficiently accurate in most applications [27]. Such moments of the distribution function are called orientation tensors of order n and can be interpreted as the two-dimensional counterparts to moments of a linear distribution function, such as the expected value (first moment) and variance (second moment). While the odd-ordered moments of a symmetric orientation distribution function defined on the unit sphere are zero, only the even moments can be evaluated. In typical applications, it is sufficient to consider the second and fourth order orientation tensors, which are defined as The evolution equation of the second order orientation tensor is obtained from the time derivative of Eq. (4) by using Eq.
(3) and integration by parts, see [27], leading tȯ Evolution equations for higher-order orientation tensors may also be derived similarly. Note that the evolution equation of an orientation tensor always requires a higher order orientation tensor. For example, the evolution equation for the second-order orientation tensor a 2 in Eq. (6) includes the fourth-order orientation tensor a 4 . Such problems are called closure problems, and as a consequence, the solution of the evolution equation of the second-order orientation tensor (Eq. (6)) requires an approximation of the fourth order orientation tensor. The fourth order orientation tensor can be approximated using a certain function of the second order orientation tensor, which is the so-called closure approximation. Numerous closure approximations for orientation tensors of different order have been proposed in the literature, each with its individual properties and accuracy. Among the possible closure approximations, the orthotropic fitted closure (ORW3) was chosen in this work to approximate the fourth-order orientation tensor by the second-order orientation tensor. It offers very good accuracy and a straightforward calculation procedure (see [36] for more details). To slow down the evolution rate of the fibers, Wang et al. [37] introduced the reduced strain closure (RSC) to yield more realistic orientation kinetics in a fiber suspension. For clarity, it should be noted that the RSC is not a closure approximation for the fourth-order orientation tensor like the ORW3 closure and can be used together with the ORW3 closure. The RSC only slows down the evolution rate of the fibers by using an empirical factor κ and certain tensors based on the eigenvalues and -vectors of the second-order orientation tensor a 2 . Equation (6) is modified by using the RSC aṡ with L i jkl = The strain rate tensor D, the spin tensor W and the equivalent strain rate |γ | appearing in Eq. (7) must be obtained from a velocity field. In this work, the velocity field is calculated using the Particle Finite Element Method (PFEM). A numerical solution procedure of Eq. (7) within the fluid dynamics framework of PFEM is discussed in Sect. 4.

Constitutive model for fresh concrete flow including fibers
This section outlines a phenomenological constitutive model for the flow of fiber-reinforced fresh concrete. An extended Bingham model is proposed to account for material anisotropy due to fiber orientation. Finally, upscaling formulas for the rheological properties of fresh concrete mixes, including fibers, are presented.

A regularized Bingham model for fresh concrete without fibers
A fresh concrete mix that consists of different aggregates, water, and cement, acts as a yield stress fluid on the macroscopic scale [34]. The material response of homogenized fresh concrete can be simplified by neglecting shear-thinning effects and using a Bingham model as a constitutive assumption. Numerically preferred is the regularized version of the where I d i jkl = (δ ik δ jl + δ il δ k j )/2 − 1 3 δ i j δ kl , δ i j denotes the Kronecker delta, μ denotes the plastic viscosity,μ = μ + τ 0 |γ d | 1 − e −m|γ d | denotes the apparent viscosity, τ 0 denotes the yield stress, m denotes a regularization parameter, and |γ d | = 2D d i j D d i j denotes the equivalent deviatoric strain rate, with D d i j = I d i jkl D kl . As seen in Eq. (8), a nearly rigid material response of stress states below the yield stress is achieved by penalizing the deformations with a large value of the apparent viscosity, controlled by the regularization parameter m. A regularization parameter tending to infinity would, in principle, yield an ideal Bingham model. However, the conditioning of the solution algorithm suffers from overly large values of m and a good trade-off between numerical stability and the error in the approximation of the Bingham model must be found. From experience and literature [39], typical values for m are between 100 and 3000. In this work, a value of 1000 is used.

An anisotropic constitutive model for fiber-reinforced fresh concrete suspensions
When fibers are added to a fresh concrete mix, the probability of solid particle collisions is increased, resulting in a larger viscosity and yield stress of the mixture. In addition, the suspension is characterized by an anisotropic response due to fiber reorientation during flow [32,40], which manifests itself through increased elongational viscosities in fiber direction and shear-thinning effects [20]. Sommer et al. [41,42] proposed a constitutive model for anisotropic fiber suspensions. In this approach, a transversely anisotropic fluid model is volume averaged according to the approach described in [27] using the second and fourth order orientation tensors. According to [41,42], assuming the same values for the in-plane and transverse shear viscosities, and utilizing the regularized Bingham model (Eq. (8)), the deviatoric Cauchy stress tensor of the fiber suspension can be defined as where μ i jkl is the fourth order viscosity tensor, η 23 is the transverse shear viscosity and R η is an anisotropy ratio, which becomes one in case of the absence of fibers. To the best knowledge of the authors, no constitutive assumption for R η is available for the case of fiber-reinforced fresh concrete. Therefore, the anisotropy ratio is assumed to take the form based on [42,43], which was derived for a highly concentrated Newtonian suspension as R η = 1 + 1 8 c √ 12c/πr 2 . As can be observed, R η can be much larger than unity for large aspect ratios, which would increase the anisotropic behavior of the material. The transverse shear viscosity must be equal to the viscosity μ in absence of fibers (c = 0). In this way, the consistency is maintained because Eq. (9) reduces to Eq. (8) when the fiber volume fraction is null. Specific constitutive assumptions for the transverse shear viscosity η 23 and the yield stress τ 0 of the fiber suspension are discussed in the following section. Similar to the viscosity, the yield stress of fiber-reinforced fresh concrete can also show an anisotropic response depending on the state of fiber orientation, as observed in [30,32]. However, due to the lack of experimental data and theoretical modeling strategies, this effect is neglected in this work. Due to major and minor symmetries, the fourth order viscosity tensor μ can be represented in Voigt notation and is used as the constitutive tangent in the solution procedure.

Rheological properties of a fresh concrete fiber suspension
The rheological properties of fiber suspensions are well understood in the case of a Newtonian fluid matrix. In a fiber suspension, hydrodynamic interactions contribute to the macroscopic rheological properties. As the number of solid particle collisions increases, the viscosity of the suspension will also increase. Significant work has been done on this topic in the case of polymer composites (see, e.g., [20,21]). Comparatively, less attention was given to the development of upscaling relations for rheological properties of Bingham fluid fiber suspensions. Férec et al. [30] developed a model for the semi-concentrated regime of Bingham rigid fiber suspensions, accounting for hydrodynamic interactions between fibers. Furthermore, Férec et al. [44] derived a model for semi-concentrated rod suspensions in a Bingham fluid using a cell model approach. These models predict similar correlations between the viscosity and the fiber content in Bingham fluid compared with a purely Newtonian fluid matrix. Despite these pioneering efforts, the model for fiber-reinforced Bingham suspensions developed by Férec et al. is not practical to use together with the fiber orientation modeling framework introduced in Sect. 2. Also, this framework does not account for the aggregate-fiber interactions, which play an essential role in fresh concrete rheology, especially when the maximum packing fraction of the solid particles is reached [31][32][33]45]. Therefore, a phenomenological approach for fresh concrete fiber suspensions is adopted. Ghanbari and Karihaloo [46] derived and validated a micromechanically motivated viscosity model for isotropic randomly oriented rigid fibers suspended in a viscous cementitious phase. According to this approach, the viscosity of a fresh concrete mix including fibers can be approximated as whereμ is the viscosity of the fluid matrix (fresh cement paste or concrete). The anisotropic model can be combined with the viscosity for isotropic randomly oriented rigid fibers in Eq. (10) by enforcing Eq. (9) to yield the same viscosity as obtained in Eq. (10) for the case of a random isotropic fiber orientation. The orientation tensors are well defined for a random isotropic orientation state using ψ = 1 4π . This yields an approximation for the transverse shear viscosity as η 23 = 15 4R η +11 μ. Several experimental studies of the slump flow test [14,31,40] have shown that the yield stress also increases with rising fiber content due to a higher probability of solid particle collisions. A fundamental investigation of the influence of rigid fibers on the yield stress of cementitious materials was reported in [32]. Based on arguments obtained by applying the contact network models for spherical inclusions, it was demonstrated and experimentally validated that the yield stress increased by several magnitudes when the maximum solid packing fraction was approached. In this context, a physical model to predict the influence of flexible fibers on the yield stress of cement pastes and mortars was developed and validated in [45]. Based on [32,45] the yield stress of a fiber-reinforced fresh concrete suspension for isotropic random fiber orientation states is given as whereτ 0 is the yield stress of the cement paste, φ f m = 4/r is the dense packing fraction of fibers, φ s is the volume fraction of solid particles (sand and aggregates), and φ sm ≈ 0.65 is the maximum packing fraction of solid particles. Note that Eq. (10) was developed for fiber volume fractions of up to c = 2 %, while Eq. (11) may also account for higher fiber concentrations. To avoid inconsistencies between viscosity and yield stress scaling when using higher fiber volume fractions, we restrict the analyses in this work to fiber volume fractions lower than c < 2 %. It should be pointed out that Eq. (10) only takes into account the increase in yield stress due to fibers with an isotropic random orientation state, but not the anisotropic effects resulting from reorientation of the fibers [30,32]. As discussed in Sect. 3.1, the effect of fiber reorientation on the yield stress is neglected in this work due to the lack of experimental data and theoretical modeling strategies.

Numerical solution procedure
This Section provides an overview of the Particle Finite Element Method (PFEM) for fluid dynamic problems. The implementation is based on a standard PFEM formulation and, therefore, is presented here in a condensed form. Further information can be found in the references provided in the individual subsections. In addition, a numerical solution strategy for the orientation tensor within the PFEM framework is proposed. The final part of this section is attributed to introducing a post-processing tool based on the maximum entropy method for reconstructing the orientation distribution function from the second-order orientation tensor [28].

Balance of momentum and mass
Fresh concrete is a heterogeneous material consisting of multiple components, like cement clinker, water, coarse and fine aggregates, and possibly different admixtures or fibers. Numerical modeling approaches, which explicitly account for the individual components, can only be realized at extremely high computational costs. As a result, fresh concrete is often assumed to act as a homogeneous mixture in numerical simulations. In this context, the balance of momentum in the time interval [0, T] of a domain Ω is given as with the density of the fresh concrete mixture ρ, the Cauchy stresses σ , the body forces ρb, and the material time derivative of the velocityv. The velocitiesṽ i are prescribed over the Dirichlet boundary d of Ω as v i =ṽ i and the tractionst i are defined on the Neumann boundary n of Ω as σ i j n j =t i . Assuming a slight compressibility of the fresh, flowable concrete and a split of the Cauchy stresses into deviatoric and hydrostatic components, the balance of mass is given in the time interval [0,T] over the domain Ω as where K denotes the bulk modulus andṗ denotes the material time derivative of the pressure. The Cauchy stresses are defined as σ i j = τ i j + δ i j p. For the constitutive assumption of the deviatoric stresses of fiber-reinforced fresh concrete we refer to Sect. 3. Note that the material time derivatives in Eqs. (12) and (13) are defined in a Lagrangian framework aṡ = ∂ ∂t .

Numerical discretization
In PFEM, the fundamental equations over 2D/3D domains are discretized using standard low order simplex finite elements. The discretized residual form of Eqs. (12) and (13) are obtained following the usual procedure, by casting the balance equations and the corresponding Neumann boundary conditions into the weak form and approximating the velocities, pressures and their virtual counterparts using shape functions as with the nodal velocitiesv, the nodal pressuresp, the mass matrix M vv , the internal force vector F v,int , the external force vector F v,ext , the velocity pressure coupling matrix G vp , the pressure stiffness matrix K pp and the stabilization matrix S pp . The stabilization matrix is used to treat the LBB criterion, which is violated for equal order finite element approximations of the velocity and pressure fields and is obtained from a direct pressure stabilization using the polynomial pressure projection technique [47]. The individual matrices and vectors of Eqs. (14) and (15) are assembled from the contributions of each element e as with the finite element assembly operator A , the number of elements NE, the shape functions of the velocity and pressure fields N v and N p , the elemental volume Ω e , the Cauchy stresses in Voigt notation σ e , the stabilization parameter α stab chosen as 0.5 in this work, the fluid viscosity μ, the discontinuous polynomialsN p [47], and m = [1, 1, 1, 0, 0, 0] T . The discrete system of equations is obtained by assembling the individual element matrices. Finally, the equations are solved using a monolithic solution procedure and an implicit backward Euler time discretization scheme for the velocity and pressure fields. In this context, the viscosity matrix is introduced as with the fourth-order viscosity tensor in Voigt notation μ e defined in Eq. (9).

The particle finite element method
The Particle Finite Element Method (PFEM) is a re-meshingbased Finite Element Method in an updated Lagrangian description [9,10]. As outlined in Sect. 4.1.2, the domain is discretized using standard low order simplex elements in 2D and 3D. The original re-meshing algorithm in PFEM was based on a Delaunay tessellation of the particles (nodes) and a free surface reconstruction algorithm based on the alpha shape method [48]. In contrast, the re-meshing algorithm used in this work was improved by using a constrained Delaunay tessellation to ensure that the free surface is modified as little as possible [26,49]. In this case, the free surface discretization is used as a constraint for the tessellation. The initial free surface is obtained from the alpha shape method and is updated by insertion and deletion of particles (nodes), followed by a local re-tessellation when the free surface discretization becomes too distorted, i.e., when surface facets become too small or nodes are too close to each other. Solution variables are then linearly mapped from the previous discretization to the new nodes. Adopting the constrained Delaunay tessellation improves the free surface recognition and volume conservation substantially due to minimal changes of the free surface. On the other hand, the procedure is limited to smooth surfaces and cannot account for excessive motions, e.g., wave splashing. Besides the constrained tessellation, mesh quality improvement is used to avoid poorly shaped elements and sliver elements [26]. In this paper, contact with the ground is modeled using no-slip conditions. Details of different contact formulations in PFEM can be found in [49,50], and more recent PFEM developments are summarized in [51]. The basic steps of the adopted PFEM algorithm and the new calculations associated with the fiber orientation model are summarized as 1. Initialize the simulation by defining t 0 and setting the time step counter to n = 0. Define the initial cloud of points and obtain the free surface discretization of the problem by the alpha shape method. 2. Start the new time step by setting t n+1 = t n + Δt. 3. Update the free surface discretization for the current time step by insertion or deletion of particles when necessary, and perform a local re-tessellation in such cases. 4. Perform the constrained Delaunay tessellation of the domain to obtain the finite element discretization. Map solution parameters from the old discretization to the new nodes. 5. Iteratively solve the system of equations, Eqs. (14) and (15), in the updated configuration at t n+1 to obtain the new velocities, pressures, and nodal coordinates for the converged solution. Within each iteration the new stresses are calculated based on the proposed constitutive model for fresh concrete including fibers according to Eq. (9). The orientation tensor is assumed constant within this calculation step and is taken from the previous time step. 6. Calculate the new fiber orientation state using Eqs. (23) and (25). The velocity field remains constant within this calculation step and is taken from the previous step 5. 7. Set n = n + 1 and go back to step 2 and repeat all steps for the new time step as long as the maximum time t max is not yet reached.
Step 6, which covers the fiber orientation state calculations, is described in more detail in the following section. Note that within the solution algorithm above, the balance of momentum and mass are explicitly coupled with the fiber orientation model. Solution of the balance of momentum and mass (Step 5) is obtained using the orientation tensor from the previous time step and the new orientation tensor is found using the updated velocities (Step 6). This leads to a robust and efficient solution scheme. In terms of computational efficiency, the number of iterations to reach convergence in Step 5 using the anisotropic Bingham model is not significantly different from that when using a standard Bingham model.
Step 6 represents an additional calculation step compared to a standard PFEM model. However, as explained in the next section, calculation of the second order orientation tensor can be performed locally for each node and does not require an expensive solution procedure. The additional computational resources associated with the fiber orientation model can be considered small, which is a considerable advantage of Folgar-Tucker fiber orientation models as compared to discrete fiber orientation models. All implementations and simulations are done based on a Matlab version of PFEM.

Algorithmic formulation of the Folgar-Tucker model
This subsection provides an efficient implementation and solution procedure of the Folgar-Tucker fiber orientation model given in Eq. (7) within a PFEM framework. When the reference frame of Eq. a 2 is chosen equal to the Lagrangian frame of the PFEM model, the time derivative of a 2 is given asȧ 2 = ∂ a 2 /∂t. Therefore, no convective terms appear in the time derivative of a 2 , and consequently, no spatial derivatives of a 2 appear in Eq. (7). Due to the absence of convective terms, Eq. (7) can be easily integrated directly. Defining discrete solution values of a 2 on integration points where the elemental velocity gradients are calculated might be troublesome during re-meshing in PFEM, when the connectivities change and integration points are lost. As an alternative, a 2 is defined on nodes using a nodebased smoothed velocity gradient to solve Eq. (7), which is inspired by Smoothed Finite Element formulations [52]. In the case of a finite element discretization based on linear simplex elements, the node-based smoothed velocity gradient reduces to the volume average of all elemental velocity gradients adjacent to a node. Note that the smoothed velocity gradient is only used to solve the Folgar-Tucker equation, and the PFEM formulation is still based on a standard finite element discretization. Further details concerning node-based Smoothed Finite Element formulations are found in [52].
As explained in the previous Sect. 4.1.3, the balance of momentum and mass and the fiber orientation model are coupled explicitly. Thus, the new fiber orientation state, based on the second order orientation tensor, is calculated using the already updated velocities, which are constant during this step. The second order orientation tensor is explicitly integrated in time by splitting the time step n into m max smaller time steps. Using Eq. (7), the updated second order orientation tensor is given for the m th sub-step as where Δt is the time step size. The superimposed bar denotes nodal quantities and the superimposed tilde are quantities obtained from the smoothed nodal velocity gradient. The deformation rate and spin tensors are defined for the substepping scheme as Equation (23) is successively solved m max times in one time step. For the numerical problems presented in this paper, m max = 1000 was found to yield a sufficient accuracy. The implementation is based upon a description in Mandel notation as outlined in [16].

Reconstruction of the orientation distribution function with the maximum entropy method
As detailed out in the previous sections, the fiber orientation state is represented solely based on the second order orientation tensor. However, the interpretation of this tensor is less intuitive and practical than that of the orientation distribution function. The orientation distribution function can be reconstructed from the orientation tensors to provide a more conclusive analysis of fiber orientation states and a more consistent comparison with experiments. This can be accomplished by truncating the spherical harmonic series as described in [27], which requires multiple higher-order orientation tensors to yield satisfying results. Alternatively, the maximum entropy method can be used to find the orientation distribution function that fits the problem, as described in [28]. The entropy S of the fiber orientation is defined as The maximum entropy method is applied by maximizing Eq.
The distribution coefficients d i cannot be calculated directly and must be obtained through an iterative procedure by min-imizing the difference between the second order tensor a 2 and the second order orientation tensor calculated using the Bingham distribution, Eq. (27). Based on the definition of the second order orientation tensor a 2 in Eq. (4), the distribution coefficients d i can be obtained from minimization of the objective function f defined in principal direction of a 2 as where λ i is the i th eigenvalue of a 2 . Note that a 2 is still obtained from the numerical simulation based on Eq. (23). The second part of Eq. (29) is the eigenvalue of the second order orientation tensor evaluated using the Bingham distribution (Eq. (27)). The minimization problem is numerically solved by an unconstrained minimization solver provided by Matlab and a discretization of the unit sphere surface as described in [28].
It should be pointed out that the reconstruction of the orientation distribution function using the Bingham distribution in Eq. (27) is done only during post-processing after the simulation has finished. The reconstruction of the orientation distribution function allows a more convenient interpretation of the results and comparison with experiments. For the interpretation of the fiber orientation states in the following sections, the orientation distribution function ψ plane on the x-y-plane and the orientation number N o are introduced as The orientation number is a measure of the fiber alignment in x-direction, which is the printing direction in the following sections. The orientation number approaches zero when the number of fibers pointing in the x-direction is small, and becomes one when all fibers are oriented in the xdirection. Accordingly, a higher orientation number means a stronger fiber alignment in x-direction. For symmetric problems w.r.t. the x-z-plane and when the sign of φ is irrelevant, the orientation distribution function on the x-y-plane ψ plane is re-defined as with φ ranging between 0 and π/2.

Numerical analysis of extrusion processes in 3D-printing of fiber-reinforced concrete
In this section, the fiber orientation state during the extrusion and deposition process in 3D-printing of fresh fiberreinforced concrete is simulated and analyzed by means of the numerical model presented in the previous sections. The numerical model is validated based on a comparison of the simulated orientation numbers (Eq. (31)) of a 3D-printed concrete layer for different extrusion nozzle diameters with experimental values from the literature [8]. In this context, the characteristics and limitations of the numerical implementation are discussed. Further parametric studies were discussed to delineate the influence of process parameters such as printing speed, printing height, and flow rate on the fiber orientation state in a single printed layer. The computational investigations aim to answer the question whether and to which extent the fiber orientation state in a printed layer can be controlled with the different process parameters to allow a more precise control of the printing process. In addition, the computational results are intended to provide novel insights into the flow mechanisms and fiber reorientation during printing. The printing process is characterized by printing one layer of fresh fiber-reinforced concrete using a circular extrusion nozzle similar to the setup presented in [8,26,54].

Description of the numerical analyses
Only a part of the extrusion nozzle was modeled in the numerical simulation. The simulation of fresh concrete flow and fiber orientation through the extruder was neglected for the current investigations and only straight round extrusion nozzles were considered to simplify the simulation and to yield a more generic analysis of different nozzle geometries. Hence, it is assumed that the fiber orientation state reaches a steadystate by the time it reaches the exit of the nozzle. This was taken into account by setting the initial condition for the fiber orientation at the inlet equal to the steady-state solution of fiber suspension flow in the pipe. The model was further simplified utilizing symmetry conditions along the printing direction, as depicted in Fig. 2. No-slip boundary conditions were applied along the ground surface. The velocity degrees of freedom belonging to nodes of the nozzle surface were prescribed to the value of the printing speed to yield a noslip boundary condition relative to the printing speed v p (see Fig. 2). The velocities of the fluid nodes at the top surface of the extrusion nozzle (fluid inlet) were set to the vector sum of the printing speed and the inlet flow velocity, which was obtained from the analytical solution of a Bingham flow in a pipe [55] and from the inlet flow rate as Q = q · d · h · v p , where q is a scaling factor and h is the printing height (see Fig. 2). The initial conditions of the fiber orientation state at the inlet were defined in the shear layer near the boundary (dark red region in Fig. 2b) as the steady-state solution and in the plug zone in the channel center (light red region in Fig. 2b) as the random isotropic fiber orientation state. The size of the shear layer and the plug zone were obtained from the analytical solution of a Bingham flow in a pipe [55]. Note that within the parameter range investigated, the influence of the anisotropic viscosity model on the flow field in a pipe is negligible compared to the flow field from the analytical solution of a standard Bingham model. Due to the prescribed steadystate fiber orientation inlet condition, the obtained results are only marginally influenced by the height of the extrusion nozzle H , which was scaled by the nozzle diameter as H = 2·d to yield the same relative height for different nozzle diameters. As discussed in [56], any other fiber orientation inlet condition other than the steady state solution, e.g., random isotropic, would make the final fiber orientation in the printed layer highly dependent on the modeling height of the extrusion nozzle H . In addition, in this case, a steady fiber orientation state in the cross-section would only be obtained after the particles have traveled through the entire extrusion nozzle, which would be numerically extremely expensive. Table 1 provides an overview of all conducted computational simulations discussed in the following sections. This includes simulations with the index S1-S5 with varying nozzle diameters for validation, simulations SN1-SN6 for discussing the influence of the number of explicit time step splits m max and S6-S12 for studying the influence of the process parameters such as the printing speed, flow rate, and printing height on the fiber orientation state in a printed layer. Unless otherwise specified, the number of explicit time step splits for integrating the evolution equation for fiber orientation state was m max = 1000, and the anisotropy ratio R η was calculated as described in Sect. 3.2.
The material and fiber flow properties were chosen asμ = 10 Pa s,τ 0 = 28.4024 Pa, φ s = 0.45, ρ = 2100 kg/m 3 , c = 0.01, r = 30, and κ = 0.2 (yielding μ = 10 Pa s and τ 0 = 300 Pa when c = 0), which were assumed based on typical values found in the literature [8,26,57], supported by a model calibration with respect to the experiments provided in [8], as discussed in Sect. 5.4. The time step size was specified as Δt = 0.004 s, and layers of 50 cm length were printed.
The simulation was initialized with an extrusion nozzle filled with printing material and a row of particles right outside the nozzle outlet (see Fig. 2). As described in Sect. 4.1, in the PFEM, the finite element mesh moves together with the material. To realize an inlet flow in a Lagrangian description, new particles must be added in a certain region during the simulation. The most practical approach was found to add these particles directly under the extrusion nozzle and perform a local re-meshing only in this region. New particles were added, when the particles within the nozzle traveled approximately one element size, which varied between the  The values in bold indicate the values that were modified in the analyses simulations according to the process parameters. In this regard, also the internal extrusion nozzle particles (red region in Fig. 2b) were moved back to their initial position in relation to the printing speed to create the necessary space for the new particles under the nozzle and to avoid problems related to particle packing in the nozzle. The solution variables of the particles were then updated according to the new positions based on the previous discretization. This procedure is independent of the PFEM-related re-meshing process. The process of particle insertion is illustrated in Fig. 3 in the vicinity of the extrusion nozzle outlet. The finite element mesh is shown in Fig. 2 for the simulation S1. It is refined at the boundary between the shear and plug zone, to accurately represent the initial fiber orientation state in the nozzle. Approximately 60,000 elements were used for the extrusion nozzle and 200,000 elements for the printed layer. The number of elements slightly varies with the nozzle geometry and printing parameters due to geometrical differences.
The discussion of the results in the following sections are organized as follows: Sect. 5.2 provides a general analysis of the fiber orientation and characteristics of the simulation based on the reference simulation S1 and serves as an introduction for the discussions of the remaining simulations. In Sect. 5.3, the influence of the number of explicit time step splits m max on the accuracy of the predicted fiber orientation states is discussed to find an appropriate value for m max based on Simulations SN1-SN6. Section 5.4 provides the analyses of varying nozzle diameter according to simulations S1-S5 and is split into two subsections. The first subsection provides an analysis of the fiber reorientation and flow mechanisms during printing for the different nozzle diameters and the second subsection is focused on a comparison with experiments for validation and a discussion about the limitations of the model. The first subsection is also intended to answer the question, whether smaller or larger nozzles are preferred to obtain a stronger fiber alignment in printing direction, which would allow a better control of the fiber orientation in printed components. Finally, Sect. 5.5 discusses the influence of the printing speed, flow rate and printing height on the fiber orientation based on simulations S6-S12 and aims to answer the question whether and to which extent the fiber orientation state in a printed layer can be controlled with the different process parameters.

Analysis of fiber reorientation and flow mechanisms of the reference simulation S1
This section covers results from the numerical analysis of the printing process regarding the fiber orientation state during extrusion according to the reference simulation S1 (see Table 1). The analysis focuses on revealing the influence of flow mechanisms during the printing process and the resulting fiber orientation states in the printed layers. The discussion of the results is intended to provide a basic understanding of these processes for the analyses of the remaining simulations in the next sections. Due to the prescription of the steady-state fiber orientation state as an initial condition at the inlet (top of the extrusion nozzle), a steady state fiber orientation state was obtained in the printed layer after only a few seconds of printing time. The simulated fiber orientation state is depicted in Fig. 4a and in Fig. 4b showing the fiber orientation along the symmetry plane for simulation S1 (parameters given in Table 1). The orientation tensors are visualized as ellipsoids and projected on a uniform grid in the post-processing step, as described in [16]. A needle-like ellipsoidal shape can be interpreted as a high probability of fibers oriented in the direction in which the needle-like ellipsoid is pointing. Figure 4a and b also show the orientation tensor component a 2,x x , which can be interpreted as a tendency of fibers orienting in x-direction, i.e., the printing direction. As can be seen, large values of a 2,x x , and needle-like ellipsoids pointing in x-direction, indicate that fibers tend to align in printing direction (x-direction) after deposition. During the extrusion, the fibers mainly oriented in the flow direction.
The steady-state fiber orientation inlet condition is characterized by a jump of the orientation tensor between the shear and plug zones, as shown in Fig. 4b). Fiber alignment in printing direction (x-direction) is most pronounced at the boundaries of the printed layer due to high shear stresses in flow direction in the shear layer near the boundaries, as well as the already high level of alignment of fibers resulting from the initial condition based on the steady-state solution of flow in the pipe. A less pronounced alignment of the fiber orientation in printing direction (x-direction) is found in the plug zone (center of the printed layer) due to a random isotropic fiber orientation, associated with lower shear stresses in the center under the extrusion nozzle. However, even the smaller shear stresses in the plug zone (center of the printed layer) are leading to a weak fiber orientation in printing direction (x-direction) already at a short distance after leaving the nozzle, as indicated by the elongation of the ellipsoids in printing direction (x-direction) in Fig. 4b at x ≈ −3 cm and z ≈ 0.5 cm.
The fiber orientation tensor ellipsoids and values of the component a 2,x x of the printed cross-section are provided in Fig. 4c) for simulation S1. Besides the characteristic jump of the orientation tensor between the shear and plug zones, an influence of shear stresses in lateral (y-) direction in the plug zone (z/h ≈ 1/2) can be observed, as is manifested by the elongation of the ellipsoids. The influence of shear stresses in lateral direction in combination with strong shear stresses in printing direction (x-direction), as observed with respect to Fig. 4b, is leading to a slight planar fiber orientation state within the xy-plane in the lower part of the plug zone Initial state at the nozzle outlet: State after few time steps: State after particle insertion: particles within the extrusion nozzle particles outside the extrusion nozzle newly inserted particles   [16]. The orientation number (Eq. (31)) was evaluated based on the averaged orientation tensor over the cross-section at the end of the simulation as N o = 0.7886, and it can be interpreted as a single value characterizing the fiber orientation state in the printing direction (x-direction). Since the orientation number can take a value between zero (no fibers are oriented in x-direction) and one (all fibers are oriented in x-direction), a value of 0.7886 means a strong alignment of the fibers in printing direction.
The analysis of Fig. 4a-c is concluded by the observation that the fiber orientation in the printed cross-section is very heterogeneous and the result of a complex interaction between the initial fiber orientation within the extrusion nozzle and the experienced flow and shear history of the material under the nozzle during extrusion. The fiber orientation within the cross-section is analyzed in more detail in comparison to results from additional computational simulations in the next sections.

Influence of the number of explicit time step splits m max
For simulation S1 (see Sect. 5.2), the number of explicit time step splits m max was 1000. This choice was accompanied by a convergence study of the parameter m max , carried out beforehand to determine the value required to obtain good accuracy. A discussion of this convergence study is given in the following. The averaged orientation number N o of the printed cross-section was used to assess convergence of the fiber orientation model with respect to the number of explicit time step splits m max . Figure 5 shows the error of the orientation number N o based on simulations SN1-SN6 assuming an isotropic constitutive model with R η = 1 and different explicit time step splits m max . The assumption of an isotropic constitutive behavior was made to obtain a mesh independent convergence study, as the influence of fiber orientation on the flow field disappears, which resulted in exactly the same mesh for all simulations at each time step. Note that for the parameter ranges studied in this work, the difference between the results of the anisotropic (here R η ≈ 1.12, calculated as described in Sect. 3.2) and isotropic (R η = 1) constitutive assumption were only marginal. The error was calculated based on N o,err or = (N o − N o,SN 6 )/N o, SN 6 with the value of the orientation number N o,SN 6 obtained from simulation SN6 using a large value of m max = 3000. As can be observed, a value of m max = 10 already leads to an averaged orientation number, which only marginally differs from the orientation number obtained for m max = 3000. Hence, it is concluded that using a value of m max = 1000 leads to a sufficient level of accuracy in the remaining simulations.

Influence of the nozzle diameter and validation
In this section, the influence of the nozzle diameter on the fiber orientation state is analyzed using the numerical results obtained from simulations S1-S5 (see Table 1), which is the focus of Sect. 5.4.1. In addition, the numerical model is validated by comparing the numerical results with experimental results from the literature in Sect. 5.4.2.

Analysis of fiber reorientation and flow mechanisms for different nozzle diameters
The following discussion is intended to answer the question, whether smaller or larger nozzles are preferred to obtain a stronger fiber alignment in printing direction and which fiber reorientation and flow mechanisms are dominating. For simulation S5, using the smallest extrusion nozzle diameter d = 1 cm, the fiber orientation state and the orientation tensor component a 2,x x are shown in Fig. 6a as a 3D detailed view and in Fig. 6b as a view on the symmetry plane of a section around the extrusion nozzle.
Comparing Fig. 6a and b with the results obtained using the larger nozzle diameter d = 3 cm in Fig. 4a and b, respectively, it is observed that a stronger fiber alignment in printing direction (x-direction) was obtained in the case with the smaller extrusion nozzle diameter. This is mainly due to substantial shear stresses in the nozzle leading to a larger size of the shear layer along the nozzle boundaries and, consequently, a stronger fiber orientation in the flow direction. The larger shear stresses in smaller extrusion nozzles are attributed to the higher extrusion velocities, which were necessary to obtain the required inlet flow rate Q = q · d · h · v p to fill the space under the nozzle. As observed from a comparison of Figs. 6b and 4b, the fibers in the plug zone slightly oriented towards the z-axis for the smaller extrusion nozzle, as is indicated by a slight tilting of the ellipsoids towards the z-axis. This is due to the larger extrusion velocities, the less constrained situation under the nozzle, and a lower ratio of nozzle diameter to printing height (d/h) in case of a smaller extrusion nozzle diameter. Based on the slight fiber orientation towards the z-axis, the effective fiber orientation in printing direction (x-direction) in the plug zone is less pronounced for the smaller extrusion nozzle.
The above observations are confirmed by inspection of Fig. 6c, showing the cross-section of the printed layer with the ellipsoidal representation of orientation tensors. Stronger fiber alignment and a more restricted (narrower) plug flow can be observed for simulation S5 with a nozzle diameter d = 1 cm as compared to Fig. 4c, which shows the simulation results using a nozzle with a larger diameter of d = 3 cm.
A more detailed presentation of the evolution of the fiber orientation state over the cross-section for the different nozzle diameters is provided in Fig. 7a and b. Figure 7a shows the fiber orientation state at the center line (z/h = 1/2) along the normalized y-coordinates (y/d) of the printed cross-section in the xy-plane (see the red lines in Figs. 4c and 6c). A stronger fiber alignment in printing direction (xdirection) is obtained for smaller extrusion nozzles at larger y-coordinates. In contrast, larger extrusion nozzles lead to a slightly higher fiber alignment for smaller y-coordinates in printing direction (x-direction), indicated by a slightly larger elongation of the ellipsoids in x-direction. This is attributed to the longer duration the material is subjected to shear and the more constrained situation under the extrusion nozzle for larger nozzle diameter to printing height (d/h) ratios. Figure 7b shows the fiber orientation state at the center line (y/d = 0) along the normalized z-coordinates (z/h) of the printed cross-section in the xz-plane (see the red lines in Figs. 4c and 6c). As already discussed with respect to Fig. 6b, smaller extrusion nozzles lead to a slight tilting of the fibers towards the z-axis and a slightly more pronounced fiber ori-entation in printing direction is obtained for larger extrusion nozzles in the plug zone (z/h ≈ 0.5). The findings from the comparative analysis above are summarized as follows: (1) Based on higher extrusion velocities and larger shear layer sizes within the extrusion nozzle, extrusion nozzles with smaller diameters lead to a stronger fiber alignment in printing direction (x-direction) along the layer edges. (2) Contrary to the first finding, due to a less constrained situation under the extrusion nozzle and higher extrusion velocities, the fiber orientation in printing direction (xdirection) within the layer center was less pronounced when using smaller extrusion nozzles.
As visually assessed above, the first effect dominates the effective fiber orientation within the printed cross-section for the examples investigated in this paper. Therefore, according to the presented results, smaller extrusion nozzles are preferred to obtain a stronger alignment of fibers in printing direction.

Comparison with experimental results for validation using different nozzle diameters
For validation and further investigation, the printed crosssection was split into two sections A (y > d/4) and B (y < d/4), as described in [8] and shown in Figs. 4c and 6c. Figure 8a depicts the averaged orientation number N o for sections A and B obtained from simulations S1-S5 in comparison to the experimental results reproduced from [8].
Larger orientation numbers can be found in Section A as compared to section B, implying that a more pronounced fiber alignment in printing direction (x-direction) is obtained at the layer edges, which is in agreement with the previous findings. Comparing the experimental and numerical results, a good correlation was found for larger extrusion nozzle diameters. For smaller extrusion nozzle diameters, specifically in section B, however, the experimentally evaluated orientation numbers show a stronger increase as compared to the orientation numbers obtained from the numerical simulations. The reasons for the observed differences between the experimental and numerical results for smaller diameters may be manifold. For example, in [8], material properties, flow rate, and geometry of the extrusion nozzle were not specified and were chosen in this work based on typical values from experience and literature. Further, in this work, the conical extrusion nozzle used in [8] was simplified as a straight pipe to simplify the boundary and initial conditions for the numerical analysis. Differences between numerical and experimental results may also be attributed to fiber dispersion, potential clustering of the fibers and boundary effects, which cannot be reproduced by the Folgar-Tucker fiber orientation model. Boundary effects imply that no fiber can appear perpendicular to a wall at a distance of half the fiber length due to geometrical reasons [34,58]. The fiber length in [8] was 6 mm, which means that within a distance of 3 mm from the wall all fibers are assumed to be aligned parallel to the wall. In all simulations, this region falls into the shear layer at the nozzle boundary and can therefore be considered to be mimicked by the initial fiber orientation state in the shear layer, which also reproduces a fiber orientation parallel to the wall. However, for smaller nozzles, boundary effects may stronger influence the results, as a larger ratio of the nozzle cross section is affected by these boundary effects. In addition, since concrete is a heterogeneous material, effects like shear induced particle migration or non-uniform structural buildup across the nozzle cross section can also influence the flow behavior [3]. From this analysis, it can be concluded that an accurate prediction of the fiber orientation during printing depends on an accurate prediction of the size of the shear layer and plug zone, which depends on the various phenomena discussed above that are very challenging to be included in a Folgar-Tucker model. Nevertheless, the results indicate that for larger extrusion nozzles, the fiber orientation state in a printed layer can be predicted precisely enough using the proposed model. ) for a different printing speeds v p (simulations S1 and S6-S8), c different flow rate factors q (simulations S1, S9 and S10) and e for different printing heights h (simulations S1, S11 and S12). Right hand side: Orientation distribution functions in the xy-plane (Eq. (32)) for different flow rate factors averaged over the printed cross-section for b different printing speeds v p (simulations S1 and S6-S8), d different flow rate factors q (simulations S1, S9 and S10) and f for different printing heights h (simulations S1, S11 and S12) Figure 8b provides the orientation distribution functions in the xy-plane (Eq. (32)), averaged over the whole crosssection, for different nozzle diameters. As can be observed, the probability of fibers having a smaller azimuth angle φ is higher for smaller nozzle diameters, which means a stronger fiber alignment in printing direction (x-direction).

Influence of the printing speed, flow rate and printing height
This section discusses the influence of the printing speed, flow rate and printing height on the fiber orientation state using the numerical results from simulations S1 and S6-S12 (see Table 1) and aims to answer the question whether the fiber orientation state in a printed layer can be controlled with the different process parameters. The numerical results are presented similarly as in Fig. 8a and b based on the orientation numbers in sections A (y > d/4) and B (y < d/4) and the orientation distribution functions in the xy-plane averaged over the whole cross-section. As indicated in Fig. 9a by the increasing orientation number with increasing printing speed v p (simulations S1 and S6-S8) in both sections A and B, larger printing speeds led to a more pronounced fiber orientation in the printing direction (x-direction). This is attributed to the higher shear stresses arising in the material during extrusion due to the higher extrusion rates, which were required to fill the space under the extrusion nozzle when using larger printing speeds. This observation is in agreement with the experimental studies in [8] and the higher probability of lower azimuth angles φ for larger printing speeds according to the orientation distribution functions in the xy-plane, as shown in Fig. 9b. Figure 9c) depicts the orientation numbers in sections A and B over the flow rate factor q (simulations S1, S9, and S10). Here, analyzing varying flow rate factors q, the situation is less clear as compared to the influence of the printing speed v p above. The evolution in section A shows a slightly decreasing trend with increasing flow rate, while the trend of the orientation number in section B is not fully clear. It should be noted that slight variations of the orientation number may also be attributed to the re-meshing procedure, which resulted in different finite element discretizations in all simulations. The decreasing fiber orientation with increasing flow rate factor q observed for section A seems to contradict the higher shear stresses in the material during extrusion when using larger flow rates. Given the higher shear stresses, more fibers are oriented in the flow direction within the extrusion nozzle. However, they were not effectively aligned in printing direction (x-direction) during extrusion. This paradox situation can be explained by the fact that the space under the nozzle is already filled when using the lowest flow rate factor q = 1.1. A further increase in the flow rate strengthens the flow in the lateral direction (y-direction), and the fibers are no longer effectively oriented in the printing direction (x-direction). This is confirmed by Fig. 9d),which shows the orientation distribution functions in the xy-plane for different flow rate factors q. A slightly lower probability is obtained for lower azimuth angles φ for increasing flow rates, which means slightly less fiber orientation in the printing direction (x-direction). Figure 9e) shows the averaged orientation number over sections A and B for different printing heights h (simulations S1, S11 and S12). Given an increasing printing height h, the available space under the nozzle and the flow rate increases based on the flow rate relation Q = q ·d · h ·v p . According to Fig. 9e), in contrast to the previous study on the influence of the flow rate factors q, the fibers are more effectively aligned in printing direction (x-direction) for increasing printing heights h. This is due to the more available space under the nozzle and less influence of the flow in the lateral direction (y-direction). Figure 9f) depicts the orientation distribution functions in the xy-plane for different printing heights h. A higher probability for lower azimuth angles φ is observed for increasing printing heights h, which means a stronger fiber alignment in printing direction (x-direction).
The investigations in this section reveal that a stronger fiber alignment in the flow direction within the extrusion nozzle, due to stronger shearing of the material, does not necessarily lead to a more pronounced fiber alignment in printing direction (x-direction) after extrusion. Only when the flow in lateral direction (y-direction) is less pronounced, this extra fiber alignment is effectively oriented in the printing (x-direction) during extrusion. Higher shear stresses within the extrusion nozzle can be obtained by using smaller nozzle diameters and/or larger flow rates, e.g., as obtained from larger printing speeds or larger printing heights. Stronger shearing under the extrusion nozzle, occurring when larger extrusion nozzle diameters are used, can also lead to a stronger fiber alignment in the printing direction as discussed in Sect. 5.4. The results discussed indicate that the fiber orientation in the printing direction can be controlled to some extent with the process parameters printing speed and printing height. Prerequisite for this is the control of the mechanisms of fiber reorientation and flow during printing discussed above.

Conclusions
A new numerical modeling approach for the flow simulations of fresh fiber-reinforced concrete with a stable solution algorithm within a Particle Finite Element (PFEM) framework has been proposed. A Folgar-Tucker fiber orientation model was adopted to calculate the evolution of the fiber orientation state. This model represents the fiber orientation state with a second-order orientation tensor. The fiber orientation model was combined with an anisotropic constitutive model, taking the current fiber orientation state of the material into account. Upscaling formulas for the rheological properties were taken from the literature. Finally, within the PFEM framework, an efficient solution algorithm for explicit coupling of the fiber orientation model and the fluid dynamics equations was proposed. The maximum entropy method was used as a post-processing tool to reconstruct the orientation distribution function from the second-order orientation tensor for a more convenient interpretation of the numerical results. The additional computational cost associated with the fiber orientation model compared to a PFEM model without fibers can be considered small, since the calculation of the second order orientation tensor is performed locally for each node and does not require an expensive solution procedure.
The numerical model was validated by comparing the simulated orientation number of a printed cross section obtained from 3D-printing processes of fiber-reinforced concrete with experimental results from the literature for different extrusion nozzle diameters. The numerical results agreed well with experimental results for larger extrusion nozzles. In comparison, slightly larger differences were observed for smaller extrusion nozzles, which can be attributed to several reasons, such as boundary effects or fiber dispersion that are not captured by the proposed Folgar-Tucker fiber orientation model, as well as by assumptions regarding the material properties, initial conditions and geometry. These effects and assumptions can influence the size of the sheared and unsheared regions during flow, which drives the fiber (re-)orientation during printing and therefore governs the quality of the fiber orientation prediction. Fiber dispersion and boundary effects may only be realistically captured by using discrete fiber modeling approaches, however, at a substantially higher computational cost.
Further parametric studies of the nozzle diameter, printing speed, flow rate and printing height were performed to analyze the fiber reorientation and flow mechanisms during printing and aim to answer the question whether and to which extent the fiber orientation state in a printed layer can be controlled with the different process parameters to allow a more precise control of the printing process. The analyses revealed that the final fiber orientation within the printed layer is a result of a "competition" between the shearing of the material in the extrusion nozzle, the shearing of the material under the nozzle, and the material flow characteristics. Stronger shearing of material under the nozzle can lead to a more pronounced fiber orientation in the printing direction, as was obtained for larger extrusion nozzle diameters. In contrast, increasing the nozzle diameter also reduced the shear stresses within the nozzle, which is contradictory to the effect of the stronger shearing of the material under the extrusion nozzle. The numerical results indicated that the shearing effect dominates the final fiber orientation state in a printed layer in the extrusion nozzle, at least for the parameter range studied in the analyses, which means that a stronger fiber alignment in printing direction was obtained for smaller extrusion nozzles. Stronger shearing within the extrusion nozzle, obtained for larger flow rates, leads to a stronger fiber alignment in the flow direction within the nozzle. However, this fiber alignment was only effectively oriented in the printing direction after deposition when the flow in the lateral direction was less pronounced, which was not the case as the flow rate increased, since the additional material could only flow in the lateral direction. On the other hand, larger shear stresses within the nozzle were also obtained using larger printing speeds and printing heights leading to a stronger fiber alignment in the flow direction within the nozzle. In contrast to the increasing flow rate, this fiber alignment was effectively oriented in the printing direction as enough space was available for the material below the extrusion nozzle and the flow in the lateral direction was less pronounced. It is concluded from the computational investigations that the fiber orientation in the printing direction can be controlled (to a certain extent) by the choice of the process parameters during 3D-concreteprinting. In particular, stronger fiber alignment in printing direction was obtained for smaller extrusion nozzles, larger printing speeds, and larger printing height.
Some aspects, such as conical extrusion nozzle shapes or multi-layer printing scenarios, have not been considered and should be investigated in future contributions. The constitutive model may be also improved by introducing upscaling relations for the rheological properties designed specifically for printable concretes. Also, a submodel which considers the influence of the fiber orientation on the yield stress is a topic of a future enhancement of the proposed model.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.