# Macro-Micro-Coupled Simulations of Dilute Viscoelastic Fluids
## 2 Governing Equations
### Numerical Considerations
#### 2.3.2 FENE Dumbbells  
Because of the nonlinearity introduced in the evolution equation, a semi-implicit, first-order algorithm is used instead of the forward Euler scheme. Briefly, a two-step operator splitting procedure is implemented, where the two steps are given by Equations (26) and (27) [29] as follows:  
\[\mathbf{\bar{Q}}\big{(}t_{j+1}\big{)}=\mathbf{Q}\big{(}t_{j}\big{)}+\bigg{[}( \mathbf{\
abla}\mathbf{u})^{\intercal}\cdot\mathbf{Q}\big{(}t_{j}\big{)}-\frac{1}{ 2 De}\bigg{(}\frac{\mathbf{Q}\big{(}t_{j}\big{)}}{1-Q\big{(}t_{j}\big{)}^{2}/b} \bigg{)}\bigg{]}\Delta t_{j}+\frac{1}{\sqrt{De}}\Delta\mathbf{W}_{j}, \tag{26}\]  
and  
\begin{table}
\begin{tabular}{l c c c c} \hline \hline
**Parameter** & **Symbol—Eqn** & **Units** & **Value** & **Ref.** \\ \hline Density & \(\rho\) & kg/m\({}^{3}\) & 1000 & \\ Solvent viscosity & \(\eta_{s}\) & Pa s & 0.001 & \\ Zero shear rate viscosity & \(\eta_{0}=\eta_{s}+\eta_{p}\) & Pa s & 0.0037 & [33] \\ Characteristic relaxation time & \(\lambda\) & s & 0.010 & [34] \\ Characteristic velocity & \(U=\sqrt{\gamma/(\rho R_{0})}\) & m/s & 0.260 & \\ Characteristic length & \(L\) & m & 0.003 & \\ Surface tension & \(\gamma\) & N/m & 0.070 & [34] \\ Initial jet radius & \(R_{0}\) & m & 0.001 & [34] \\ Reynolds number & \(Re=\rho UL/\eta_{0}\) & Dimensionless & 78 & \\ Ratio of viscosities & \(\beta=\eta_{s}/\eta_{0}\) & Dimensionless & 0.270 & \\ Deborah number & \(De=\lambda U/L\) & Dimensionless & 0.800 & \\ Ohnesorge number & \(Oh=\eta_{0}/\sqrt{\rho\gamma R_{0}}\) & Dimensionless & 0.038 & \\ FENE parameter & \(b\) & Dimensionless & \(10^{5}\) & \\ Wavenumber & \(k\) & Dimensionless & 0.8 & \\ \hline \hline \end{tabular}
\end{table}
Table 1: Simulation parameters.  
\[\left[1+\frac{1}{4 De}\frac{\Delta t_{j}}{1-\frac{Q^{2}\left(t_{j+1} \right)}{b}}\right]\mathbf{Q}\big{(}t_{j+1}\big{)} = \mathbf{Q}\big{(}t_{j}\big{)}+\frac{1}{\sqrt{De}}\Delta\mathbf{W}_ {j}+ \tag{27}\] \[\frac{1}{2}\left[(\boldsymbol{\
abla}\mathbf{u})^{\intercal} \big{(}t_{j+1}\big{)}\cdot\bar{\mathbf{Q}}\big{(}t_{j+1}\big{)}+(\boldsymbol{ \
abla}\mathbf{u})^{\intercal}\big{(}t_{j}\big{)}\cdot\mathbf{Q}\big{(}t_{j} \big{)}\right.\] \[\left.-\frac{1}{2 De}\left(\frac{\mathbf{Q}\big{(}t_{j}\big{)}}{1- \frac{Q^{2}\left(t_{j}\right)}{b}}\right)\right]\Delta t_{j}.\]  
Equation (27) can be written as  
\[Q_{j+1}^{3}-P\,Q_{j+1}^{2}-b\left(1+\frac{1}{4 De}\Delta t\right)Q_{j+1}-b\,P=0, \tag{28}\]  
where \(P\) is the magnitude of the right-hand side of Equation (27). This cubic equation has one real root in the interval \(\left(0,\sqrt{b}\right)\)[29, 35]. The approach is to find the magnitude of \(\mathbf{Q}\big{(}t_{j+1}\big{)}\) using the polynomial's solution in this interval; then, its direction is found by solving Equation (27) for each vector component.  
Note that Equation (27) requires the velocity gradient tensor, \((\boldsymbol{\
abla}\mathbf{u})^{\intercal}\big{(}t_{j+1}\big{)}\), at time \(t+1\), which is unknown for non-visometric flows. In our computations, \(\mathbf{u}\big{(}t_{j+1}\big{)}\) is found by a second-order extrapolation following [35],  
\[\mathbf{u}\big{(}t_{j+1}\big{)}=2\mathbf{u}\big{(}t_{j}\big{)}-\mathbf{u} \big{(}t_{j-1}\big{)}. \tag{29}\]