# Numerical Study of Viscoplastic Flows Using a Multigrid Initialization Algorithm
## 2 Numerical Methodology
### Discretization of the Equations  
The domain is divided into a number of control volumes (CVs) using a Cartesian grid of equally spaced horizontal and vertical grid lines as shown in Figure 2. We denote the unit vectors in the \(x\) and \(y\) directions by \(i\) and \(j\), respectively. In finite volume method, a discrete approximation of the continuity and momentum equations is obtained by integrating each equation over every cell. The resulting algebraic expressions involve the values of the unknowns \(\mathbf{u}\) and \(p\) at the centroid of each cell and at the centers of neighboring CVs. The surface flux integrals are evaluated separately on each face using various interpolation schemes. Figure 2 shows a control volume P and its neighbors, S, E, N and W. The letters P, S, E, N and W also denote the position vectors of the centers of the respective CVs. Figure 1: Papanastasiou regularization for the Bingham model: stress magnitude (**left**) and apparent viscosity (**right**) as functions of the magnitude of the shear rate tensor \(\dot{\gamma}\). Integrating Equation (12) over the volume \(\Delta V\) of the computational cell, in the absence of any source term, yields the following form:  
\[\oint_{\Delta V}\frac{\partial}{\partial x}(\rho u\phi)dV+\oint_{\Delta V}\frac{ \partial}{\partial y}(\rho v\phi)dV=\oint_{\Delta V}\frac{\partial}{\partial x }\Big{(}\Gamma\frac{\partial\phi}{\partial x}\Big{)}dV+\oint_{\Delta V}\frac{ \partial}{\partial y}\Big{(}\Gamma\frac{\partial\phi}{\partial y}\Big{)}dV \tag{13}\]  
By integrating the governing equation over a control volume and applying Green-Gauss theorem, the semi-discretized form of the equation, at its nodal point P, can be written as follows :  
\[(\rho uA\phi)_{e}-(\rho uA\phi)_{w}+(\rho vA\phi)_{n}-(\rho vA\phi)_{s}=\left( \Gamma A\frac{d\phi}{dx}\right)_{e}-\left(\Gamma A\frac{d\phi}{dx}\right)_{w} +\left(\Gamma A\frac{d\phi}{dy}\right)_{n}-\left(\Gamma A\frac{d\phi}{dy} \right)_{s} \tag{14}\]  
where A is the cross-sectional area of the control volume face. With the aim of obtaining an useful form of the discretized equation, the diffusion coefficients \(\Gamma\), the mass fluxes \(\rho u\) and \(\rho v\) and the gradients \(d\phi/dx\) and \(d\phi/dy\) at the cell faces \(e\), \(w\), \(n\) and \(s\) are required. At the control volume faces an approximate distribution of properties between nodal points is used. Following well-established practice, linear approximations seem to be the obvious and simplest way of calculating interface values and the gradients. This practice is called central differencing. It is convenient to define two variables \(F\) and \(D\) to represent the convective mass flux per unit area and diffusion conductance at cell faces. For x direction it is written as:  
\[F =\rho u \tag{15a}\] \[D =\frac{\Gamma}{\delta x} \tag{15b}\]  
By considering Equation (15), the integrated general conservation Equation (13) can now be written as:  
\[\Delta y(F_{e}\phi_{e}-F_{w}\phi_{w})+\Delta x(F_{n}\phi_{n}-F_{ s}\phi_{s})=\\ \Delta y(D_{e}(\phi_{E}-\phi_{P})-D_{w}(\phi_{P}-\phi_{W}))+ \Delta x(D_{n}(\phi_{N}-\phi_{P})-D_{s}(\phi_{P}-\phi_{S})) \tag{16}\]  
We assume that the velocity field is known, which takes care of the mass fluxes values \(F\). To solve Equation (16), we need to calculate the transported property \(\phi\) at each cell face using appropriate numerical schemes. In the discretization of divergence  
Figure 2: Two-dimensional orthogonal grid. N, S, E and W correspond to neighbor cells of cell P; n, s, e and w denote cell P faces; \(\Delta x\) and \(\Delta y\) are cell P dimensions in the x and y spatial coordinates; \(\delta y_{n}\), \(\delta y_{s}\), \(\delta x_{e}\) and \(\delta x_{w}\) correspond to cell-center to cell-center distances from cell P to neighbor cells. terms, many schemes are available in literature (first-order upwind, second-order upwind, central differencing scheme, QUICK, etc). Higher-order schemes have been widely used in applications involving orthogonal and uniform meshes. In this study, [PERSON] et al.