Viscoelastoplastic classification of cementitious suspensions: transient and non-linear flow analysis in rotational and oscillatory shear flows

Rheological modeling of special concretes such as ultra-high performance concretes (UHPCs), self-compacting concretes (SCCs), or environmentally friendly concretes with low clinker contents but containing many fine particles increases model complexity because high amounts of colloidal and fine particles and complex chemistry result in strongly non-Newtonian rheological behavior. Straight-forward viscoplastic rheological models such as the well-known Bingham- or Herschel-Bulkley approach do not fit the flow behavior on a large scale at and lead to boundary value problems in numerical flow simulations. Increased viscoelastic properties due to chemical admixtures are not described by pure viscoplasticity analysis methods. To fill this gap, our research presents viscoplastic and viscoelastic linear and non-linear rheological characterization methods for common and strongly non-Newtonian cement pastes and combines the methods for a comprehensive viscoelastoplastic modeling of cementitious building materials. We illustrate and discuss the benefits and boundaries of each measurement technique in dependence of cementitious paste complexity. Three solid volume fractions from ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} = 0.45 to 0.55 and three different flowability ranges, adjusted through varying superplasticizer amounts, were investigated in respect of steady-state and transient yield stress, viscosity, and non-linear viscoelastic material properties. Our results reveal that with increasing solid volume fraction and polymer amount, viscoelastoplastic modeling describes rheological flow on a large flow scale more appropriately than common viscoplastic methods. The results show a new approach for a full quantitative and qualitative rheological classification of cementitious building materials and thus improve the understanding of densely packed colloidal cementitious suspensions. The findings serve as prospective guideline to model complex cementitious building materials both at low and high shear rates, and to find an appropriate rheological description at the boundaries of flow.


Introduction
Concrete rheology describes concrete flow dependent on external shear forces and internal material structuring effects arising from physical and chemical interactions within the cementitious suspension. Cementitious suspensions are two-phase suspensions with colloidal and non-colloidal cement particles and a carrier liquid. Above the percolation threshold, cement particles flocculate due to colloidal and interparticle forces and build up a network (Barnes 1999;Firth und Hunter 1976). The particle network strength and thus its resistance against deformation and breakage depends on attractive and repulsive interparticle forces (Kee and Man Fong 1994) (Lowke 2013)) and hydration nucleation products that build up over time. The paste mixture and the ions within determine interparticle forces and hydration speed and thus the particle network strength (Flatt und Bowen 2007). Particle distances decrease with increasing solid volume fraction, causing attractive forces and frictional forces during shear (stated already by Krieger und Dougherty (1959) for cement pastes by Coussot et al. (2002)). The suspension starts to flow as soon as viscous forces exceed interparticle forces, which is defined as yield stress 0 (Barnes 1999;Genovese 2012). During shear, particle networks are increasingly destroyed. Below the critical shear rate and at rest, structural build-up of the cementitious particle network takes place dependent on attractive forces and hydration nucleation. Thixotropy is the fully reversible and thus solely flocculation-dependent structural build-up, which determines several workability properties in cementitious building materials such as additive manufacturing. Precise characterization of the yield stress 0 , viscosity , and structural build-up is essential for workability predictions, processing adjustment, and efficient numerical modeling of concrete flow.
However, varying concrete and cementitious paste properties and abundant investigation and calculation methods make the rheological material classification difficult. Inappropriate rheological set-up and analysis methods lead to inaccuracies in the rheological description (Wallevik et al. 2015), and the choice of appropriate parameter calculation methods depends on the overall material behavior (e.g. shear-thinning or shearthickening) and the focus of processing. Whilst steady-state material characterization is needed to describe fast pumping, and extrusion in additive manufacturing requires the parameter for full structural build-up during rest, transient processing description requires shear gradient and flow property analysis over a wider shear range (Wallevik und Wallevik 2017).
General rheological analysis of non-Newtonian materials ranges from Bingham's ideal viscoplastic model (Bingham 1916) to Maxwell's ideal elastic, ideal viscous, and viscoelastic classification and various classifications in between. Figure 1 illustrates two simplified constitutive models: Fig. 1a depicts an ideal viscoplastic Bingham model with a dashpot and a block. Figure 1b illustrates a rheological model as a combination of viscoelastic deformation prior to flow and viscous or viscoplastic flow once the yield point is exceeded (adapted from Souza Mendes and Thompson 2012). The complexity of non-Newtonian material analysis was investigated by Coussot et al. (2002) and reviewed extensively by Malkin (2013). In concrete rheology, straightforward viscoplastic models adequately describe concrete flow. However, special concretes such as UHPC, self-compacting concrete, or environmentally friendly concretes with low clinker content exhibit more complex rheological behavior especially at very high or low shear rates. Containing high colloid and fine particle content and chemical additives, they become strong non-Newtonian materials with discontinuous rheological properties, which is why the rheological parameters yield stress 0 and viscosity become functions of time and shear. In this case, state-of-the-art rheological modeling does not depict the material properties. On the contrary, in polymer sciences, more sensitive viscoelastic material characterization using oscillatory rheometric methods describes the full range of transient and non-linear material behavior. Since Pipkin's extensive research of viscoelastic materials (John et al. 1986), viscoelastic matter is described qualitatively through Lissajous-Bowditch curves that display stress as a function of deformation and frequency and by calculating elastic and viscous parts and their change as a function of strain and time t . Two mean characterization approaches, i.e., Fourier transformation (FT), introduced in Wilhelm et al. (1998), and stress decomposition using Chebyshev polynomials, proposed by Ewoldt et al. (2008), allow sensitive qualitative and quantitative material characterization for strongly non-Newtonian, nonlinear viscoelastic behavior. Recently, Ewoldt and McKinley extended the traditional two-dimensional Pipkin-diagram to a three-dimensional structure-implementing analysis in order to characterize thixo-elastoviscoplastic (TEVP) material behavior in Ewoldt und McKinley (2017), and also compared simple shear with viscoelasticity analysis (Ewoldt et al. 2010).
Certain researchers investigated viscoelastic properties in cement and concrete as well (Haist 2009;Nehdi and Al Martini 2007;Roussel et al. 2012;Struble et al. 2000), and approaches for comparing calculated yield stresses and oscillatory viscoelastic analysis can be found in Ukrainczyk et al. (2020) and Yuan et al. (2017), for example. Conte and Caouche analyzed cement pastes beyond the viscoelastic regime in Conte und Chaouche (2016) using Lissajous-Bowditch figures and quantitative techniques. However, transient and non-linear rheological analysis for different cementitious building materials is not state-of-the art in concrete sciences, as was recently stated in Roussel et al. (2019). Also, a comparative approach between traditional dynamic steady-state and complex viscoelastic material analysis is Fig. 1 Idealized rheological models: a Type I ideal viscoplastic Bingham model with a Newtonian damper and a St. Venant-element, b Idealized viscoelastoplastic model with a spring stiffness for elasticity prior to viscous or viscoplastic flow still lacking, even if sensitive rheological techniques allow understanding and adjusting complex cementitious materials. A clear understanding however will support prospective complex cementitious building material development and bring forward computational fluid dynamics in this field, as already successfully implemented for soft gels, wormlike micelles (see Kim et al. 2013) or blood aneurisms (see Elhanafy et al. 2019).
In this publication, we aim to phenomenologically model complex cementitious pastes as viscoplastic and viscoelastic suspensions to benefit from a combined viscoelastoplastic modeling approach. We discover paste-dependent appropriate calculation methods and suitable material parameters, and also benefits and boundaries of different analysis setups. A generalized property map with material-relevant parameters shall give an overview of the whole range of cementitious paste flow properties and prospectively shall serve as guideline for precise rheological analysis of complex cementitious mixtures.

Steady-state modeling
In "true yield stress" analogies, cementitious building materials are considered as viscoplastic suspensions with a finite yield stress and Newtonian or non-Newtonian behavior beyond yield. ̇− -flow curves easily describe steady-state phenomenological flow behavior. The most common rheological scalar constitutive equations for cementitious building materials are the Bingham model, which describes Newtonian fluids with a yield stress ( = B +̇ ; with B = Bingham yield stress in (Pa), = plastic viscosity in (Pas), and ̇ as shear rate in (s −1 )) and, more sufficiently, the Herschel-Bulkley model, which, once a yield stress is exceeded, describes shear rate-dependent shear-thinning or shear-thickening behavior: with 0,HB = Herschel-Bulkley yield stress in (Pa), k as flow coefficient in (Pas), n as Herschel-Bulkley-index with n < 1 for shear thinning behavior, n = 1 as ideal viscous behavior, and n > 1 for shear thickening behavior. Viscosity curves describe yield stress-independent viscosity functions: the cross model (Eq. (2)) is a common viscosity model, originally invented for the characterization of uncrosslinked polymers: (1) = 0,H−B + k̇n with (̇) as shear-dependent viscosity, 0 as zero shear-rate viscosity in (Pas), ∞ as infinite shear-rate viscosity in (Pas), and c and p as material fitting constants. Advantageously, this model can describe the viscosity in very slow, transient, and fast flow conditions through the implementation of a zero-and infinite shear-rate viscosity as boundaries, but has the drawback of not considering a yield stress in its mathematical function nor elastic recovery after stress relief; also, it does not contain shear-thickening behavior or viscosity bifurcation description (Haist 2009).

Transient and discontinuous modeling
Below the critical shear rate, agglomeration and flocculation processes lead to an increase in microstructural strength and thus a shear stress increase at constant shear rates.
Material models for the implementation of structural build-up, viscosity bifurcation, and thus discontinuous flow exist in different disciplines, e.g. Bénito et al. (2008), Le-Cao et al. (2020), Mujumdar et al. (2002), Souza Mendes (2009), Souza Mendes and Thompson (2013, and Zarei und Aalaie (2020). Often, phenomenological expressions are combined with microstructural formulations to explain the increase in shear stress through agglomeration processes (Kee and Man Fong (1994), Le-Cao et al. (2020), Phan-Thien und Mai-Duy (2017)). In cement and concrete sciences, Roussel proposes a phenomenological equation for the stress based on structural processes: with a structural parameter ( ) , dependent yield stress 0 , and a common Herschel-Bulkley equation beyond the startup of flow. The structural build-up parameter is defined as where Θ is the material-dependent flocculation time, and and m are (de-)flocculation-dependent fitting parameters (Roussel 2006).
Equation (3) and Eq. (4) depict an indirect phenomenological fitting approach; the structural build-up parameter has no clear physical foundation. This was critically reviewed by de Souza Mendez et al. in Souza Mendez und Thompson (2012), who then proposed a viscosity-based steady-state model with a critical bifurcation shear rate ̇c rit and defined dynamic and static yield stress: ,0,S − ,0,ḋ e −̇∕,c + ,0,Ṡ + k̇n −1 + ∞ where 0,S is the static yield stress, 0,d the dynamic yield stress, ̇c the critical shear rate, k the consistency index, and n the power-law index. The proposed Eq. (5) avoids a fitting-based structural parameter , but has seven variables of which at least four parameters need to be fitted ( k , n , 0,d , 0,S ) and three parameters should be calculated in advance ( 0 , ̇c, ∞ ) to reduce the fitting error. Furthermore, Eq. (5) only describes shear-thinning material behavior, wherefore not all cementitious materials can be adequately described. Therefore, Zarei et al. propose a case-dependent discontinuity function in (Zarei und Aalaie 2020). Single functions f I and f II distinguish the viscosity evolution below and above the critical shear threshold ̇c rit : The functions f I and f II contain rheological models that best fit the experimental rheological behavior, e.g., the cross-model (Eq. (2)) or the de Souza Mendez function (Eq. (5)). Zarei et al. propose different functions to correctly predict shear-thickening rheological behavior with similarities to Eq. (2) in Zarei und Aalaie (2020): In Eqs. (7) and (8), five resp. six parameters need to be known, from which only two parameters (c and n I in Eq. (7) and c and n II in Eq. (8)) are fitted if the boundaries of the viscosity function and the maximum investigated shear rate ̇m ax are known. Equations (7) and (8) enable the rheological f low modeling of cementitious pastes over different ranges in a discontinuous way. Further kinetic, phenomenological, and structural-based flow descriptions have been the subject of research both in complex fluid and concrete science. A literature review is beyond the scope of this paper, but can be found elsewhere (e.g., (Jiao et al. 2021;Le-Cao et al. 2020;Mewis und Wagner 2009a, 2009bRuckenstein und Mewis 1973;Souza Mendes 2011;Tattersall und Banfill 1983)).

Viscoelastic rheological characterization
Compared to dynamic flow rheology, oscillatory measurements allow a sensitive material characterization of viscoelastic material properties especially at very low deformations. At a constant angular frequency , the amplitude of the load, i.e., strain amplitude 0 or stress amplitude 0 , is increased with time. The system response at each measuring point (for strain amplitude tests, the stress response , and vice versa) is determined as a function of angular frequency and time. In strain-controlled amplitude tests, the oscillatory strain is ( , t) = 0 sin( t) , with 0 as strain amplitude, as frequency in (rad/s), and t as time in (s). The strain rate ̇ is 90° out-of-phase of the inserted strain : ̇(t) = 0 cos( t) . The shear response in fully elastic materials is in phase with the applied strain, which corresponds to the Hooke elasticity law (t) = G * (t) , where G * as complex shear modulus is the proportionality constant in (Pa) between (t) and . The stress response in fully viscous materials is in phase with the strain rate ̇ . The rheological model corresponds to the Newtonian law (t) = * ⋅ (t) with the complex viscosity * . The shear modulus G * decomposes into a viscous and elastic shear modulus, where the storage modulus G ′ corresponds to the elastic retention and the loss modulus G ′′ as complex component corresponds to the viscous dissipation, both in (Pa). Nearly all materials are viscoelastic, which leads to a partial stress response in phase with the strain and a partial stress response in phase with the strain rate ̇ , or calculated as (t) = 0 (G � sin t + G �� cos t) . The phase shift angle (°; rad) indicates the phase shift of the system response compared to the initial oscillatory load and is calculated as tan(δ) = G �� G � . is a simple indicator for the viscoelasticity of a material: The larger the phase shift angle , the more viscous the material. Figure 2a shows the normalized schematic system input as function of and t and resulting stress for three different input strains. In Fig. 2b, the timeand frequency-resolved γ − τ-relationship is plotted, also referred to as Lissajous-Bowditch curves (LB curves).
The enclosed area of one LB cycle stands for the dissipated energy within (Ewoldt 2009). In the ideal-elastic case, the ellipse becomes a line. In the ideal viscous case, the ellipse is a circle. For viscoelastic material, the phase shift angle can be interpreted as the tilt angle of the ideal circle to an ellipse. As soon as the induced strain results in particle destruction, non-sinusoidal stress output follows a sinusoidal input deformation (see Fig. 2a and b, transient and non-linear curves). LB curves then appear backwardand forward-tilted (see Fig. 2b, transient and non-linear plot). Mathematical methods, e.g. Fourier transformation, describe the non-sinusoidal outputs: In order to analyze and interpret the system response and thus the material behavior, the non-harmonic periodic signal is decomposed into several sinusoidal, i.e., harmonic, oscillations (Mezger 2016). In non-linear rheology, the Fourier transformation provides harmonic decompositions from which rheological quantities can be obtained, e.g., higher-order moduli, but to which no clear physical meaning can be assigned when considered alone. (Hyun et al. 2011).

Linear viscoelastic material analysis (inter-cyclic analysis)
Linear viscoelastic analysis provides meaningful material information in the non-destroyed regime. The amplitude sweep captures the change in storage and loss modulus of the first harmonic decomposition with increasing 0 (also referred to as inter-cyclic material analysis). Figure 3 shows the schematic first-order storage and loss moduli during a small amplitude oscillatory shear (SAOS) protocol. The − G diagram with deformation on the x-axis and the resulting first-order storage and loss moduli on the y-axis classify the linearviscoelastic regime (LVE), the transient yielding zone, and non-linear material behavior after the cement paste starts to flow.
Within the LVE regime, the storage modulus and loss modulus possess almost constant values. If the storage modulus G ′ in the LVE region is greater than the loss modulus G ′′ , the substance is in the solid or gel state because the elastic component predominates. Conversely, if G ′′ > G ′ , the material is in the liquid state.
In the LVE range, the storage modulus G ′ , the loss modulus G ′′ , and the corresponding provide sufficient information about viscoelasticity. The yield point at the linearity limit L is determined as the value above which the storage modulus drops significantly (Mezger 2016). The yield stress F at the flow deformation F , also referred to as gel point, is the shear stress at the intersection of the G ′ and G ′′ -curves. A direct calculation of the yield stress F , corresponding to the yield stress 0 , is not possible in a strain-controlled oscillatory amplitude test because F generally is beyond the linear range and contains a non-harmonic, non-sinusoidal system response. The region between and F is the yielding zone. The closer the yield index is to 1, which is the ratio of F and , the more brittle material fails. Beyond F , the material is assumed to flow.

Intra-cyclic transient and non-linear viscoelastic analysis
Beyond the LVE regime, higher harmonics of the Fourier Transformation need to be considered. Fourier-transform analysis or stress decomposition with Chebyshev polynomials enables a mathematical analysis for every LB curve of each given strain amplitude 0 (also referred to as intracyclic material analysis). Fourier-transformation rheology (FT-rheology) in general, but especially for polymers, was intensively investigated by Wilhelm et al. in Hyun et al. (2011) and Wilhelm et al. (1998), applied by various researchers and comprehensively reviewed by Hyun et al. in [21]. In FT-rheology, the non-sinusoidal oscillations (see Fig. 5) are decomposed into a row of harmonic, sinusoidal viscous, and elastic parts: with 0 as strain amplitude and G ′ n and G ′′ n as odd storage and loss components of the odd higher harmonics n. The higher order harmonic components represent the nonlinearity. An intensity ratio I 3∕1 between the third and first harmonic displays the non-linearity material properties: during the linear viscoelastic response, I 3∕1 ≈ 0.
Besides FT-rheology, harmonic analysis using stress decomposition and Chebyshev polynomials was investigated by Ewoldt et al. (Ewoldt 2009;Ewoldt et al. 2008). The general assumption is a stress decomposition into an elastic and a viscous stress: Orthogonality allows the individual polynomials, and thus the elastic and viscous properties of one order, to be considered independently of those of the other orders. The coefficients e n and v n describe the elastic and viscous behavior of the nth order (Ewoldt et al. 2007) whereas often, the third-order coefficients e 3 and v 3 classify the material sufficiently into Ewoldt et al. describe further non-linear material parameters based on LB curves: The elastic material behavior is displayed by a tangent module G 0,M (with M as minimum strain) at = 0 and a secant module G 0,L (with L as large strain) at = 0 . These modules are directly related to the Chebyshev coefficients e 3 and v 3 (Fig. 4): (10)

Material modeling and material matrix
Viscoplastic material parameters presented in "Viscoplastic rheological characterization" together with viscoelastic and (t) illustrated as smoothed normalized values material parameters presented in "Viscoelastic rheological characterization" allow a precise material characterization. Depending on the point of interest, flow curves or intercyclic material analysis describe suspensions on a macroscopic scale. Shear rate-dependent structure analysis or intra-cyclic viscoelasticity analysis methods provide more sophisticated material information on a microscopic scale. The relevant rheological parameters are collected in a rheological property map in Fig. 5a: Depending on cementitious paste parameters such as solid volume fraction and polymeric content (through the addition of superplasticizers), elastic and plastic indices prior to flow and viscoelastic deformation and viscoplastic flow parameters during load define the overall viscoelastoplastic material behavior. Figure 5b presents a macroscopic spatial approach for the definition of viscoelastoplasticity for cementitious pastes through the definition of phase shift angle , non-Newtonian index n , and the yield index.
The experimental program presents the material characterization of different cementitious pastes for all methods from Eq. (1) to Eq. (16) and finally maps viscoelastoplasticity acc. Figure 5. To analyze the effect of paste composition, nine cementitious paste series were investigated using dynamic and oscillatory shear rheology at increasing load and analyzed with different mathematical models in order to classify their rheological behavior, schematically shown in Fig. 6.
As cementitious suspensions satisfy the full range of behavior from shear-thinning to shear-thickening dependent on their solid volume fraction and added chemicals, investigations were conducted for three different solid volume fractions (ϕ = 0.45, ϕ = 0.48, and ϕ = 0.55). Flowability was investigated and compared through mini slump flow measurement according to Norm DIN EN 1015-3 and analytical yield stress conversion 0A,R acc. Roussel et al.: where 0 is the yield stress in (Pa), the density of the tested material, V the tested volume, and R the final slump flow radius. For the sake of flowability comparability, three flowability values were defined by increasing PCE dosages (corresponding to slump flow values of 200 mm ± 5, 250 mm ± 5, and 300 mm ± 5) for each solid volume fraction, respectively. The mixture compositions are given in Table 1. For the experimental investigations, Ordinary Portland Cement (OPC) CEM II 42.5 A-LL N (Heidelberg Cement, Germany), demineralized water with temperature adjusted for a constant paste temperature of 20 °C, and a PCE superplasticizer (BASF, Germany) were prepared. For each test series, 0.5 l of paste was prepared by mixing water and dry cement using a hand mixer with a 4-bladed screw for 90 s at 1700 rpm. PCE was dosed 30 s after mixing had started. The paste was left at rest until 12:30 after water addition. An external pre-shear of 30 s with the hand mixer led to a pronounced de-flocculation and erased the effect of structuration history. Rigorous research into the effect of pre-shear time on rheological parameters was published by the author et al. in Thiedeitz et al. (2020). Rheometric measurements then were started directly after external pre-shear 15 min after water addition. All cementitious pastes were investigated using linear and non-linear rotational and oscillatory rheometry (schematically depicted in Fig. 6). The rheological measurements were conducted with an Anton Paar Rheometer MCR 502 (Anton Paar, Germany). Parallel plates with a diameter of 50 mm and a serrated surface to prevent wall-slip were used, and the gap was set to h = 1 mm. Raw data handling for the different measurement techniques is given subsequently.

Rotational dynamic shear rheometry
Dynamic steady shear analysis was conducted using rotational decreasing step-rates from ̇ = 80 s −1 to ̇ = 0.02 s −1 , with each step having a duration of 6 s (see Fig. 7). The shear profile was chosen in order to grasp the rheological response of cementitious pastes over a large applied range of shear rates. Prior to the step-rate sequence, a pre-shear rate of γ˙ = 40 s −1 with a duration of 30 s took place to erase heterogeneous placement effects into the rheometer (see Fig. 7a and c). Figure 7a presents the decreasing step-rate protocol with time. A rheometer measures the resulting torque after an induced rotational speed . The shear rate ̇ at the radius r is then calculated as a function of ω as ̇= r H , and the shear stress τ is calculated from torque T as (R) = 2T R 3 (Souza Mendes et al. 2014). Thus, towards the outer radius R , shear rates ̇(r) increase and also vary over the height H (see Fig. 7b). The larger the gap (especially from H > 1 mm), the greater the risk of inhomogeneous shear behavior. For Newtonian fluids, calculation errors from raw data to rheological data have been investigated extensively for different shear field assumptions and conversion factors. Assuming steady shear conditions, Mezger proposes an effective average shear rate calculation of ̇ at r = 2 3 R , which leads to conversion formulas for averaged shear rate ̇m and shear stress m of ̇m = 2 3̇( R) = 4nR 3H and m = 2 3 (R) = 4T 3 R 3 . This conversion approximation was analyzed for several benchmark fluids and cementitious paste in a Round Robin test including the authors, published in Haist et al. (2020). It was shown that especially for cementitious pastes, the  chosen test geometry and conversion function affect the rheological data, but are comparable with Newtonian fluids as long as steady shear conditions are assumed. The gap height of 1 mm and serrated plate surfaces prevent both plug flow and wall slip, wherefore the aforementioned conversion equations were chosen, as also indicated in the shear-rate conversion in Fig. 7a and the shear stress conversion of raw data in Fig. 7c. All flow curves reached equilibrium measuring conditions within one rate step within steady-state flow zone (see Fig. 7d, left). Steady-state flow represented the shear range without material restructurations. The critical shear rate ̇c rit was calculated as the shear rate beyond which the applied shear led to increasing shear stress (or, without conversion, torque) within one constant rate step (see Fig. 7d, right). Flocculation below critical shear can affect the shear field homogeneity due to possible shear banding effects (intensively discussed in Ovarlez et al. (2009) andFall et al. (2009) and schematically shown in Fig. 8).
For the presented research, a shear rate correction due to varying shear fields was not conducted, as a distinguished shear field localization is not possible, and the gap height of 1 mm was assumed to be sufficient to reduce shear field errors to a minimum. The average shear stress for each rate step in equilibrium condition was calculated from the equilibrium of each rate step; structural build up within one rate step below critical shear was included into transient yield stress analysis as time-dependent parameter. After raw data conversion, flow curves were analyzed regarding their steady-state and transient viscoplastic flow properties (Eqs. (1) to (9)). Model regression was conducted using python protocols with optimization algorithms.

Oscillatory viscoelastic rheometry
Viscoelastic analysis in the linear, transient, and non-linear regime was conducted using Large Amplitude Oscillatory Shear (LAOS) rheometry with parallel plates, as previously presented by the authors in (Prakesch and Thiedeitz 2021). A strain-controlled protocol with a strain amplitude increase from 0 = 0 to 100% and constant frequency at = 10 rad/s was used, as illustrated in Fig. 9a. The pre-shear time prior to the oscillatory measurement equaled the pre-shear from dynamic measurements (not depicted in Fig. 9a), which ensured a fixed deflocculation state of each paste. Fig. 7 Raw data handling in rotational parallel plates rheometry: a Dynamic shear profile with decreasing step-rates, b geometry-dependent shear field assumption, c flow curves for all test series, and d steady-state and slow flow regime of c Generally, the rheometer software conducts oscillatory cycles with one oscillatory input pair of , 0 until a steady material state is reached; subsequently, a set number of cycles is oversampled, as explained in detail (Ewoldt et al. 2010;Hyun et al. 2011). With cement paste being an agglomerating and chemically hardening material, steadystate regulation of the rheometer software was limited to ten oscillatory cycles. The resulting data were oversampled for each strain amplitude 0 , the results shown exemplarily for the test series 0.55-200 in Fig. 9b. To calculate linear viscoelastic parameters, the freely available software MITlaos, provided by Ewoldt et al. (2007), was used. It contains an optimization algorithm for the Fourier transformation of input raw data pairs ( 0 ; 0 ), which results in smoothed time-resolved strain-stress-cycles as exemplarily shown in Fig. 9c and d. All raw data cycles were analyzed regarding linear and non-linear viscoelastic parameters as explained in "Viscoelastic rheological characterization" using the software MITLaos, python, and MATLAB protocols.

Viscoplastic rheological modeling
All flow curves of dynamic rotational rheometry are plotted (log-log) in Fig. 10. Figure 10a−c show the flow curves for each solid volume fraction, distinguished by their flowability value. Changing non-Newtonian behavior is observed with increasing solid volume fractions and independently of different flowabilities: whereas the = 0.45 fraction generally behaves in a shear-thinning manner, = 0.48 approaches Bingham-like flow. Pastes with the highest solid volume fraction = 0.55 exhibit shear-thickening behavior. Viscosity bifurcation below critical shear increases with increasing solid volume fraction. Viscosity curves are given in Fig. 10d−f. They clearly show differing non-Newtonian behavior at high shear rates. Pronounced discontinuity due to structural build-up below critical shear is shown for each viscosity curve. Viscosity increase at very low shear rates is pronounced especially for the highest solid volume fraction.
Viscoplastic modeling of the phenomenological flow curves is shown in Table 2 and Fig. 11. Steady-state description beyond critical shear for shear stress (Herschel-Bukley, Eq. (1)), and viscosity (Cross, Eq. (2)) is presented in Fig. 11 with the calculated material parameters given in Table 2. Figure 11 shows the results for the whole shear range with the magnified excerpts for steady state. With increasing solid volume fraction , the non-Newtonian index n increases independently of flowability. Flowability, indeed, affects yield stress and viscosity ranges: With increasing superplasticizer amount, yield stress 0,H−B , 0 , and ∞ decrease. On the contrary to the Herschel-Bulkley model, the Cross model describes viscosity at the boundaries of flow (zero shear-rate viscosity 0 and infinite shear-rate viscosity ∞ ), which is a benefit concerning numerical flow simulation. However, it is only feasible for steady-state shear-thinning viscosity description, wherefore it is presented solely for the test series 0.45 (all) and 0.48-200 and additionally for 0.48-250 to illustrate application limits (Fig. 11d−e). The zero shear-rate viscosity 0 decreases with increasing flowability as does the infinite shear-rate viscosity ∞ , whereas ∞ increases with increasing solid volume fraction. Due to the differences between 0 and ∞ , c also decreases with increasing flowability. An increasing p value indicates a low change of viscosity with increasing shear rate, wherefore the applicability for Newtonian or Bingham-like flow is not given: p approaches infinity. The test series 0.48-250 thus depicts the boundary of the model. Below critical shear, material behavior is not described. For pastes with low interparticle forces and thus little or no restructuring effects, the mathematical description is sufficient. As soon as secondary effects as interparticle forces lead to higher network structures, low shear rheology is not depicted.
Structural build-up included yield stress analysis as described in Eqs. (3) and (4), and discontinuous viscosity modeling as described in Eqs. (7) and (8) is presented in Table 3 and Fig. 12.
As Table 3 shows, the initial structural build-up parameter generally increases with increasing solid volume fraction and stiffness and is greater for visibly more pronounced structural build-up below the critical shear (see Fig. 12). is the flocculation parameter, which generally increases as the shear stress increases more quickly and thus stands for shows that the zero shear-rate viscosity 0 , the viscosity at critical shear ̇c rit , and infinite shear-rate viscosity ∞ aligned with the viscosity fitting parameters c , n I , and n II , describe viscosity in both shear zones below and beyond critical shear, and, thus, allow a viscosity characterization for shearthickening material behavior and give a regression function for low shear zones.
Summarized, dynamic simple shear rheology in combination with viscoplastic material models grasps the rheology of cementitious suspensions at lower solid volume fractions and without consideration of low-shear zones. At higher solid volume fractions, rheometry and raw data conversion get more unsecure, and also mathematical modeling contains challenges. At low or no shear, rheology is weakly described and might contain errors due to secondary effects such as interparticle forces. Discontinuity functions can display viscosities at the boundaries of shear but also contain uncertainties due to probably heterogeneous shear fields and thus raw data conversion errors. Viscoelastic modeling, thus, presents another approach to describe rheology of complex suspensions in linear, transient, and non-linear state.

Linear viscoelastic analysis
Qualitative and quantitative macroscopic viscoelastic analysis for all test series is shown in Table 4 and Figs. 13 and 14. Table 4 gives all relevant viscoelastic parameters: The storage modulus G ′ , loss modulus G ′′ , phase shift angle , and complex viscosity * are given at the linearity limit L . The yield stress L at the linearity limit is calculated according to the Hookean model. The gel point F is the deformation at the first measurement point at which G ′ <G ′′ . The linearity limit L was defined as the deformation with a storage modulus deviation of 10% from the plateau, which was adapted from existing standards. The yield index was calculated to determine the elasticity of the material. With decreasing yield index, elasticity increases. The qualitative linear-viscoelastic analysis for the first harmonic storage and loss modulus and phase shift angle is given in Fig. 13.
At low deformations, the storage modulus G ′ exceeds the loss modulus G ′′ . All suspensions are gel-like with a pronounced elastic part that increases with increasing solid volume fraction . G ′ and G ′′ first show maxima and then drop sharply, which is caused by initial percolation of the particle network and cannot be prevented due to constantly evolving material changes. The phase shift evolution indicates viscoelastic material properties at increasing strain load. In Fig. 13a, the strong increase of δ corresponds to brittle failure. The evolution of δ changes significantly with increasing ϕ : A constant is followed by a increase at higher 0 for   Fig. 14 shows first-order storage moduli (see Fig. 14a−c and phase shift angular evolution (see Fig. 14d−f), for all test series distinguished by flowability value, respectively. With increasing and flowability, viscoelasticity develops from brittle failure and thus strong viscoplasticity (see 0.45-200) to large elastic response and thus more pronounced viscoelasticity (see 0.55-300).
With increasing , decreases indicating a more pronounced elasticity, independently of the paste flowability. Figure 14d−f shows that the phase shift angle evolution illustrates viscoelasticity and its change with increasing load clearly. The linear viscoelastic limit can be calculated clearly as the end of phase shift angle linearity. Subsequently, paste properties in the transient yielding zone are clearly distinguished by the phase shift angle comparison. With increasing flowability, G ′ decreases (despite higher solid volume fractions). This is explained due to the superplasticizer amount and its structure. Hot has determined (Hot et al. 2014) that the retained elastic strength of the particle network is much more extensive with high polymer content. With an increasing polymer amount, particle distances increase, which reduces direct interparticle forces (Flatt und Bowen 2007). Thus, with increasing polymer content at high solid volume fractions, the overall elastic deformation  at the shear surface of particles exceeds rigid particle surface interactions, and elastic retention and shear-thickening rheological behavior occur. Strong interparticle forces between densely packed particles can accommodate greater stresses than the polymers, but might fail at much lower deformation if no elastic forces are incorporated. This can explain the most brittle failure for the test series 0.45-300, as the particle distances are high due to enough excess water and additional steric repulsion due to PCE. The lowest yield index and thus least rigidity is observed for the highest solid    volume fraction, as the particle network is strong and retains elastic energy due to the polymers.

Non-linear viscoelastic analysis
Following the macroscopic viscoelasticity analysis of the oscillatory first harmonics, non-linear intra-cyclic analysis characterizes the viscoelastic material properties within one load cycle and with increasing load more sophistically. Qualitative LB curves and quantitative Fourier and Chebyshev coefficients were analyzed for each test series. For each load Fourier coefficients and Chebyshev parameters acc. Equations (9), (15), and (16) describe the ellipse shapes and their corresponding material behavior mathematically from low to large deformations. LB-curves for each test series are given for three characteristic points within the linear viscoelastic (LVE), transient, and non-linear regime, respectively, in Fig. 16, with consecutive quantitative description in Table 5. For each flowability value and solid volume fraction, normalized LB curves are plotted in the linear viscoelastic (red), transient (blue), and non-linear (green) regime. LVE parameters were analyzed shortly below the linearity limit, as indicated Fig. 15b, to reduce the effect of noise (see raw data handling, Fig. 9). For test series with the lowest solid volume fraction, i.e., 0.45-200 and 0.45-250, the calculation of Fourier coefficients still led to inhomogeneities, which is why this paste composition can be seen as a boundary for nonlinear LAOS viscoelastic analysis. The transient material behavior above the linearity limit was investigated between L and F . Non-linear material analysis was performed above the gel point. For the sake of comparability, the deformation at the amplitude of 0 = 80% was analyzed for all test series. Fourier coefficients were calculated up to the 7th order for all calculations and plots. To understand viscoelastic characterization by LB curvatures, in Fig. 16, all LB curves are plotted dimensionless normalized by the maximum value (i.e., strain, stress). Thus, solely the curvature shape, not the surface area, should be compared within one test series. Generally, test series with dense solid volume fractions possess more solid-like behavior with narrow ellipses. Low solid volume fractions possess wider ellipses (and thus more viscous material behavior). With increasing amplitudes, viscosity increases. All test series show intra-cyclic nonlinear material behavior. As early as in the linear viscoelastic regime, slight deviations from the elliptical shape are visible, which increases for the transient and non-linear regime. The deviation from ellipsoid shape is more pronounced with increasing solid volume fraction. For all 0.45 test series, the deviation is quite low, whereas the 0.55 shows a strong increase in stress towards the maximum amplitude in the transient regime (backward-tilted shape, depicting shearthickening behavior) and a decrease of stress towards the maximum amplitude in the non-linear regime (forwardtilted, depicting shear thinning behavior). Finally, the LB curves are clearly deformed, and the viscous material properties predominate.
Relative intensity I 3∕1 and Chebyshev parameters S, T, e 3 and v 3 are given in Table 5 for the transient and nonlinear regime. I 3∕1 increases with increasing ϕ and increasing flowability. In the linear viscoelastic regime, neither Fourier coefficients nor Chebyshev polynomials are meaningful, and parameter analysis acc. Table 4 is sufficient. As elastic and viscous moduli stay more or less constant in the linear viscoelastic regime, Chebyshev parameters approach the value 0. The general ellipsoid shape is similar for the all test series -however, absolute values and the inter-cyclic evolution change greatly, which is adequately described through the Chebyshev constants. Generally, beyond linear viscoelasticity, the elastic shear stiffening index S is positive and the viscous shear thickening index T is negative. e 3 is positive (elastic shear-stiffening), and v 3 is negative (viscous shearthinning). Within a ϕ series, e 3 decreases with increasing flowability value; v 3 increases towards 0. The Chebyshev parameters increase positively with increasing slump flow value. An exception is the series 0.45-300. As the paste is extremely fluid, secondary effects such as sedimentation or increased wall slip might lead to contrary results, which should be investigated in further analysis.
As LB analysis indicates, Chebyshev parameters describe non-linear rheology in a sophisticated way. Thus, S,T , e 3 , and v 3 are plotted for each load cycle in Fig. 17. Figure 17a shows the inter-cyclic evolution of the strain stiffening ratio S , Fig. 17b shows the inter-cyclic evolution of the shear thickening ratio T , Fig. 17c shows the intercyclic evolution of the strain stiffening ratio e 3 , and Fig. 17d shows the intercyclic evolution of the shear thickening ratio v 3 . For clear presentation of material properties, all graphics are given in log-lin scale for the load range from 0 = 10 −1 − 10 2 . Below γ 0 = 10 −1 , Chebyshev parameters are close to 0 due to linear viscoelasticity.   In the transient regime, all solid volume fractions show an increase of the elastic stress response in the maximum deformation range, which corresponds to a positively increasing elastic shear stiffening ratio S plotted in Fig. 17a. With an increasing solid volume fraction, the positive increase of S is more pronounced and takes until a larger 0 , which is larger for higher ϕ values. The viscous shear thickening ratio T is consistently negative and drops significantly further with increasing . Combining the LB curves with the quantitative Chebyshev parameters, it can be seen that beyond the gel point F , the slope of the elastic shear stress drops and the LB curves change their shape significantly, as they deviate from the ellipse. The indices S and T are negative beyond that region. Thus, the viscous properties dominate and the material flows. Intracyclic elastic shear stiffening and viscous shear thinning description clearly illustrates the viscoelasticity in the yielding zone between the linearity limit and the gel point. In the nonlinear regime, all elastic stress curves are sharply decreasing near the maximum amplitude, which corresponds to negative S and T values. The increase in the width of the elastic LB curve over the measurement process shows the reduction of the ratio of storage modulus and loss modulus inter-cyclically. Intra-cyclically, the material behavior changes to elastic shear softening ( S decreases below 0) and viscous shear thinning ( T increases its magnitude towards 0). This results in an increasing dominance of viscous effects and an inter-cyclic increase in flowability. A possible explanation for this observation is that the deformation applied above the gel point exceeds the interparticle interactions to such an extent that the "yield point" of all elastic components of the suspension is exceeded and the bonds of the particle network are irreversibly destroyed. A permanent deflection of the particles relative to each other occurs, and the suspension flows. In combination with the analysis of first-order moduli G ′ and G ′′ , the superplasticizer-dependent reduction of stiffness and rigidity of cement pastes is described.

Viscoelastoplastic modeling
The presented analysis methods help to describe large-range rheological flow behavior of different cementitious suspensions and reveal the strong paste composition effect on viscoelastoplastic properties. Dependent on the solid volume fraction and polymeric content, fundamentally different Fig. 17 Non-linear material analysis through Chebyshev coefficients for all test series: a strain stiffening ratio S, b shear thickening ratio T, c strain stiffening ratio e 3 , and d shear thickening ratio v 3 flowability (viscoplastic) and deformability (viscoelastic) properties appear. Dynamic rheological analysis of the test series illustrates clear viscoplastic parameters for macroscopic shear-thinning or shear-thickening behavior and viscosity bifurcation; oscillatory rheological measurements are able to depict more elaborate transient inter-and intracyclic material behavior and viscoelastic parameters. While the first-order oscillatory parameters , * , and the yield index can be used to characterize the overall viscoelastic behavior, non-linear analysis gives insight into intracyclic shear-thinning or shear-thickening material behavior and thus knowledge especially regarding the effect of chemical admixtures. Figure 18 illustrates macroscopic viscoplastic and viscoelastic results for all test series. Paste composition effects are illustrated in Fig. 18a−c: With increasing , the non-Newtonian index increases independently of polymeric content. With increasing PCE dosage, both viscoplastic yield stress 0,H−B and viscoelastic complex viscosity * decrease. Not illustrated here, but analyzed before, is the effect of polymeric content on further rheological parameters such as zero and infinite shear-rate viscosity and critical yield stress. Figure 18d and e show the combination of viscoplastic non-Newtonian index n and viscoelastic parameters phase shift angle and yield index (as description of viscoelastic range). The results clearly indicate an increasing viscoelasticity with increasing non-Newtonian flow. Figure 18f maps the cementitious pastes regarding their viscoelastoplastic range as indicated in Fig. 5b. The results reveal that, if oscillatory and dynamic rheometric characterization methods are combined, they enable a material characterization both for viscoplastic steady-state flow and for admixture-dependent viscoelastic rheology.

Conclusions
The experimental investigations have shown that cementitious paste deformability and flowability are paste composition dependent with varying viscoelastoplastic behavior. Pastes with adjusted flowability values may behave similarly in one flow range but differ completely below or beyond. Therefore, the same rheological constitutive model and mathematical parameter calculation might fit one test series but falsify the rheological characterization of another series. This results in different appropriate rheological analysis methods. Particularly in case of viscosity bifurcation (and thus structural build-up at low shear), the rheological characterization of cementitious building materials through (traditional) dynamic flow curves is of limited use: Secondary flow effects during rheometric measurements can lead to unknown flow fields and thus errors during raw data conversion. Furthermore, mathematical regression methods at the boundaries of flow (no and low or very high shear ranges) incorporate potential calculation errors. Steady-state rheological evaluation in an explicit flow equation is only physically justified for a true viscoplastic (e.g., Bingham or Herschel-Bulkley) fluid or for a clearly identified load range. These boundary material analyses can be scoped with oscillatory analysis methods that clearly identify viscoelastic parameters in the linear viscoelastic range. In recent research, dynamic and oscillatory parameters, especially the yield stress, have generally been compared in order to shift rheological parameters. However, we consider that it is more beneficial to combine the methods and both identify viscoelastoplastic material behavior over a large range of shear and take advantage of sophisticated non-linear measurement techniques in the transient material states. The main research findings are.
1. Dynamic phenomenological analysis of complex non-Newtonian cementitious pastes requires a distinct choice of the mathematical regression model. Depending on the point of interest, large-range flow analysis or boundary rheological parameters can be identified. However, transient or discontinuous material properties during shear or the choice of an insufficient rheological model can lead to tremendously wrong material characterization. 2. A main benefit of distinct mathematical viscosity modeling is the clear analysis of zero and infinite shear-rate viscosity, which are mandatory input parameters for the starting point of flow calculation, e.g., flow simulations using computational fluid dynamics (CFD). A combination of viscoplastic simple shear and viscoelastic oscillatory methods can provide reliable material parameters for both no and high shear rates. 3. A phenomenological unifying approach of combined dynamic and oscillatory linear and non-linear rheological analysis provides more information regarding rheological flow behavior of cementitious suspensions over a broad range. A combination of the presented methods does not lead to comparative (as suggested in previous research) but to complementary material parameters that help to fingerprint the viscoelastic or even elastoviscoplastic deformation and flow. 4. The rheological property map gives an indication about the effect of solid volume fraction and flowability index on wide-range rheological flow behavior. From a phenomenological point of view, this approach helps to characterize cementitious building materials according to their most appropriate rheometric investigation and constitutive equation. Nevertheless, the bifurcation points, i.e., the dynamic critical shear ̇c rit , elastic linearity limit L , and gel point F depend on interparticle forces of cementitious suspensions, which themselves are a function of chemical and physical composition of the colloidal particle network. A comprehensive modeling approach for cementitious suspensions would need to add a microstructure dimension that includes further microstructural parameters.
Open questions mean that further research is essential: 1. Broadening of the rheological property map over the material scale (investigation of effect of various other additives or admixtures) and over the material analysis scale (phenomenological rheological description) 2. Broadening of the rheological property map from a purely phenomenological to a microstructure-implemented rheological description. A more general approach to classifying a broad range of cementitious building materials and linking flow behavior to mixture composition, flow history, and processing is needed for the prediction of complex concrete flow and holistic material characterization. The presented research serves as basis for further microstructural and kinetic analysis. 3. Applicability of the chosen rheological models in further flow analysis, e.g., through linking to numerical fluid dynamics, needs to be investigated.
In summary, the research presents a holistic material characterization approach of complex non-Newtonian cementitious building materials. Established and state-ofthe-art analysis methods were analyzed and critically discussed. It was found that especially for the case of complex non-Newtonian behavior, more than standard mathematical models need to be chosen for material characterization. The additional viscoelastic analysis described rheology in the low-and no shear rate zone, transient material behavior, and inter-cyclic non-linearity in a more elaborated way than common analysis methods. Finally, a rheological map with different rheological models and hence calculation methods was presented that will to serve as aid for the prospective analysis of complex cementitious building materials.