# Numerical Simulation of KBKZ Integral Constitutive Equations in Hierarchical Grids
## 3 Numerical Method
### HiGTree/HiGFlow System
#### 3.1.3 Calculation of the Extra-Stress Tensor \(\mathbf{\tau}(\mathbf{x},t_{n+1})\)  
We follow the methodologies described in [PERSON] et al. (2018) and [PERSON] et al. (2019) to calculate the extra-stress tensor \(\mathbf{\tau}(\mathbf{x},t_{n+1})\). The constitutive Equation (3) can be written as follows:  
\[\mathbf{\tau}(t_{n+1})= \int_{-\infty}^{t-t_{c}}M(t_{n+1}-t^{\prime})H(I_{1},I_{2})\mathbf{ B}_{t^{\prime}}(t_{n+1})dt^{\prime} \tag{11}\] \[+\int_{t-t_{c}}^{t}M(t_{n+1}-t^{\prime})H(I_{1},I_{2})\mathbf{B}_ {t^{\prime}}(t_{n+1})dt^{\prime}\,\]  
where \(t_{c}=t\) if \(t<s_{c}\) or \(t_{c}=s_{c}\) if \(t\geq s_{c}\). The parameter \(s_{c}\) (\(s_{c}\) is a time interval) depends on the relaxation parameter \(\lambda_{ref}\). This methodology is called \(s\)-approach and is described in more detail in [PERSON] et al. (2018). Now, let \(t^{\prime}_{j}=0,1,\cdots,N\), be \((N+1)\)-points in the interval \([t_{n+1}-t_{c},\ t_{n+1}]\), where \(N\) is a fixed number. Then, the integral equation can be written as:  
\[\mathbf{\tau}(t_{n+1})= \int_{-\infty}^{t-t_{c}}M(t_{n+1}-t^{\prime})H(I_{1},I_{2}) \mathbf{B}_{t^{\prime}}(t_{n+1})dt^{\prime} \tag{12}\] \[+\sum_{j=0}^{N-2}\int_{t^{\prime}_{2j}}^{t^{\prime}_{j+2}}M(t_{n +1}-t^{\prime})H(I_{1},I_{2})\mathbf{B}_{t^{\prime}}(t_{n+1})dt^{\prime}\,\]  
where \(t^{\prime}_{0}=0\) or \(t^{\prime}_{0}=t_{n+1}-t_{c}\). We consider \(\mathbf{B}_{t^{\prime}}(t_{n+1})=\mathbf{B}_{t-t_{c}}(t_{n+1})\) for \(t^{\prime}<t_{n+1}-t_{c}\) and, therefore, the first integral becomes:  
\[\int_{-\infty}^{t-t_{c}}M(t_{n+1})H(I_{1}(\mathbf{B}_{t-t_{c}}(t_{n+1})),I_{2 }(\mathbf{B}_{t-t_{c}}(t_{n+1})))\,\mathbf{B}_{t-t_{c}}(t_{n+1})dt^{\prime} \tag{13}\]  
which can be solved without any further issues. Regarding the integrals within the summation operator in Equation (12), we use the method of undetermined coefficients (with a second-order quadratic formula) for their calculation (for details, see [PERSON] et al. (2018)). In the following sections, we describe the method used to compute the tensor \(\mathbf{B}_{t^{\prime}(t_{n+1})}(t_{n+1})\) and how the points \(t^{\prime}_{j}(t_{n+1})\) are calculated. * Discretization of the time interval \([t-t_{c},\ t_{n+1}]\) One of the key issues of the _deformation fields method_ is how the integration nodes \(t-t_{c}=t^{\prime}_{0}<t^{\prime}_{1}<\cdots<t^{\prime}_{N}=t_{n+1}\) are distributed over the interval \([t-t_{c},\ t_{n+1}]\), because such distribution can affect the accuracy of the results when solving complex flows. In [PERSON] et al. (2019), the authors presented one discretization using a function that allowed them to determine the distribution of the time-integration points, which showed excellent results in some of the specific flow cases studied (such as extensional flows). However, care must be taken if we plan to generalize these resultsto more complex flows. In this work, we decide to use the more generic methodology presented in [PERSON] et al. [21]. We consider time-dependent flows so that the integration nodes are calculated using a _geometric progression_ at time \(t_{n+1}\) as follows:  
1. Set \(t_{0}^{\prime}=t-t_{c}\) and \(t_{N}^{\prime}=t_{n+1}\);
2. Using \(t_{c}\), where \(t_{c}=t_{n+1}\) if \(t<s_{c}\) or \(t_{c}=s_{c}\) if \(t\geq s_{c}\) make \(t_{N-j}^{\prime}=t_{N}^{\prime}-\delta t\ q^{j}\), \(j=1,2,\cdots,N-1\), where \(q=(t_{c}/\delta t)^{1/N}\), \(\delta t\) is the time-step. * Computation of the Finger tensor \(\mathbf{B}_{t^{\prime}(t_{n+1})}(\mathbf{x},t_{n+1})\) One of the difficulties in the numerical simulation of viscoelastic flows using integral constitutive models is how to calculate accurately the strain history. In finite elements, this can be accomplished by a particle-tracking method based on the velocity field (see [12]), but, here, a different approach is taken.