Comptes Rendus Mécanique

Modeling of the vane test using a power-law ﬂuid and model order reduction techniques: application to the identiﬁcation of cement paste properties


Introduction
Concrete casting is an active research topic. Faster and more reliable casting processes are sought on a daily basis, without compromising the final structural properties. However, the material's behavior and properties are not mastered yet, and are currently the subject of many ongoing research works, both experimentally and theoretically [1,2]. Because of the material complexity as well as the 3D flow involved in concrete casting, simulation becomes necessary to understand the physics involved in the process. Some works aimed at modeling the fresh concrete and cement paste behavior using finite element methods [3] with an elastic-viscoplastic approach, while others resorted to data-driven modeling using neural networks and artificial intelligence [4,5].
In this work, we target the identification of material parameters, cement paste parameters for instance, using two different vane tests, to identify the flow consistency K and the flow behavior index n involved in a power law. Model reduction techniques are used to create a manifold of solutions and identify the missing properties from experimental data. The application shown here is performed for a cement paste, however, any power-law fluid can be characterized using the same method developed in this paper, by comparing the obtained experimental results to the model simulation as illustrated in this work.
To identify the material properties by comparing experimental results to simulations, we need to compute many direct solutions of the problem at each iteration of the optimization algorithm used to minimize the difference between the experimental results and the simulations, e.g. [6]. For instance, using a gradient-based simulation requires a minimum of two simulations per parameter at each iteration to identify the sensitivities of the material properties [7]. Therefore, using classical simulation techniques, simulation-based optimization algorithms become highly time consuming. The chosen material for this application is a fresh cement paste, a fluid with particle suspension which is commonly modeled using a non-Newtonian power law [8,9]. This work models the fresh cement paste as a fluid, without a yield stress, as opposed to previous approaches considering the paste as a solid subject to yielding. In fact, non-Newtonian powerlaw modeling is commonly used for fresh cement paste modeling and yields better results than the solid model counterparts [10][11][12][13][14][15]. Indeed, the latter hypothesis leads to predications contradicting experimental data [16]. This is discussed further in Section 4. This is where the use of model order reduction techniques becomes an appealing approach, as successfully used in various previous works [17]. In fact, one may wish to add all the material and process parameters as extra coordinates of the problem. For instance, in a vane test of a power-law fluid, the rotational velocity ω, the fluid flow consistency K and the flow behavior n can be considered as extra parameters of the problem. Such simulation will lead to a multidimensional (6D) problem, leading therefore to prohibitive computational time when using classical simulation techniques. However, using model order reduction techniques, one may solve the 6D problem by a sequence of resolutions in lower dimensionality spaces [18]. Therefore, we choose to use the Proper Generalized Decomposition (PGD) to solve the problem as a sequence of 3D × 1D solutions. Due to the high nonlinearities involved in the process, the use of the PGD is combined to a linearization process as well as the use of a reduced basis constructed "offline" and used "online" to identify the required properties. An "offline" solution can be computed on a cluster and contains all the possible variabilities of the required parameters. It therefore constitutes a "Modern Chart" or an abacus. The "online" exploration of the chart is performed fast, in real time, on light platforms such as tablets or smartphones [19]. Such an approach is possible since the "online" phase consists of only straightforward arithmetic operations to reconstruct the full 3D solution of the required quantity of interest of the problem.

Modeling of the vane tests
In this section, we detail the modeling of the two different vane tests required for the identification of the two material parameters K and n. We choose to model the vane test using a solid vane illustrated in Figure 1(a) and a slotted vane illustrated in Figure 1(b). Both vane models contain four blades each, the in-plane thickness of the vanes is 0.7 mm. The vanes rotate inside a cylinder of diameter d = 100 mm and filled to a height h = 115 mm. The outer dimensions of the blades are the same; one vane is slotted. The ratio of the vane total height to the vane total diameter is r = 2. The fact of having a slotted vane will change the torque on the vane, for the same rotation speed, leading to a system of two equations with two unknowns while identifying the properties as shown in Section 5.
The rotation speed is imposed on the axis of the vane. Since cement paste is a fluid with particle suspension characterized by a shear-thinning behavior, we model the fluid using a nonlinear power-law pseudoplastic behavior. The apparent viscosity η of the fluid therefore writes [17,20]: where K and n are the flow consistency and flow behavior index respectively, D eq is the equivalent strain rate tensor defined by: while D is the strain rate tensor defined by: v being the velocity of the fluid (u, v, w). For n = 1, the fluid becomes Newtonian, with constant viscosity independent of the velocity. For n < 1, the fluid exhibits shear-thinning behavior. In the simulated problem, we consider non-slip conditions at all the fluid-solid interaction surfaces, as well as a free surface on the top of the fluid domain. Therefore, the boundary conditions write:    v = 0 at z = 0 and at x 2 + y 2 = 0.05 m v = ω × r at the contact with the vane and rotating shaft σ · n = 0 at z = 0.12 m, where ω is the imposed constant rotational speed and r the in-plane radius measured from the zaxis: r = x 2 + y 2 . σ is the stress tensor and n the outward normal to the free surface. Considering the relatively low imposed rotational speed ω, the deformation of the solid vane is neglected. Replacing (1) in the steady-state conservation of momentum equation of the studied fluid, one may write: where P is the pressure in the fluid. In this work, we solve the problem using the penalty formulation. In fact, neglecting density changes, the conservation of mass writes: Considering λ small enough, the penalty formulation writes [21]: Replacing (7) in (5), the final governing equation of the problem writes: Considering is the strain rate vector defined by: and multiplying by a test function v * and integrating by part, the weak form of the problem becomes: where * is the strain rate vector of test function v * and C is the 3D constitutive matrix defined by [21]: All the surface integrations vanish in our problem since we only have non-slip conditions or free surfaces. To solve the weak form depicted in (10), we consider the material parameters K and n as extra coordinates of the problem, as well as the imposed rotational velocity ω appearing in the boundary conditions in (4). This leads to a 6D problem. Using the PGD, we attempt to solve the problem as a sequence of 3D × 1D problems where the 3D domain is the physical (x, y, z) domain, whereas the 1D domains correspond to K , ω and n. One may imagine to decompose the 3D physical domain into a sequence of 2D × 1D domains [22][23][24]. However, since this work focuses on the nonlinear aspect of the simulation performed using the PGD, we will keep the 3D physical domain as one domain and focus on the treatment of the nonlinearity involved in this problem.
The meshing of the 3D physical domain is performed using GMSH, an open-source meshing software available online [25]. The mesh of the solid vane test is made out of 44,235 tetrahedral linear elements, as illustrated in Figure 2, whereas the mesh of the slotted vane has 41,711 tetrahedral linear elements and is shown in Figure 3.
In Section 3, we illustrate the PGD formulation and the methodology used to solve the multidimensional problem.

PGD formulation and solution of the problem
In this section, we illustrate the construction of the solution of the multidimensional (6D) problem depicted in (10). Considering the problem to be highly nonlinear in v and n, we start by a linearization process to construct the PGD multidimensional solution [26][27][28]. Since the nonlinear term η depends on n, we cannot easily introduce n as an extra coordinate of the problem using the PGD simulation. In fact, each operator of the weak form of the problem must be written as a sum of products of operators in each dimension. For instance, if L is a linear operator, we should be able to write it as [26]: where A i j are linear operators, N a finite number, M is the number of coordinates in the problem and c j the problem coordinates. However, due to the presence of n as an exponent in the viscosity operator, such notation is impossible for a realistically small value of N , if n is one of the coordinates c j . Instead, we will solve a 5D PGD problem for different snapshots, each snapshot being a simulation for a given value of n. This amounts to solving a 5D problem ten times, for n = 0.1, 0.2, . . . , 1, which provides 10 snapshots. The quantities of interest (the torque on the system for instance) shall be computed by interpolation between the snapshots. Such an approach is successfully used in many model order reduction works to construct the reduced basis [29][30][31].

Solving the linear problem for n = 1
To solve the 5D problem, we start by a linearization process. First of all, we solve the problem for n = 1: the Newtonian fluid case. In this case, Equation (8) is reduced to: ∇ 2 designating the Laplacian operator. Since K is constant, Equation (13) develops in Cartesian coordinates to: C. R. Mécanique -2021, 349, n 3, 501-517 where u, v and w are the three Cartesian components of the velocity vector v. To solve (14) using the PGD, we suppose v can be written in a separated form as: which can also be written in a compact form as: We also define v as a series known until an order o with: and we define v o+1 by: In what follows, we try to compute the enrichment functions R, S and T . By replacing v o+1 into the weak form of the problem depicted in (10) and considering n = 1, one can reach the separated form of the problem. To solve this problem using the PGD, one must write the test function v * also in a separated form [32]. Thus, we assume the test function can be written in a separated form as: This assumption is successfully used in various publications addressing different physical problems [6,19]. By replacing (18) and (19) into the weak form of the problem (10) and using the boundary conditions (4) one reaches the final problem to solve. However, the resulting problem is nonlinear in R, S and T . To overcome this difficulty, we use a fixed point rank-one update greedy algorithm outlined as follows: (1) We consider S and T known, therefore S * = T * = 0. The resulting equation is a 3D problem where R(x, y, z) is the only unknown. Once (R, S , T ) converge, v o+1 is known and we may search for v o+2 . The iterations continue until the convergence of the residual of the weak form in (10).
To solve the problem in R, we are using the mesh depicted in Figures 2(b) and 3(b) for the solid vane test and the slotted vane test, respectively. To solve the problem in K and ω, we are using two 1D meshes. The first one meshes K ∈ [1, 5000] Pa·s and the second one meshes ω ∈ [0.001, 0.1] rad/s. The iterative PGD algorithm converges at N = 7 within a few minutes on a normal portable PC. The velocity profiles in both mixers are shown in Figure 4 for ω = 0.0505 rad/s and K = 1 Pa·s. Once the linear Newtonian problem is solved, we start solving the nonlinear problem as explained in Section 3.2.
While solving for R, and since we are using the penalty formulation, the final finite element matrix problem has the following form: where v d is the discretized velocity field v, M η is the matrix resulting from the integration of the shape functions of v multiplied by the Laplacian operator, while M λ is the matrix resulting from integrating the shape functions on the operator resulting from the penalty method and multiplied by 1/λ. Since λ should be small enough, the penalty matrix M λ plays a dominant role in the problem solution, and if M λ is regular, v d = 0 becomes a trivial solution. It follows that M λ should be singular to obtain a meaningful solution [21]. However, integrating the penalty operator classically with FEM conforming elements always gives regular matrices. To render M λ singular, one should lower the quadrature order while integrating the penalty operator [33]. The chosen value of λ for the simulations in this work is set to 10 −5 .

Solving the nonlinear problem for n < 1
To solve the nonlinear problem, we start by solving the problem for n = 0.9, considering the approximation of D eq from the velocity fields computed for n = 1. Therefore, the problem depicted in (10) becomes linear, using: and therefore η ≈ K 2 (D : D) @n=1 0.9−1 . Eventually η is now a function of K , ω and the position vector x = (x, y, z). To solve for fluid velocity v, η should be written in a separated form. Therefore, we construct the full η tensor and then write it in the form: the index 1 on η designating the iteration number. One approach would be using a high-order singular value decomposition known as HOSVD [34,35], η being a 3D tensor η = f (x, ω, K ). The other approach, adopted by the authors, is constructing η K (x, ω) for each value of K in the mesh and separating η K (x, ω) into a sum of products of functions of x and ω using a simple Singular Value Decomposition (SVD) such as: Constructing the separated form of η @n=0.9 1 from the η @n=0.9 for all values of K can be easily done by multiplying each η @n=0.9 K 1 by an indicator function and summing all the obtained functions using: where p is the number of nodes in dimension K , f i (K ) is an indicator function equal to 1 on node i and 0 elsewhere. Such an approach has the advantage of high fidelity reconstruction of the separated form of the viscosity when compared to the HOSVD, which would definitely generate information losses.
Replacing the separated form of η @n=0.9 1 in the weak form of the problem (10) gives a linear problem that we can solve using the same technique explained in Section 3.1, since the resulting problem is linear. Thus, we obtain v @n=0.9 1 , a first iteration in the resolution of the velocity field at n = 0.9.
Once the velocity field v @n=0.9 1 is obtained using η @n=0.9 1 , we compute η @n=0.9 2 using v @n=0.9 1 . Later on, using η @n=0.9 2 we solve for v @n=0.9 2 and then compute η @n=0.9 3 , etc. The iterative process continues until reaching convergence on η @n=0.9 i , for instance using the L ∞ norm: At this stage, v @n=0.9 i is considered the solution for v at n = 0.9, the second snapshot of the problem: the first snapshot corresponding to the Newtonian solution with n = 1. The process is repeated for n = 0.8 starting with an approximation of D eq using v @n=0.9 i , and so on, for all the required snapshots.

Results and discussion
The non-Newtonian simulation for a single value of n is computed within an hour, whereas the full 6D manifold is computed within a few hours on a normal portable PC, core i7 with 8 GB of RAM. The results are two generic abaci (or two modern charts) containing the velocity fields v of all the possible combination of the chosen parameters K , ω and n, for both the slotted and solid vane tests. Figure 5 illustrates the velocity fields on an horizontal section in the solid vane test for n = 1 and 0.5 at K = 5000 Pa·s and ω = 0.0514 rad/s. Figure 6 illustrates the velocity fields on a horizontal section in the slotted vane test for the same values of the parameters. Figures 7(a) and (b) illustrate the velocity magnitude on a plot along the y-axis of the slices shown in Figures 5 and 6 respectively. We can clearly see that the velocity profiles drop drastically with decreasing values of n and thus increase the gradients, for the same value of K = 5000 Pa·s.
The domain mesh contains 14,096 nodes in the 3D physical domain for the solid vane mesh, and 14,563 nodes for the slotted vanes domain. One PGD converged solution takes 85 seconds on average to converge in v(x, y, z, K , ω) for a given value of n. The PGD convergence is achieved when the residual of the partial differential equation R PDE is reduced by a factor larger than PGD = 10 6 [32]:    The number of iterations of the nonlinear algorithm, after initialization with the solution for higher values of n, gradually increases from two iterations at n = 0.9 up to seven iterations of the fixed point algorithm at n = 0.1. The nonlinear loop is assumed to converge once the apparent viscosity changes by the infinity norm become negligible, as illustrated in (25).
Once the velocity fields are computed, we calculate the torque. For the sake of simplicity, we computed the torque on the outer solid cylinder. The torque T is therefore determined as: where A † corresponds to the area of the outer cylinder and the disk support on the bottom of the fluid domain. Moreover, r is the radius measured from the central axis of rotation to the differential area dA. Since v is computed as a function of ω, K and n, the result is therefore a torque manifold T (ω, K , n). Leveraging the separated variables representation of v, the torque T can be written as: where f n ni is an indicator function equal to 1 on node ni and 0 elsewhere (ni goes to 10, the number of simulations performed for different n values). Figure 8 illustrates the variation of the torque T as a function of the rotational speed ω for different values of n, choosing K = 1072 Pa·s in the solid vane test, whereas Figure 9 illustrates the same results for the torque in the slotted vane test. The result shows a linear variation for the Newtonian case n = 1, whereas the variation becomes highly nonlinear with decreasing values of n.
Moreover, Figures 8 and 9 show a difference in the torque experienced on the outer support for the same combination of parameters. The same torque shapes can be seen in both figures, while the values of the torques are lower in the slotted vane case. Such behavior is illustrated experimentally in previous publications, e.g. [16]. However the authors' interpretation-modeling the cement paste as a yielding solid-does not convincingly explain the lower torque experienced   Using our torque manifold, we can also illustrate the torque for a given value of the rotational speed ω. Figure 10 shows the torque as a function of K and n for an imposed vane angular velocity ω = 0.01 rad/s. Figure 11 illustrates the torques in the slotted vane for the same combination of the parameters.
One may notice the high dependency of the torque on the chosen values of the parameters K , n and ω. It can also be seen, by comparing Figures 10 and 11, that lower torque values are computed in the slotted vanes.

Properties identification
Once the torques are available from the post-processing of the 6D abaci as explained in Section 4, we move on to properties identification from experimental measurements. One may note that all the sensitivities with respect to the parameters are readily available. For example, if one wishes to compute the sensitivity of the torque with respect to the imposed angular velocity for example, one needs to differentiate the functions of ω: Experimental measurements are found in the literature for the chosen dimensions of the simulated vane tests in [16] at ω exp = 0.01 rpm. If we designate the experimental measurements by M solid as the torque experienced at steady state in the solid vane test and M slotted as the torque experienced at steady state in the slotted vane test, one may write a system of two equations: The above system can be solved using Newton's algorithm. Designating by X = (K , n) and the vector of cost functions F as: Figure 12. The pressure field experienced by the vanes for the identified non-Newtonian power-law parameters.
one iteration of the Newton algorithm can be written as: After convergence, the Newton algorithm gives the values of n = 0.1883 and K = 1465.4 Pa·s for the cement paste tested in [16].

Validation
To further validate the parametric simulation performed using the PGD, computational fluid dynamics (CFD) simulations were performed using the commercial software ANSYS Fluent 19. The CFD solver uses the cell-centered finite-volume method to discretize the continuity and linear momentum equation. The paste is modeled as a non-Newtonian power-law fluid after adopting the identified values of K and n using the PGD method, as described in the previous section. Following a mesh sensitivity analysis, the final mesh consists of around one million tetrahedral cells refined near the vanes. A single rotating reference frame is used with a rotational speed of 0.01 rpm. The SIMPLEC scheme is adopted for velocity-pressure coupling with a least-square cell-based scheme for the gradient spatial optimization. The pressure discretization is performed adopting the Presto scheme and a third-order MUSCL method is used for the momentum equations. A wrapped-face gradient corrector is used to enhance the convergence. Figures 12(a) and (b) illustrate the pressure fields while Figures 13(a) and (b) illustrate the shear stress field experienced by the vanes for an angular velocity ω = 0.01 rpm. As expected, the highest shear stress appears on the tip of the blades, creating therefore the highest torque. We also note that the torque due to shear stresses is dominant as compared to the torque due to normal pressures.
The fluid velocity fields are illustrated in Figure 14 for both the solid and slotted vanes. We can clearly identify a sharp decrease in the fluid velocity fields magnitude near the blade. Such phenomenon is explained by the low value of the fluid behavior index n approaching 0.18 for the fresh cement paste studied in [16].
The total torque identified by the CFD simulation on the commercial software were 1.91 mN·m and 1.73 mN·m for the solid and slotted vanes respectively. The identified torques are within 16 and 22% error respectively from the experimental values identified in [16].

Conclusion
In this work, we illustrate a novel approach for using the PGD in the simulation of multidimensional nonlinear problems. The approach relies on the use of the separated representation to illustrate the nonlinear terms of the equation before starting the PGD algorithm. The application shown corresponds to the simulation of vane tests conducted on a power-law fluid with two different geometrical configurations. The results illustrate the steady-state velocity profiles in 3D for different combinations of the chosen parameters: the imposed angular velocity ω, and the two non-Newtonian power-law parameters K and n.
A post-processing of the 6D abaci results in 3D torque manifolds: T (ω, K , n). The torque results replicate experimental observations such as a smaller torque incurred by slotted vanes as compared to solid vanes, for the identified parameters. Using experimental results from the literature, we were able to identify the non-Newtonian fluid properties through the simple resolution of a nonlinear system of two equations with two unknowns. The experimental data used for this work correspond to cement pastes, however, the same approach can be applied to identify the properties of any non-Newtonian power-law fluid.