In order to analyze the viscous coupling, the first benchmark used was the layered Poiseuille flow, focusing on the kinematic viscosity ratio. In the present case, the setup consists of a two-dimensional flow between infinite parallel plates where a central layer of fluid $N$ with thickness $2a$ is sandwiched between layers of fluid $w$ (wetting), which is in contact with the solid phase. The total thickness of the channel is $2b$, as illustrated in Fig. below.
The fluids were considered identical except for their kinematic viscosity that were set as $\nu_{n}$ and $\nu_{w}$, for $n$ (non-wetting) and $w$, respectively. A constant force $G$ is applied throughout the channel, pointing in the z-direction. Simulations were performed across viscosity ratios $M=\frac{\nu_{n}}{\nu_{w}}$ ranging from $10^{-6}$ to $10^{6}$, using three different mesh refinements. The results were compared with the analytical solution given by:
$$
u_{z}^{ana}\left(y\right)=\begin{cases}
A_{n}y^{2}+C_{n} & \text{for } 0\leq|y|\leq a, \\
A_{w}y^{2}+B_{w}y+C_{w} & \text{for } a\leq|y|\leq b,
\end{cases}
$$
where the coefficients are:
$$
A_{n} = - \frac{G}{2\rho_n\nu_n}, \quad A_{w} = - \frac{G}{2\rho_w\nu_{w}}, \quad B_{w} = 2\left(A_{n}M\gamma-A_{w}\right)a, \quad C_{n} = \left(A_{w} - A_{n}\right)a^{2} - B_{w}\left(b-a\right) -A_{w}b^{2}, \quad C_{w} = -A_{w}b^{2} - B_{w}b,
$$
having $\gamma =\frac{\rho_n}{\rho_w}$ and $\rho_w=\rho_n =1$.
The results obtained using LBPM are shown in the Fig. below. Notice a significant error in the fitting between the numerical results and the analytical solution, even for small viscosity ratios.
This error can be attributed to the relaxation time $\tau$ interpolation in the interface region, once $\nu \propto (\tau - 0.5)$. Currently, in the Color model, LBPM sets this interpolation as:
$$
\tau = \tau_{n} + \frac{(1-\phi)}{2}(\tau_{w} - \tau_{n}).
$$
This leads to a finite interface with varying viscosity, which is not ideal for describing problems where the scale of the interface transition is not relevant. A possible solution to this problem is using the following interpolation:
$$
\tau = \tau_{n} + \frac{\tanh{\left(\phi\right)}}{2}(\tau_{w} - \tau_{n}),
$$
which create an approximatelly step function transition of viscosity across the fluid interface. As depicted in Figure below, this will result in a better agreement with the analytical solution. Notice an excellent agreement for viscosity ratios up to $10^{6}$ and extendable to higher values in this problem, but the number of time steps to reach the steady-state regime is very high and may be infeasible depending on your computational capacity.
Figure below shows the error $e$ for both interpolations, where $e$ is obtained through:
$$
e = \sqrt{\frac{\sum_{y=-b}^{b}\left[ u_{z}(y)-u_{z}^{ana}(y)\right]^{2}}{\sum_{y=-b}^{b}\left[ u_z^{ana}(y)\right]^{2}}}
$$
Consequently, we suggest changing the correlation in order to improve the viscous coupling description.
