## 3 Numerical Algorithms and Computational Methodologies
### LCR Technique  
The Log-Conformation-Reformulation (LCR) approach is implemented for numerical stabilisation in the _RheoMultiRegionFoam_ solver. In the LCR approach, the polymeric-stress tensor, \(\mathbf{\tau}\), is replaced with the logarithm of the conformation tensor, \(\mathbf{D}\), where:  
\[\mathbf{\tau}=\frac{\eta_{p}}{\lambda}(\mathbf{D}-\mathbf{I}). \tag{9}\]  
The process to recast the viscoelastic-stress constitutive-equations in terms of the logarithm of conformation tensor (log-conformation) requires that the stress tensor, \(\mathbf{\tau}\), be initially replaced (in the viscoelastic-stress constitutive-equations) by the conformation tensor, \(\mathbf{D}\), via Equation (9). This initial process leads to the alternative constitutive equations for the conformation tensor:  
\[\frac{\partial\mathbf{D}}{\partial t}+(\mathbf{U}\cdot\
abla)\mathbf{D}-(\mathbf{D}\cdot\
abla \mathbf{U}^{T}+\
abla\mathbf{U}\cdot\mathbf{D})=-\frac{1}{\lambda}f_{R}(\mathbf{D}), \tag{10}\]  
Figure 2: Computational mesh for the double-pipe geometry.  
where \(f_{R}(\mathbf{D})\) is a polynomial function in the conformation tensor [21]. Since \(\mathbf{D}\) is positive definite, it can be diagonalised that:  
\[\mathbf{D}=\mathbf{R}\cdot\mathbf{\Lambda}\cdot\mathbf{R}^{T}, \tag{11}\]  
where \(\mathbf{R}\) is an orthogonal matrix containing, as its columns, the eigenvectors of \(\mathbf{D}\), and \(\mathbf{\Lambda}\) is a diagonal matrix with corresponding eigenvalues of \(\mathbf{D}\). The second step in the LCR technique is based on the realisation that it would be more efficient to solve constitutive equation in terms of the logarithm of \(\mathbf{D}\) (i.e., introduce and solve for, \(\mathbf{\Theta}=\log\mathbf{D}\), in the constitutive equations) rather than solving directly for \(\mathbf{D}\). In particular:  
\[\mathbf{\Theta}=\log\mathbf{D}=\mathbf{R}\cdot\log\mathbf{\Lambda}\cdot\mathbf{R}^{T}. \tag{12}\]  
Following [52; 53], the velocity gradient, \(\
abla\mathbf{U}\) can be decomposed as:  
\[\
abla\mathbf{U}=\mathbf{\Omega}+\mathbf{B}+\mathbf{N}\cdot\mathbf{D}^{-1}, \tag{13}\]  
where both \(\mathbf{\Omega}\) and \(\mathbf{N}\) are anti-symmetric and \(\mathbf{B}\) is a symmetric traceless tensor which commutes with \(\mathbf{D}\). Equation (10), therefore, reduces to:  
\[\frac{\partial\mathbf{\Theta}}{\partial t}+(\mathbf{U}\cdot\
abla)\mathbf{\Theta}-(\mathbf{ \Omega}\mathbf{\Theta}-\mathbf{\Theta}\mathbf{\Omega})-2\mathbf{B}=\mathbf{R}\bigg{[}\frac{1}{ \lambda}(\mathbf{\Lambda}^{-1}-1)\bigg{]}\mathbf{R}^{T}. \tag{14}\]  
In a two-dimensional flow, the eigen-decomposition of \(\mathbf{D}\) can be expressed as:  
\[\mathbf{D}=\mathbf{R}\begin{bmatrix}\lambda_{1}&0\\ 0&\lambda_{2}\end{bmatrix}\mathbf{R}^{T}, \tag{15}\]  
where \(\lambda_{1}\) and \(\lambda_{2}\) are eigenvalues of \(\mathbf{D}\), and also:  
\[\mathbf{R}^{T}(\
abla\mathbf{U})\mathbf{R}=\begin{bmatrix}m_{11}&0\\ 0&m_{22}\end{bmatrix},\quad\mathbf{B}=\mathbf{R}\begin{bmatrix}m_{11}&0\\ 0&m_{22}\end{bmatrix}\mathbf{R}^{T},\quad\Omega=\mathbf{R}\begin{bmatrix}\delta&0\\ 0&-\delta\end{bmatrix}\mathbf{R}^{T}, \tag{16}\]  
with:  
\[\delta=\frac{(m_{12}+m_{21})}{(\lambda_{2}-\lambda_{1})}.\]