# Numerical Study of Viscoplastic Flows Using a Multigrid Initialization Algorithm
## 1 Introduction  
Computational fluid dynamics has firmly established itself across a wide range of applications and research areas. Subsequently, it has spread throughout engineering industry. This has gone hand in hand with major improvements in high-performance computing hardware and treatment of more sophisticated mathematical models. Notably, numerical simulations, along with mathematical analysis of viscoplastic flows, has assisted in better understanding the behavior and fundamental properties of fluids exhibiting a yield stress in different applications relevant to natural and engineering sciences [1, 2, 3, 4]. Viscoplastic flows account for an important branch of non-Newtonian fluid mechanics, as a considerable proportion of fluid flows occurring both naturally and industrially is known to exhibit yield stress. Suspensions of particles in a fluid matrix, such as slurries, gels, and nanocomposites, typically display viscoplasticity, and the study of their behavior has led to a large literature base. These fluids behave as solids as the shear rate tends to zero, and as liquids beyond a certain critical shear stress level, i.e., the yield stress. The fluid flow field is therefore divided into unyielded and yielded fluid regions. A comprehensive review of viscoplastic fluids has been carried out by [PERSON] [5]. The most widely known model for inelastic viscoplastic fluids was first proposed by [PERSON] [6], which describes a linear interrelation between shear stress and shear rate during flow. The constant ofproportionality of this linear relationship is often called the plastic viscosity of the yielded fluid, \(\mu\). Thus, the [PERSON] model [6] can be expressed, in tensorial form, as follows:  
\[\left\{ \begin{array}{ll}\dot{\gamma}=0&\text{if }|\tau|\leq\tau_{y}\\ \tau=\left(\mu+\frac{\tau_{y}}{|\dot{\gamma}|}\right)\dot{\gamma}&\text{if }|\tau|>\tau_{y} \end{array} \right. \tag{1}\]  
where \(\tau\) is the stress tensor and \(\dot{\gamma}\) is the shear rate tensor, \(\dot{\gamma}=\
abla u+(\
abla u)^{T}\), \(\
abla u\) being the velocity gradient tensor. \(|\tau|\) and \(|\dot{\gamma}|\) are evaluated based on the second invariants of the stress and the rate of strain tensors respectively, as:  
\[|\tau|=\left[\frac{1}{2}\tau:\tau\right]^{1/2} \tag{2}\]  
\[|\dot{\gamma}|=\left[\frac{1}{2}\dot{\gamma}:\dot{\gamma}\right]^{1/2} \tag{3}\]  
Later, more general forms of the Bingham model were adopted in an attempt to capture the non-linear flow behavior beyond the critical shear stress. These slightly more complex empirical models include, among others, Herschel and [PERSON], Casson, and Robertson-Stiff models. The resulting constitutive equations combine the behavior of solids in the so-called 'rigid' regions and of non-Newtonian liquids in the 'flow' regions through discontinuous equations. Thus, simulating viscoplastic flows is likely to encounter many difficulties due to the non-differentiability of the constitutive equations and the indeterminacy of the stress tensor below the yield stress. In most cases, it is necessary to determine the location and shape of the interface between the yielded and unyielded regions at which the flow must switch from one branch of the constitutive law to the other. In Equation (1), the yield surface is represented by the locus of points where \(|\tau|=|\tau_{y}|\). Since performing numerical simulations of a yield stress material flow is not a straightforward task, two main families of solution approaches were suggested over the past decades, the regularization method and augmented Lagrangian algorithm. The former approach includes methods which approximate Equation (1) by one regularized and smooth constitutive equation, which is well determined regardless of the shear rate magnitude. The regularized equation treats the whole material domain as a fluid of variable viscosity and locally assigns a large but finite value of viscosity to the unyielded regions.