A simulation-based approach to evaluate objective material parameters from concrete rheometer measurements

This work is licensed under the Creative Commons Attribution 4.0 License Appl. Rheol. 2019; 29 (1):130–140 Research Article Florian Gerland*, Alexander Wetzel, Thomas Schomberg, Olaf Wünsch, and Bernhard Middendorf A simulation-based approach to evaluate objective material parameters from concrete rheometer measurements https://doi.org/10.1515/arh-2019-0012 Received Jul 23, 2019; accepted Sep 09, 2019 Abstract: Modern concretes such as ultra-high performance concrete (UHPC) show excellent strength properties combinedwith favorable flowproperties.However, the flow properties depend strongly on process parameters during production (temperature, humidity etc.), but also change sensitively even with slight variations in the mixture. In order to ensure desired processing of the fluidlike material and consistent process quality, the flow properties of the concrete must be evaluated quantitatively and objectively. The usual evaluation of measurements from concrete rheometers, for example of the ball probe system type, does not allow the direct determination of the objective material parameters yield stress and plastic viscosity of the sample. We developed a simulation-based method for the evaluation of rheometric measurements of fine grained high performance concretes like self-compacting concrete (SCC) and UHPC. The method is based on a dimensional analysis for ball measuring systems. Through numerical parameter studies we were able to describe the identified relationship between measuring quantities and material parameters quantitatively for two devices of this type. The evaluation method is based on the Bingham model. With this method it is possible to measure both the yield stress and the plastic viscosity of the fresh sample simultaneously. Device independence of the evaluation process is proven and an application to fiber-reinforced UHPC is presented.


Introduction
With the increasing economic relevance of modern concretes such as ultra-high performance concrete (UHPC) or shotcrete, the demand for precise characterization of the rheological properties of fresh concrete mixtures is getting more attention. This is necessary to ensure consistent process quality even if the rheology of the fresh concrete undergoes changes. Such changes occur as a result of batch changes, natural fluctuations in raw material properties, scalability of the mixture quantities produced or environmental influences such as fluctuating temperature or humidity [1].
The measurement of the rheological behavior of cement based systems is part of standard techniques. For the evaluation of the consistency of the concrete there are several standard test setups (DIN EN 12350-5, DIN EN 1015-3, DIN EN 12350-12) leading to the characterization of the flow behavior by measuring the spread diameter of a certain amount of concrete after a certain time and is then classified from C0-C3 and F1-F6 with increasing flowability (DIN EN 206-1/DIN 1045-2).
A weakness that is common to this classification is its inadequate objective measurability. The flow behavior can only be roughly predicted. In addition, this classification does not allow the fresh concrete properties to be characterized quantitatively to the extent that the flow properties can be used for flow simulations, for example of mold filling processes. Furthermore, experiments such as the slump flow measure can only provide limited comparability due to measurement inaccuracies and operator influences.
Rotational rheometers (e.g. cylinder test specimens according to DIN EN ISO 3219) perform a constant shear rate resulting in a constant shear stress in the gap, the relationship of which can be directly deduced from the viscosity of the sample material by a simple slit flow. Due to the very narrow gap of typically 1 mm, these devices are unsuitable for concretes, especially when fiber-reinforcement is applied. Modifying the geometry to better suit the requirements for concrete-like suspensions can be done but leads relevant drawbacks, e.g. that an analytical description is not possible anymore, as investigated by Ferraris et al. [2] and Feys et al. [3]. Alternative approaches to avoid the difficulties associated with the small gap size in typical rotational rheometers have been presented by several authors. The underlying idea is to use device geometries that are large compared to the granular constituents of the concrete, so that the movement of the concrete can be understood as a continuum, even though they do not produce a simple flow state. High performance mixer [4], rotation rheometer with mixer configurations [5] or ball probe rheometers for coarse grained concrete are prominent examples. A ball probe rheometer has been proposed by Fleischmann et al. [6,7] for measurements on selfcompacting concretes. These kind of devices measure the relative motion (i.e. velocity) between the gauging head and the sample container as well as the fluid force acting on the gauging head similar to rotational rheometers. The principal construction is illustrated in Figure 1. This type of instrument is comparable to falling-ball systems. These instruments are designed to determine Bingham material parameters and move a ball vertically through the fluid at a constant velocity, driven by gravity [8]. Ball probe rheometers also move a ball through the fluid. However, the velocity can be varied by changing the power of the motor. However, due to the more complex flow conditions within such a device, an analytical solution and therefore a direct return of the objective material properties is not known. Instead, the analysis of the measured values based on simplifications like a stratified flow of a Bingham plastic and is done by derivation of representative quantities. A main challenge is the definition of a representative shear rate, because the shear rate on the surface of the sphere differs locally and depends on the rotation velocity of the sphere. Fleischmann's method predicts a linear affine relationship of the torque (expressed by a representative shear stress) from the rotational velocity (expressed by a representative strain rate) for the Bingham plastic. The slope of the curve is determined solely by the plastic viscosity, while the axis intercept is determined solely by the yield stress. In the measurement series presented by the author, at the same time, a nonlinearity of the measurement curves in the range of low shear rates is clearly recognizable in which the model underestimates the measured values. Different measuring geometries (e.g. different ball diameters or cylinders) led to different absolute values of the material parameters, which proves that these cannot be reliably determined with the evaluation method. A transferability of the measured values between devices with different geometries is thus possible to a very limited extent.
In order to gain a deeper insight into the most often complex flow in concrete rheometers of different types, various authors have applied numerical flow simulations. Oliva et al. [9] investigated the flow in a vane type concrete rheometer to improve its geometry. Zhu et al. [10] investigated the flow in a device of the same design using a twodimensional Smoothed Particle Hydrodynamics method and were able to show that the characteristics of the flow are influenced by the dimensionless Bingham number. Le et al. [11] used numerical simulations to develop an improved method for determining material parameters based on parameter studies of the torque-speed relation of devices with Tattersall geometry. The authors investigated samples with known Herschel-Bulkey material parameters and were able to prove that both slope and axial section of the measurement curves are dependent on plastic viscosity, yield stress and yield index. The authors tried to fit a model equation in the parameter studies, which indicates the relationship of the torque-speed relation of the material parameters. However, they reported weaknesses of the model in some of the investigated substances, which they attributed to an insufficient approximation approach. Wallevik et al. [12] also report in a review on the complex de-pendency of the torque-speed relation on the material parameters.
In this paper, we provide a method to overcome the reported issues by a systematic analysis of the dependencies between the physical quantities involved and by the support of numerical modeling of the flow in the concrete rheometer. We develop an evaluation methodology to determine objective material parameters such as plastic viscosity and yield stress from the output quantities of concrete rheometers with complex flows like the ball probe type without a necessity to define any representative shear rate. We prove the device-independence by comparison of the results from two different devices, eBT2 and Viskomat NT (both from Schleibinger Testing Systems Teubert u. Greim GmbH).

Dimensional Analysis
For concrete rheometers based on the principle of the balltype measuring system, a dimensional analysis is carried out in order to reveal fundamental relations. A torque is measured as a function of a predefined rotational speed.  geometric quantities could be added formally. However, the ones mentioned here are the essential ones, that sufficiently characterize the flow problem. The presented dependence provides the dimensionless numbers presented in Eq. (1) to (7), namely the dimensionless torque Π 1 , Reynolds number Π 2 = Re, Bingham number Π 3 = Bm, as well as the geometric dimensionless numbers Π 4 , . . . , Π 7 .
The relationship between the dimensionless numbers is , whereby the function itself is unknown and will be approximated by the results of a numerical parameter study.

Numerical model
The concrete is modeled as continuous, homogeneous fluid. For the solution, the finite volume method is applied within the framework of OpenFOAM ® v1712. Due to the incompressibility of the fresh concrete, the mass balance is formulated with the fluid's velocity v We assume that the concrete obeys the Bingham flow model, i.e. that it has a yield point and a constant plastic viscosity. In principle, it would also be possible to use more complex flow models involving shear thinning (e.g. Herschel-Bulkley). However, due to the low rotational speed in the concrete rheometer and the resulting low maximum shear rate in the flow, no shear thinning can be expected and the additional material parameter could not be determined from the measurements. The low rate of change of the rotational speed allows a quasi-steady state treatment of the flow¹, so that the momentum balance results in: where T is the stress tensor and p represents the pressure.
The following relationship applies to the stress tensor with the shear rate dependent viscosity η(˙), the strain rate tensor and the shear rate˙: The Bingham model readṡ The viscosity function is defined as which results in for the Bingham material in case of |τ| > τy. This flow function results in diverging viscosity for low shear rates liṁ →0 η Bm (˙) = ∞, which makes it unsuitable for the numerical realization. Therefore, instead the Papanastasiou regularization is applied [13,14] η with the regularization parameter α, which limits the viscosity to liṁ  range of material parameters investigated, it was found that α ≥ 3000 s is sufficient as convergence of the recorded torques is reached in all cases. A convergence plot is presented in Figure 3 for both, Viskomat NT and eBT2. The diagram shows the evaluated torque depending on the regularization parameter normalized by the evaluated torque in case α = 300000 s. The study was performed for the highest Bingham numbers appearing in the numerical study by selecting the highest yield stress and the lowest plastic viscosity as well as the lowest rotational speed. It must be noted, that for vanishing rotational velocity, the Bingham number strives for infinity and no shearing of the fluid appears at all. For that case, the Papanastasiou regularization parameter must be infinite in order to represent the Bingham material appropriately. However, the measuring devices used in this study are limited in the lowest rotational speed at which they can work. Therefore, we selected the lowest rotational velocity investigated numerically to 0.002 rad/s for both devices. This does not negatively influence the deviation of neither yield stress nor plastic viscosity, as the yield stress mainly influences the height of the overall curve (i.e. not only in the region of very low rotational speed but as well in the higher range) and the plastic viscosity mainly influences the slope of the M-Ω curve.
The influence of the concrete's free surface is negligible with regard to the torque measurement within the examined parameter range. This was proven by numerical two-phase simulations with free surface compared to single-phase simulations. This is as well supported by the topology observed in the experiments. A significant deviation of the topology from the flat liquid level can be observed neither with silicone oils nor with concrete.
For the eBT2, the sample container rests while the sensor block rotates. Therefore, the no-slip condition v = 0 applies to all walls of the sample container, whereas the varying velocity on the sphere with the distance from the axis of rotation r is v = Ωw×r, where w is the unit vector of the rotation axis. The solution was performed in the noninertial rotating reference frame and the fictitious forces were taken into account². In contrast, for the Viskomat NT the sample container rotates and the solution has been carried out in the non-rotating inertial system. Therefore, the no-slip condition for the wall velocity is v = Ωw × r. The velocity of the sphere is v = 0. At the interface between concrete and air, zero shear (∂v/∂n = 0) was applied on both devices assuming negligible topology changes (v · n = 0) as stated above.
The 3d-geometry of the eBT2 was discretized by a total of 303200 cells, of which 95% are hexahedral. In order to adequately resolve the stresses at the surface of the spheres, these are each discretized with about 18600 cells, see Figure 4. The Viskomat NT was meshed with a total of 448556 cells (95% hexahedral). There are 2144 cells on the surface of the sphere. The amount of cells on the surface of the sphere is smaller here, because the volume of the sphere relative to the sample volume is much smaller than for eBT2 and the gap between the sphere and the wall is relatively large, so that smaller gradients appear. For both geometries convergence was proven in a mesh independence study.
The parameter study of the Bingham parameters was carried out in the same range for both devices. The range of material parameters examined was based on typical fiber-reinforced UHPCs: the plastic viscosity was varied between 5 and 500 Pa s, the yield stress between 0 and 50 Pa. The rotational speed was varied between 0.002 and 0.2 rad/s for eBT2 and 0.002 and 0.12 rad/s for Viskomat NT according to the measuring profiles of the instruments. Altogether, the number of data points (M; Ω; τy; µ) for each sensor of the eBT2 was 378 resp. 294 data points for the Viskomat NT.

Experiments
The measurements were carried out with two different instruments: eBT2 and Viskomat NT, both by Schleibinger. Both instruments were equipped with spherical measuring bodies and have a design which corresponds to the characteristics described in Section 2.1 and illustrated in Figure 1. The characteristic dimensions of both instruments are presented in Table 1. The eBT2 is equipped with two measuring bodies (MK105 and MK145) at different distances from the central axis, which are offset by 180 ∘ . In this instrument, the specimens rotate while the sample container rests. On the Viskomat NT, however, only one sample body is mounted, which rests while the sample container rotates around its axis.
For practical determination of the flow curve, a symmetrical ramp profile of the rotational speed with linear acceleration and deceleration was applied. The measured data was each recorded in a single rotation of the sensors. For validation of the numerical model two silicone oils were used, which behave Newtonian without yield stress: PDMS³ 12500 and PDMS 100000, with a density of 970 kg m −3 . We measured the dynamic viscosity of both oils in a rotational rheometer AR-G2 by TA Instruments. For PDMS 12500 this is 12.3 Pa s and for PDMS 100000 97 Pa s, each without shear thinning in the range of interest.
In a next step the behavior of a fine-grained fresh UHPC formulation was tested purely as well as after addition of fibers in rod-and ring-shape. The concrete formulation is based on the standard formulation investigated within the German Priority Programme 1182 [15]. The ratio of the solids among each other, the solids concentration and the proportion of liquid components in the mixture were kept constant. In case of fiber reinforcement, highstrength micro-steel fibers were used. The concrete was produced using an intensive mixer R05T (Eirich), which has a maximum mixing volume of 40 L. The steel fibers were added with an inclined vibrating chute to the mixture which ensured an optimum separation and addition of the fibers. Evidently, in no mix fiber agglomerations could be determined, which suggests a very homogeneous mixture. Two batches were produced and measured simultaneously in both devices in order to be able to prove the device-independence of the measurement results as well as the presented evaluation procedure. The mixes were tested using a linear acceleration and deceleration with a maximum rotational speed of 0.18 rad/s in eBT2 and 0.11 rad/s in Viskomat NT, both of which are within the simulated range described above. The measurements were performed approximately 10 min after the end of the mix. In order to avoid effects such as thixotropy and inhomogenization, the sample material had been pre-sheared before starting the measurements. However, since in some experiments a hysteresis of the torque-velocity relation is still recognizable, for all concrete probes only the deceleration curve was evaluated.

PDMS (Newtonian Behavior)
In order to validate the numerical model, the experimental measurements of the two silicone oils are compared to the corresponding simulation results. The measurements were carried out on two days with two repetitions each. The low scatter between the measurement series shows a very good reproducibility. The results of the validation are shown in Figures 5 and 6. A good agreement of the simulations with the measurement results can be seen for both instruments.

Functional Relationship
The aim is to deduce the Bingham material parameters from the M-Ω measurement curves. For this purpose, a functional relationship of the form M = f (Ω, µ, τy) is to be found, which provides the Bingham material parameters of the sample by parameter fitting to the measurement data (M, Ω). The basis to find this relationship is founded on the dimensional analysis in Section 2.1 and the numerical study. In order to determine the functional relationship between the dimensionless torque and the Bingham number Π 1 = f (Bm), the data points determined in the numerical study are presented in the Π 1 -Bm diagram in Figure 7. The data points of all measuring bodies are each described by a continuous curve, which can be approximated by functions of the form with sufficient accuracy. This function type has been selected due to its ability to represent the whole measuring range with as few coefficients as possible. The coefficient a 0 describes the pure viscous part of the material characteristics. It can be determined directly by evaluation on Newtonian material by evaluating the constant expression Π 1 (Bm = 0 for Newtonian material). a 1 and a 2 describe the slope respectively the beginning of the slope of the function. The coefficients a 0 to a 2 depend on the geometry of the device and thus are functions of the dimensionless numbers Π 4 to Π 7 , wherein Π 7 is the significant influencing factor, because the ball speed is decisively influenced by the lever. This is also indicated by the fact that the coefficients for the two measuring bodies of the eBT2 are not identical, although both differ significantly only in Π 7 . However, the geometric dependencies are not further investigated since the coefficients a 0 to a 2 can be determined directly from the parameter studies. The dimensionless torque in the operative rotational speed range is no function of the Reynolds number (Re ≪ 1, creeping flow). The coefficients identified by parameter fitting can be found in Table 2. Note that the highest Bingham numbers appearing in this study result for very low rotational speeds. This is due to the physically reasonable selection of material parameter range. However, in the range of very low rotational speed, the measured torque is low compared to the medium and higher rotational speed range independent of the material parameters. The rise of the dimensionless torque visible in Figure 7 is caused by its inverse proportionality to the rotational speed (cf. Eq. (1)). Therefore, the fitting of the coefficients of Eq. (17) to the simulation data has been optimized for the lower to medium Bingham numbers with respect to the curvature at the beginning rise of the graph and to the slope thereafter.
By reinserting the dimensioned variables into the correlation function Eq. (17), the following relationship is obtained which is presented exemplary in Figure 8 along with the simulation data for the Viskomat NT. Note that for Newtonian fluids (τy = 0 Pa s) Eq. (18) reveals a linear proportionality between the rotational velocity of the sphere and the torque. This is consistent with the well-known Stokes relationship for a sphere in a free-streaming Newtonian fluid F = 3πηDU∞ [16]. For illustration, the resulting force of a Newtonian fluid acting on the ball probe with Ωl correlating to U∞ and Newtonian viscosity η equals µ, is F = M/l = a 0 µD. However, the proportionality coefficients are higher than 3π ≈ 9.42 for all devices due to the presence of the container walls as well as the circular path, which both cannot be neglected.
On the simulation data presented in Figure 8, we can clearly recognize a nonlinearity of the curves for non-Newtonian fluids, which intensifies with increasing yield stress (i.e. with lowering Bingham number). For that reason, the model coefficient a 2 will always be less than 1.0.   Table 2. The comparison shows selected material parameter combinations (i.e. Bingham materials) for the whole range of rotational velocity investigated. (eBT2, MK145) In case of a 2 = 1.0 Eq. (18) would be an affine-linear relation between torque and rotational velocity starting from an offset-value linearly related to τy, which is neither valid for our investigation nor for those of other authors [11,12].
As Figure 8 demonstrates, the selected approximation function Eq. (18) with only three instrument-constants allows a sufficient representation of the torque-velocity relation for the entire measuring range. However, it must be noted that the model function for vanishing rotational velocities has the limit value lim Ω→0 M(Ω) = 0, which does not correspond to the expectation. Rather, we would expect a constant, positive and solely yield stress-dependent limit value. This circumstance, nevertheless, is of no concern for the evaluation of measurement data, since the measuring instruments do not record any data for infinitely small rotational velocities. A more complex model would therefore not offer any increased benefit and, in the worst case, would involve the risk of overfitting.
In the torque-speed curves, a nonlinearity occurring with increasing yield stress is clearly recognizable. It has been excluded that this is caused by the regularization of the Bingham material model by checking the convergence for large α, as shown in Figure 3. The model does not include shear thinning (e.g. the Herschel-Bulkley law) or thixotropy, thus the shape of the curves must have other reasons. Rather, the cause is to be found in the size of the unsheared fluid portion varying with the shear rate. This domain is defined by η Bm (˙)˙≤ τy, which means the portion of fluid, where the shear stress remains below the yield stress. It is exemplary shown in Figure 9 for τy = 50 Pa and µ = 50 Pa s for Ω = 0.01 rad/s and Ω = 0.2 rad/s. The lower the unsheared fluid portion is in the overall volume, the weaker is the nonlinearity of the torque-speed relation. The interface between sheared and unsheared fluid domain behaves similar to a wall for the sheared fluid. With rising rotational velocity, i.e. with rising shear stress, the interface moves away from the probes. This is similar to an increase of the wall distance and results in a lower rise of the shear rates at the probe. The increase of the shear rate drops down with the decrease of unsheared fluid volume until most of the fluid is sheared. From that point on, the increase of shear rate at the probe is constant and therefore the relationship between rotational velocity and torque is linear proportional.
Eq. (18) can thus be used to determine the material parameters µ and τy from the recorded M-Ω curves by parameter fitting for unknown material samples.

Application to UHPC
In Figure 10 the measurement results obtained in eBT2 for the total of four measurements from two batches of the fiber-free UHPC are presented.
Other than with the silicone oils, recognizable deviations between the measurements are noticeable, which indicates a slightly lower repeatability. This is due to slight variations in the formulations of both batches, the strongly time-dependent change in the flow properties associated with the high reactivity of the concrete, and generally the limited validity of the continuous model for fresh concrete. The curve shown in those figures depicts the model approximation of the measured data for the material parameters µ = 86 Pa s and τy = 5 Pa, which represents an average of all measurements. These parameters represent the measurement results of both sensors equally well. Concrete of the same batch was measured simultaneously in Viskomat NT. The measurement results are displayed in Figure 11.
For the fiber free UHPC, fitting Eq. (18) to the measurement data of Viskomat NT results in µ = 98 Pa s and τy = 4 Pa. This deviation is in the order of magnitude of the scatter of the measurements (cf. Figures 10 and 11) and thus demonstrates the device independence of the measurement results.  Furthermore, Figure 11 contains two additional measurements, each of the same basic mix, but with additional 2.0 vol-% of steel fibers. In one batch, common rodshaped fibers with a length of 10 mm and a diameter of 0.4 mm were added. In the other batch, fibers of the same dimension were added, which have been reshaped into rings with a diameter of approximately 3.5 mm. It should be noted that such geometry has currently only academic value but can be understood as an example when dealing with innovative fiber geometries in the future. However, we can demonstrate that this deformation can have a drastic influence on the flow properties of the concrete. It can be seen that the approximation function with the material parameters τy = 16 Pa, µ = 357 Pa s for the concrete with rod-shaped fibers and τy = 11 Pa, µ = 228 Pa s for the concrete with ring-shaped fibers fits very well. Also the slight curvature, which is recognizable in the measurement data at low rotational velocities, is captured by the model. The increase in plastic viscosity and yield stress due to the addition of fibers is shown in Figure 12. It is evident that ring-shaped fibers in our investigations can reduce both the increase in yield stress and in plastic viscosity by about one third compared to rod-shaped fibers.

Conclusion
We presented a method which completely takes into account the complex flow behavior in concrete rheometers of the ball measuring system type. We assumed Bingham material behavior, which can represent our measurements well. For the investigated material and the investigated load range (shear rate range) it is therefore a suitable material model. By measuring the concrete in two rheometers of different geometries, the reliable transferability of the measurement results between measuring instruments is demonstrated and the objectivity of the evaluated material parameters is underlined. With the knowledge of objective parameters it is now possible to carry out flow simulations e.g. for mold filling processes. These can be used to define quantitative limits of the parameters, which must be reached by a batch not to be rejected. Material parameters can now be measured quantitatively with the aid of the devices immediately before the mold filling is carried out.
The simulation of the flow shows a nonlinear torquevelocity relation depending on the influence of the yield stress due to the decreasing size of the non-sheared fluid region with the rotational velocity. Models that simplify the flow considerably (e.g. [7]) cannot take this into account, which leads to incorrect material parameters for concretes with pronounced yield stress.
The method can also be used for other sample body geometries (cylinders, mortar paddles) or other rheometer types with flow conditions that are not described analytically or only approximatively (vane or similar), but also in situ experiments by evaluating the torque in the mixer are conceivable, so that a subsequent measurement would not be necessary. The assumed material model can be adapted to the expected material behavior and is not limited to Bingham plastic. By selecting a more complex material model (e.g. Herschel-Bulkley), one should however consider that the highest shear rates occurring in this investigation were in the maximum order of 10 s −1 , depending on the material parameters. If the material shows effects like shear thinning at higher shear rates, these can only be measured in a concrete rheometer of the type discussed here, if correspondingly higher shear rates predominate. This could be achieved by smaller gaps (unfavorable since the continuum hypothesis could be violated) or by higher rotational velocities.
For fiber-reinforced concrete, flow-driven fiber orientation can occur and thus influence the flow properties of the concrete itself [17]. The method presented can be extended by a suitable fiber orientation model as presented in [18].
In the acceleration curves of concrete and mortar measurements, a super elevation of the torque is often recognizable, which can be attributed to thixotropic behavior. By considering a thixotropy model like [19], this can also be considered in the numerical model. By evaluating the entire experiment (acceleration and deceleration characteristic), the thixotropic behavior could be additionally quantified and its influence on the mold filling simulated.