# A Backwards-Tracking Lagrangian-Eulerian Method for Viscoelastic Two-Fluid Flows
## 3 Numerical Method
The trajectory \(\mathbf{x}(t)\) for \(t\in I_{n}\) is first calculated by solving (8) backwards in time. Subsequently, when the trajectory is known, the viscoelastic stress \(\boldsymbol{\tau}(t_{n},\mathbf{x}(t_{n}))\) is interpolated from the viscoelastic stresses from the solution of the previous simulation step. The Lagrangian constitutive Equation (7) is then solved forwards in time and the viscoelastic stress is thus obtained at the final location \(\mathbf{x}(t_{n})\), i.e., at the predefined location. Generally, the choice of locations for the Lagrangian nodes is somewhat arbitrary. However, by choosing the locations of the Eulerian grid nodes, the connectivity and structure of the octree grid may be utilized when the stress at time \(t_{n}\) is interpolated to the position of a Lagrangian node \(\mathbf{x}(t_{n})\). The stress at the node at time \(t_{n}\) is obtained through bilinear or trilinear interpolation, respectively, for two-dimensional (2D) and three-dimensiona (3D). Furthermore, the viscoelastic stress contribution to the discretized momentum equation may be calculated directly from node stresses. Therefore, no unstructured interpolation of viscoelastic stresses is necessary. The calculation of this contribution is described in detail further on in this section. Figure 2 shows a schematic description of the steps that are involved in the algorithm. In summary, the performed steps are:  
1. Calculate the Lagrangian node trajectory by solving (8) backwards in time, starting at the predefined location \(\mathbf{x}(t_{n+1})\). 2. Interpolate the stress \(\boldsymbol{\tau}(t_{n},\mathbf{x}(t_{n}))\) to the Lagrangian node from the known node values at time \(t_{n}\). 3. Solve (7) forwards in time along the trajectory \(\mathbf{x}(t)\), \(t\in I_{n}\). The ordinary differential Equations (7) and (8) may be solved with any appropriate choice of ODE solution algorithm. In this work, the fourth order Runge-Kutta method, which is commonly referred to as the RK4 method [33], is used. Furthermore, \(N_{\text{sub}}\) equally sized substeps of length \(\Delta t_{\text{sub}}\) are used, being defined, such that \(N_{\text{sub}}\Delta t_{\text{sub}}=\Delta t=t_{n+1}-t_{n}\). In the current work, \(N_{\text{sub}}=3\) is used. When solving (7) and (8), local quantities that are stored on the Eulerian grid are required along the trajectory of the Lagrangian nodes. More specifically, the velocity \(\mathbf{u}\) is needed for the backwards-tracking and the velocity gradient \(\
abla\mathbf{u}\) is needed to solve the constitutive equation. The velocities are stored at the Eulerian cell centers. When the local velocity is required, the cell centers containing the Lagrangian node are identified and the velocity is interpolated to its location while using bilinear or trilinear interpolation, respectively, for 2D and 3D. When the velocity gradient is required, it is calculated from the bilinear or trilinear interpolation formula using the same interpolation basis. For the interpolation of a field \(\phi\) to a location within a box spanned by the corners \((x_{i},y_{j},z_{k})\), \(i,j,k\in\{0,1\}\), the trilinear interpolation formula in order to calculate the interpolant \(\hat{\phi}\) reads  
\[\hat{\phi}=\sum_{i=0}^{1}\sum_{j=0}^{1}\sum_{k=0}^{1}c_{x}^{i}c_{y}^{j}c_{z}^{k }\phi_{ijk}, \tag{10}\]  
where \(\phi_{ijk}\) is the value of \(\phi\) at the corner \((x_{i},y_{j},z_{k})\) and with the coefficients  
\[c_{l}^{i}=\left\{ \begin{matrix}1-\bar{x}_{l}/\Delta x_{l},&i=0\\ \bar{x}_{l}/\Delta x_{l},&i=1\end{matrix} \right.\quad,\quad l=x,y,z, \tag{11}\]  
where \(\bar{x}_{l}\) is the \(l\)th coordinate of the interpolation position relative to the lower corner of the box and \(\Delta x_{l}\) is the size of the box that coordinate direction.