Monolithic Finite Element Methods for the simulation of thixo-viscoplastic flows

This talk is concerned with the application of Finite Element Method (FEM) and Newton-Multigrid solver to simulate thixotropic flows.The thixotropy phenomena are introduced to yield stress material by taking into consideration the internal material micro structure using a structure parameter. Firstly, the viscoplastic stress is modified to include the thixotropic stress dependent on the structure parameter. Secondly, an evolution equation for the structure parameter is introduced to induce the time-dependent processof competition between the destruction (breakdown) and the construction (buildup) inhabited in the material. Substantially, this is done simply by introducing a structure-parameter-dependentviscosity into the rheological model for yield stress material, as for instance the Houska model based on a viscosity approach for the Bingham model [2].The modified viscoplastic stress w.r.t. the structure parameter which is integrated, in quasi-Newtonian manner or lagrangian multiplier manner, into the generalized Navier-Stokes equations and the evolution equation for the structure parameter constitutes the main core of full set of mod-eling equations, which are creditable as the privilege answer to incorporate thixotropy phenomena.The nonlinearity of the problem, related to the dependency of the diffusive stress on the material parameters, is treated with generalized Newton’s method w.r.t. the Jacobian’s singularities having a global convergence property. The linearized systems inside the outer Newton loops are solvedusing the geometrical multigrid methods with a Vanka-like smoother taking into account a stable FEM approximation pair for velocity and pressure with discontinuous pressure and biquadratic velocity spaces.We analyze the accuracy, robustness and efficiency of the Newton-Multigrid FEM solver [1] throughout the solution of thixotropic flow problems of benchmarking character in channel and Couette device [3].

In quasi-Newtonian modeling approach for thixo-viscoplastic flows, an extended viscosity µ(·, ·) is used for the generalized Navier-Stokes equations [10]. As for instance [13]: where k is the regularization parameter. The generalized Navier-Stokes equations and the evolution equation for the structure parameter constitute the full set of modeling equations which is given as follows: where u denotes velocity, p the pressure, λ the structure parameter, F and G the nonlinear functions for buildup and breakdown of material micro-structure. D II = 1 2 D(u) : D(u) is the second invariant of the strain rate tensor D(u).

FINITE ELEMENT DISCRETIZATION
To derive the variational form for thixo-viscoplastic flows, we consider the spaces T := L 2 (Ω), V := (H 1 0 (Ω)) 2 , and Q := L 2 0 (Ω) associated, respectively, with the corresponding L 2 -norm, ||·|| 0 , H 1 -norm, ||·|| 1 , and L 2 -norm, ||·|| 0 . Letũ := (λ, u, p) ∈ T ∩ H 1 (Ω) × V × Q, andṽ := (ξ, v, q) ∈ T × V × Q be a test function. The weak formulation for the thixo-viscoplastic flows reads: where a λ (ũ)(·, ·), a u (ũ)(·, ·), and b(·, ·) are given as follows The finite element approximations of the problem (5) have to take care of its saddle point character, due to the bilinear form (8). Furthermore, since thixo-viscoplastic flows are usually slow, the only remaining issue is the control/continuity of the bilinear form (6) in the norm of space T. We opt for higher order stable pair bi-quadratic for velocity and piece-wise linear discontinuous for the pressure, Q 2 /P disc 1 , and higher order quadratic for structure parameter Q 2 with the appropriate stabilization terms [10,15]. Indeed, let the domain Ω be partitioned by a grid K ∈ T h which are assumed to be open quadrilaterals such that Ω = int k∈T h K . For an element K ∈ T h , we denote by E(K) the set of all 1-dimensional edges of K. Let E i := k∈T h E(K) be the set of all interior element edges of the grid T h .
We define the conforming finite element spaces T h ⊂ T, V h ⊂ V, and Q ⊂ Q h such that: Q h = q h ∈ Q, q h|K ∈ P disc Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain The approximate problem reads: The stabilization term jũ(·, ·) is given as follows The stabilization (13) is consistent, control the convective terms and makes the coercivity and continuity match in T h associated with the norm |||·|||, where

GENERALIZED DISCRETE NEWTON
We use the Newton method to approximate the nonlinear residuals.
denote the residuals for the system (12). The nonlinear iteration is updated with the correction δũ,ũ k+1 =ũ k + δũ. Then, the Newton linearization gives the following approximation for the residuals: The Newton's method iterations, assuming invertible Jacobians, are given as follows: The damping parameter ω l ∈ (0, 1) is chosen such that The damping parameter is not sufficient for the convergence of this type of highly nonlinear problem, mainly due to the presence of Jacobian's singularities related to the problem or simply by being out of the domain of Newton's convergence [8,10]. We use a generalized Newton's method which consists of using approximate Jacobians far away from the quadratic convergence range or close to singularities and accurate Jacobians in the quadratic region of convergence in an adaptive way [8].
Indeed, based on a priori analysis of Jacobians property. Let the Jacobian be written as follows: The Jacobian (20) is splitted into a direct sum of corresponding operators with different properties. The parameter δ l ∈ (0, 1) is solely dependent on the rate of actual residual convergence [8]. It is worth mentioning that the operator-related damped Jacobian method (20) is related to the continuous Newton's method. The Jacobian approximation is only dependent on the rate of the actual residual convergence R l / R l−1 . This generalized Newton's method assures a global nonlinear convergence [8].
Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain

MONOLITHIC MULTIGRID LINEAR SOLVER
To develop an appropriate linear solver, we segregate the Jacobian as follows which is a saddle point problem. Then, the resulting linear system is treated with a Multilevel Pressure Schur Complement (MPSC) approach with Vanka-like smoother i.e.
In (22), we solve exactly on real element, K, and perform an outer Gau-Seidel iteration [5]. We use standard geometric multigrid solver for linearized system with standard Q 2 and P disc 1 restriction and prolongation operators. The combination of a stable finite element approximations, Q 2 /P disc 1 , for Stokes problem together with multigrid results in a high numerically accurate, flexible, and efficient FEM-multigrid solver.

THIXO-VISCOPLASTIC FLOW IN LID-DRIVEN CAVITY
Lid-driven cavity flows represent an academic common standard benchmark for incompressible CFD codes. Therefore, we present the corresponding results for Newtonian, viscoplastic, and thixoviscoplastic flows. Furthermore, this problem is accepted as a test configuration to check points wise mesh convergences despite the lack of regularity due to the pressure singularity in the corners of upper-lid. From thixotropic collection models given in Table 1, we use Houska's material model (m = 1).

Newtonian flow in lid-driven cavity
The global accuracy of the approximation which consist of the L 2 -norm of the velocity is investigated using the kinetic energy. On the other hand, the point wise accuracy is investigated using the velocity magnitude at vertical center-line beside the primary vortex and the lower left secondary vortex. In order to check the solver convergence, we list in Table 2 the kinetic energy, 1 2 Ω ||u|| 2 dx, and Newtonmultigrid iterations, the number of nonlinear iteration versus the average number of multigrid sweeps (N/M), w.r.t. mesh refinement for an increased Reynolds numbers Re = 1000, Re = 5 000, and Re = 10 000. The starting solution for any level is the interpolated one from one level coarser. Table 3 shows the primary vortex and lower left secondary vortex w.r.t. mesh refinement for Reynolds numbers Re = 1000, Re = 5 000, and Re = 10 000. Moreover, we provide in Figure 1 the stream function contours for the mesh refinement level 9 and the velocity magnitude at vertical center-line w.r.t. mesh refinement for Reynolds numbers Re = 1000, Re = 5 000, and Re = 10 000. As can be seen in Table 2, grid independent results are achieved for the kinetic energy, as well as for Newton-multigrid solver. It is worth mentioning that for the coarser levels few extra nonlinear iterations are required in contrast to finer mesh due to the decrease of interpolation error in the starting solutions.

Viscoplastic flow in lid-driven cavity flow
We acheived point wise convergence under mesh refinement for Newtonian flow. Moreover, no further improvement by increasing the resolution beyond mesh refinement level 8. Now, we investigate the impact of the regularization parameter in quasi-Newtonian modeling approach for viscoplastic flow. Figure 2 shows the boundary limit of the numerical approximation of the rigid zone w.r.t. regularization parameter k. Clearly, the relative convergence of the boundary limit of the rigid zone w.r.t. regularization parameter k is obtained. Furthermore, there is an optimal regularization K L ≈ 2 L8 ≈ 10 3 from from which there is no further accuracy improvement in capturing the rigid zone by increasing the regularization parameter k. In Figure 3, we use the optimal choice of the regularization parameter and mesh refinement level to predict the relative position of the rigid zone to stream function contours for an increased non-thixotropic yield stress parameter τ 0 = 1, τ 0 = 2, τ 0 = 5, τ 0 = 10, τ 0 = 20, and τ 0 = 50. Furthermore, we provide the corresponding Newton-multigid data in Table 4 which depicts the number of Newton-multigrid iterations, i.e. the nonlinear number of iterations and the average number of multigrid iterations (N/M), w.r.t. different regularization parameters k and mesh refinement levels L. The solutions are calculated using the continuations process w.r.t. regularization k.
From Table 4, we conclude the Newton-multigrid solver is mesh refinement independent. Clearly, the nonlinearity of the problem is increased by increasing the non-thixotropic yield stress parameter τ 0 , But, the slightly increases in the nonlinearity w.r.t. mesh refinement is due to the continuation process w.r.t. regularization parameter k used to obtain the solutions.

Thixo-viscoplastic flow in lid-driven cavity flow
Armed with the knowledge of the point wise mesh convergence of viscoplastic driven cavity flow. Indeed, we obtained the point wise convergence of the boundary limit of the rigid zone w.r.t. the regularization parameter k. Furthermore, there is pair (K, L) ≈ (2 L8 , L8) regularization and mesh refinement level beyond which no further resolution improvement is possible. Now, we are ready the investigate thixo-viscoplastic driven cavity. Figure 4 sets out the relative position of the rigid zone to stream function contours for an increased thixotropic yield stress parameter τ 1 . The solutions are calculated with the resolution barrier pair (K, L). In Table 5, we summarize the number of Newton-multigrid iterations, nonlinear number of iterations and the average number of multigrid iterations (N/M), w.r.t. different regularization parameters k and mesh refinement levels L for thixo-viscoplastic flow for different values of thixotropic yield stress parameters τ 1 . The solutions are calculated using the continuations process w.r.t. regularization k. Table 5 shows the mesh refinement independent of Newton-multigrid solver. Indeed, the nonlinearity of the problem is increased by increasing the thixotropic yield stress parameter τ 1 , But, the slightly increases in the nonlinearity w.r.t. mesh refinement is due to the continuation process w.r.t. regularization parameter k used to obtain the solutions.

SUMMARY
We presented a Newton-multigrid FEM solver for the quasi-Newtonian modeling approach for thixotropic flows. Based on a two-fields Stokes solver, we used higher order stable Q 2 /P disc 1 FE approximations for velocity and pressure and higher order Q 2 FE approximation for the structure parameter field with appropriate stabilization term. The combination of a stable finite element approximations, Q 2 /P disc 1 , for Stokes problem together with multigrid results in high numerically accurate, flexible and efficient FEM-multigrid solver. The nonlinearity is handled with generalized Newton's method w.r.t. the Jacobian's singularities having a global convergence property. For the numerical investigations; we used lid-driven cavity benchmark to find out the optimal setting, mesh refinement, and regularization. Indeed, we achieved a point-wise mesh convergence as well as a resolution barrier, (k, L) regularization mesh refinement level, beyond which no further resolution's improvement is possible. Furthermore, the solver shows a mesh refinement independency. For viscoplastic and thixo-viscoplastic solutions, we used the discrete continuation process w.r.t. regularization which might be integrated continuously within the solver in a black box manner.