# Different Formulations to Solve the Giesekus Model for Flow between Two Parallel Plates
## 3 Different Formulations to Obtain the Fully Developed Laminar Flow with the Giesekus Model
### High-Order Simulation (HOS)  
For the numerical simulation of the Giesekus fluid flow, a high-order method is adopted. In order to eliminate the pressure term in the Navier-Stokes Equation (2), the vorticity-velocity formulation is adopted. Thus, the vorticity component in the \(z\) direction, \(\omega_{z}\), can be written as:  
\[\omega_{z}=\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}. \tag{33}\]  
The equation system to be solved is given by:  
\[\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=0, \tag{34}\]  
\[\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}=- \frac{\partial\omega_{z}}{\partial x}, \tag{35}\]  
\[\frac{\partial\omega_{z}}{\partial t}+\frac{\partial(u\omega_{z})}{\partial x }+\frac{\partial(v\omega_{z})}{\partial y}=\frac{\beta}{Re}\left(\frac{ \partial^{2}\omega_{z}}{\partial x^{2}}+\frac{\partial^{2}\omega_{z}}{ \partial y^{2}}\right)+\frac{\partial^{2}T_{xx}}{\partial x\partial y}+\frac {\partial^{2}T_{xy}}{\partial y^{2}}-\frac{\partial^{2}T_{xy}}{\partial x^{2 }}-\frac{\partial^{2}T_{yy}}{\partial x\partial y}, \tag{36}\]  
and the Giesekus constitutive equation:  
\[\mathbf{T}+Wi_{\mathbf{T}}^{\
abla}+\frac{\alpha_{G}WiRe}{(1-\beta)}( \mathbf{T}\cdot\mathbf{T})=\frac{2(1-\beta)}{Re}\mathbf{D}. \tag{37}\]  
In this equation, the log-conformation method [17; 18] is adopted to overcome the high Weissenberg number problem--HWNP. Using this technique, a conformation tensor \(\mathbf{A}\) is adopted. The relation between \(\mathbf{T}\) and the conformation tensor \(\mathbf{A}\) is given by  
\[\mathbf{T}=\frac{(1-\beta)}{ReWi}(\mathbf{A}-\mathbf{I}), \tag{38}\]  
and  
\[\mathbf{\Psi}=\log_{a}(\mathbf{A}). \tag{39}\]  
The equation to be solved using this technique is given as follows:  
\[\frac{\partial\mathbf{\Psi}}{\partial t}+\
abla\cdot(\mathbf{u}\mathbf{\Psi} )=(\mathbf{\Omega}\mathbf{\Psi}-\mathbf{\Psi}\mathbf{\Omega})+2\mathbf{B}+ \frac{1}{Wi}a^{-\mathbf{\Psi}}(\mathbf{I}-a^{\mathbf{\Psi}})[\mathbf{I}+ \alpha_{G}(a^{\mathbf{\Psi}}-\mathbf{I})], \tag{40}\]  
where \(a=e\) is the Euler's number. All of the spatial derivatives are approximated by fifth- and sixth-order compact finite differences [19]. Time derivatives are discretized using a classical fourth-order Runge-Kutta scheme [20]. The Poisson equation is solved using a multigrid full approximation scheme (FAS) [21]. The calculation of the vorticity on the wall is performed according to [22] using compact high-order finite difference approximations. The boundary conditions adopted are Newtonian Poiseuille profile at the channel entrance (left boundary), non-permeability and no-slip conditions at the walls (upper and lower boundaries), and Neumann boundary conditions for the velocity at the channel exit (right boundary). The problem is solved as an unsteady problem and a very long channel is adopted to avoid the influence of the left boundary. The simulation is carried out until the maximum difference between the vorticity in two consecutive time steps is lower than \(10^{-9}\). The calculation sequence to obtain the main flow components by solving Equations (34)-(36) and (40) is given by:  
1. Apply a time integration for the vorticity and the extra-stress tensors \(\mathbf{\Psi}\) (Runge-Kutta method);
2. Calculate the extra-stress tensor components through the log-conformation method;
3. Calculate the right-hand side of the Poisson equation given by Equation (35);Calculate the velocity \(v\) by solving the Poisson equation--Equation (35);
* Calculate the velocity \(u\) using the continuity equation--Equation (34);
* Update the vorticity \(\omega_{z}\) at the walls;
* Apply a filter after the last step of the time integrator. The filtering strategy adopted in the last step is a 6 th order compact filter given by [23]. The filter is applied at the end of each time step integration. It consists of recalculating vorticity distribution through a tridiagonal system to eliminate the spurious oscillations that can appear in the numerical solution.