Nonlocality of the slip length operator for scalar and momentum transport in turbulent flows over superhydrophobic surfaces

The incompressible Navier-Stokes equations are given below and consist of continuity and momentum transport.

Here, $\rho$ is density, $p$ is pressure, and $\nu$ is kinematic viscosity. The scalar transport equations are

The Reynolds-averaged Navier-Stokes equations (RANS) are obtained by averaging the Navier-Stokes equations in homogeneous directions. In our problem, these are the streamwise ($x_{1}$) and spanwise ($x_{3}$) directions, as well as time.

The rightmost term is the divergence of Reynolds stress and is unclosed. It is often modelled with an eddy viscosity $D_{v}$, such that the original equation simplifies to

In the Reynolds-averaged scalar transport equation, the divergence of turbulent flux is modelled analogously, with eddy diffusivity $D_{c}$.

We perform direct numerical simulations of turbulent channel flow with homogenized slip and pattern-resolved boundary conditions. The three-dimensional Navier-Stokes and scalar transport equations are discretized and solved by code originally developed by Bose et al. (2010) [30] and modified for superhydrophobic textured boundaries by Seo and Mani (2016) [27]. Simulations are performed at $Re_{\tau}\approx 197.5$ and $Pr=1$. The computational domain is $2\pi h\times\pi h\times 2h$ in the streamwise ($x_{1}$), spanwise ($x_{3}$), and wall-normal ($x_{2}$) directions, where $h$ is the channel half height. We use a second-order finite difference scheme on a staggered mesh, with uniform mesh spacing in the $x_{1}$ and $x_{3}$ directions, and stretched mesh in the $x_{2}$ direction. The grid resolution is $\text{max}(\Delta x_{1}^{+})=6.4$, $\text{max}(\Delta x_{3}^{+})=3.2$, $\text{min}(\Delta x_{2}^{+})=0.5$, and $\text{max}(\Delta x_{2}^{+})=6.4$. Additionally, grid resolution is chosen such that at least eight grid points resolve the texture width ($\Delta x_{1}^{+}\leq w^{+}/8$ and $\Delta x_{3}^{+}\leq w^{+}/8$). Second-order Adams-Bashforth and Crank-Nicolson schemes were used for time advancement [31].

For the gas-liquid interface, a Navier slip boundary condition was implemented for the homogenized slip cases and a patterned free-shear and no-slip condition was implemented for the pattern-resolved cases [12]. We assume infinite surface tension, zero Weber number, and flat interfaces for the pattern-resolved cases. Both the top and bottom wall are treated as superhydrophobic surfaces in all simulations. Simulations were conducted for square posts of texture wavelengths $\lambda^{+}=77.6$, $\lambda^{+}=103.4$, and $\lambda^{+}=155.1$, and corresponding homogenized slip lengths. In outer units, the listed pattern lengths correspond respectively to $\pi h/8$, $\pi h/6$, and $\pi h/4$. Streamwise-aligned ridges with texture wavelength of $\lambda^{+}=155.1$ were also simulated. All simulations were with solid fraction of $\phi_{s}=w^{2}/\lambda^{2}=1/9$. See Table 1 for a summary of cases.

To evaluate eddy diffusivity and eddy viscosity operators, we utilize the macroscopic forcing method [28]. Here, we shall present our methodology with the scalar transport equations, but this analysis is also applied to the momentum transport equations.

We apply forcing $s_{\text{MFM}}$ to the scalar field such that the mean scalar gradient takes on the Dirac delta ($\delta$) distribution.

In this form, the eddy diffusivity is a second-order tensorial kernel $D_{c,jk}$, which can account for both anisotropy and nonlocality. Using a Kramers-Moyal expansion [35], we can further decompose this eddy diffusivity into a leading-order term that characterizes local behavior and higher-order terms that characterize nonlocal behavior.

The leading-order term, $D_{c,jk}^{0}(\textbf{x})$ is the zeroth moment of the eddy diffusivity kernel, and can be taken as a local approximation of the eddy diffusivity operator [32]:

Because we are simulating a turbulent channel, we will focus on the components corresponding to the dominant scalar gradient $\frac{\partial\bar{c}}{\partial x_{2}}$. Thus,

And the leading order approximation is

Similar analysis is carried out for the momentum transport equation. The generalized eddy viscosity operator and relevant leading order approximation are

We further extend the macroscopic forcing method to evaluate a generalized slip kernel of the form

In the Stokes flow limit ($Re\ll 1$), the velocity slip length kernel simplifies to a $\delta$ function, corresponding to purely local behavior as represented by the standard homogenized slip length boundary condition [14, 15] and by the leading order moment. We will evaluate the kernel at finite $Re$ to quantify any possible nonlocality that arises.

At this point, we would like to note that previous MFM analysis has only been performed on flows where the boundary conditions are already consistent with macroscopic forcing. The question arises of how to define a macroscopic boundary condition for our forced equations that is consistent with the pattern-resolved slip boundary condition in the DNS equations. We relax the no-slip boundary condition on the solid surfaces of the pattern geometry to allow for a mean, forced value $s_{c}$ or $s_{v}$ for scalars and momentum, respectively.

To motivate our consideration of nonlocal effects on slip length, we shall first consider the problem of scalar mixing. The scalar transport equations are linear, and so present a simplified baseline from which to analyze nonlocality. For this problem, we shall choose temperature as our scalar, but the results shown here are generally applicable to any passive scalar. Past work has shown that patterned slip surfaces increase scalar transport of many interfacially driven transport phenomena [36, 37, 38].

Consider the problem of turbulent channel flow. Both top and bottom walls of the channel are superhydrophobic surfaces with square post geometry, of solid fraction $\phi_{s}=1/9$ and texture wavelength $\lambda^{+}=155.1$. The bottom wall is held at a temperature of $+T_{\text{wall}}$ and the top wall is held at a temperature of $-T_{\text{wall}}$. For simplicity, we choose $T_{\text{wall}}=1$. We will compare the scalar mixing of cases with and without volumetric heating $Q$.

We can fit a slip length for this problem formulation as below. For temperature, this slip length has also been defined as the Kapitsa length [38, 39, 40].

While in the low thermal Peclet limit, slip length is shown to solely depend on the pattern geometry [14]. In the finite thermal Peclet limit, it is expected to additionally depend on the underlying flow and molecular diffusivity [41, 40]. But it should not be affected by the presence of uniform heat generation. To test this premise, we choose $Q=0.025$, which produces a mean temperature profile in which $\frac{\partial T}{\partial x_{2}}$ varies significantly in the near-wall region of the top wall.

We can then calculate the slip length for each of these cases. The average slip length of the zero $Q$ case is $23.3$. For the nonzero $Q$ case, we will consider the slip length of the bottom and top walls separately. The bottom wall has an average slip length of $23.2$ and a negligible deviation of $0.16\%$ from the fitted value based on the zero heating case. The top wall, which is where see significant variance in temperature gradient, has average slip length $21.5$ and a substantial deviation of $8\%$.

This is a remarkable result. Two walls, from the same simulation and with the same boundary conditions, exposed to statistically identical flows, have markedly different slip lengths. In other words, the local slip length formulation is not robust to variations in mean temperature gradient. Thus, we are motivated to utilize MFM to measure the effects of nonlocal gradients. In the MFM equations, the heating term is absorbed into the forcing, and so we can measure a slip operator which is independent of any volumetric heating.

Figure 5a shows the generalized slip length kernel as calculated by MFM. We observe that this kernel can be approximately represented as a $\delta$ function at $x_{2}=0$ and a decaying exponential that represents nonlocal effects of the SHS patterning. From this kernel, we proceed to quantify nonlocality in two ways. First, by fitting a decaying exponential to the slip length kernel (Figure 5b), we can define a length scale $\ell^{+}$ that corresponds to the rate of decay.

The decay length scale is $\ell^{+}=20.6$. We shall later compare this value to results from varying pattern geometries and momentum transport (Table 2).

The second method by which we propose to quantify nonlocality is to integrate over the slip length kernel and determine the percentage of the kernel that is influenced by local $y_{2}$. That is, how much the nonlocal velocity gradients contribute to slip velocity, compared to the contributions of local velocity gradient at the wall. For this case, the local contributions are $\int_{0}^{0^{+}}\beta^{+}_{T}dy_{2}=22.16$, the nonlocal contributions are $\int_{0^{+}}^{2h}\beta_{T}^{+}=9.01$, and the total integral is $\int_{0}^{2h}\beta^{+}_{T}dy_{2}=31.17$, for a nonlocality percentage of $29\%$ (Table 3).

We note that neither the local contributions of the slip kernel nor the total contributions of the slip kernel are equal to the measured slip length of the zero heating case, $23.3$. In Section III.1, we found that this slip length is not generalizable to nonzero heating. In contrast, the slip kernel presented here is a unified slip length model which is applicable to both zero and nonzero heating, and should give consistent results on both the top and bottom walls. We will presently investigate the performance of this nonlocal slip kernel in Section III.3.

Figure 6a shows the leading order eddy diffusivity and Figure 6b shows the full eddy diffusivity kernel for $\lambda^{+}=155.1$. A purely local eddy diffusivity kernel would only be active on the diagonal, where $x_{2}$ and $y_{2}$ values are equal. In contrast, there is significant nonlocality present in our patterned SHS case, represented by the spread in nonzero values of the kernel (Figure 7).

To assess the impact of the measured nonlocalities, we shall utilize the eddy diffusivity and slip length operators, as calculated by MFM, to solve the Reynolds-averaged scalar transport equations (Equation 7). Solutions were calculated with leading order eddy diffusivity $D_{T,12}^{0}$ and either local slip length $\beta_{T}^{0}$ or nonlocal slip kernel $\beta_{T}(y_{2})$. We would like to emphasize that the same eddy diffusivity and slip length operators are used for both the zero and nonzero heating cases.

For both cases, using only the leading order eddy diffusivity and local slip length results in large error in all regions of the domain. The maximum error between the DNS and Reynolds-averaged solutions is a difference of $0.07$ for the zero heating case and $0.26$ for the nonzero heating case, both at the boundaries. The introduction of a nonlocal slip length kernel shows dramatic improvement in the near wall region. The maximum error is reduced to $0.04$ for the zero heating case and $0.08$ for the nonzero heating.

In the zero heating case, the RMS errors are $0.03$ for the solutions with local operators and $0.01$ for the solutions with nonlocal operators. In the nonzero heating case, the RMS errors are $0.16$ for the solutions with local operators and $0.02$ for the solutions with nonlocal operators.

Our measurements of leading order eddy viscosity $D_{v,2121}^{0}$ for square posts of $\lambda^{+}=155.1$ and the corresponding $b^{+}_{\text{eff}}=20.8$ reveal that the eddy viscosity operator for SHS is distinctly different from that of a smooth wall (Figure 9), particularly in the inner zone of the flow. Specifically, the eddy viscosity for both pattern-resolved and homogenized slip SHS in this zone is larger in magnitude than that of smooth walls, indicating that high mixing is pushed closer to the wall due to the slip behavior near the surface. This enhanced mixing in the inner zone is consistent with observed extension of the log-layer towards the wall and thus shortening of the buffer layer leading to positive $\Delta U$ (momentum deficit) in the mean velocity profile.

Moreover, our results show that the pattern-resolved SHS results have substantially higher mixing, i.e., higher eddy viscosity, near the wall, with considerable momentum mixing pushed towards the viscous sublayer ($x_{2}^{+}<5$)when compared to the homogenized slip SHS. This observation is in spite of the fact that the pattern-resolved system has an equivalent effective slip length to the uniform slip case. While we expect that the usage of a homogenized slip length model in place of the more computationally expensive pattern-resolved SHS is not well-founded for our largest choice of $\lambda^{+}$, we find that this observation remains true for smaller $\lambda^{+}$ as well (Figure 9c). The significant difference in mixing behavior close to the wall may indicate that the presence of the pattern geometry is critical for all patterned SHS, regardless of the magnitude of $\lambda^{+}$.

Figure 10 shows the full eddy viscosity kernel for $\lambda^{+}=155.1$. As with the eddy diffusivity kernel for scalar transport, there is significant nonlocality. The generalized slip length kernel is plotted in Figure 11a for square posts of $\lambda^{+}=77.6$, $\lambda^{+}=103.4$, and $\lambda^{+}=155.1$, and streamwise-aligned ridges of $\lambda^{+}=155.1$. As with eddy viscosity, the trends of the slip length kernel for momentum reflect what we previously observed with scalar transport. We once more observe a local $\delta$ function and a nonlocal effect that resembles a decaying exponential.

Decay length scales are given in Table 2. We find that $\ell^{+}$ grows with larger texture size $\lambda^{+}$. That is, larger texture sizes show more nonlocal effects. This likely contributes to why the homogenized slip length model performs poorly for $\lambda^{+}>\mathcal{O}(10)$ [27].

In the limit of small pattern size, the kernel is purely local, so we expect $\ell^{+}/\lambda^{+}=0$. In the limit of large pattern size, $\ell^{+}$ is limited by the channel half height and $\lambda^{+}$ is very large, so we again expect $\ell^{+}/\lambda^{+}=0$. Thus, there is some form of nonmonotonic behavior in the ratio of $\ell^{+}/\lambda^{+}$ as $\lambda^{+}$ increases. In the square post cases, $\ell^{+}/\lambda^{+}=0.06$ for $\lambda^{+}=77.6$, $\ell^{+}/\lambda^{+}=0.08$ for $\lambda^{+}=103.4$, and $\ell^{+}/\lambda^{+}=0.09$ for $\lambda^{+}=155.1$. For the streamwise ridge case, $\ell^{+}/\lambda^{+}=0.09$.

The nonlocality percentages of the momentum transport cases can be found in Table 3. Again, we find that larger texture size $\lambda^{+}$ corresponds to more significant effects arising from nonlocality. The nonlocality percentage is noticeably affected by the difference between pattern geometries of square posts and streamwise-aligned ridges. This behavior is perhaps predictable from close examination of Figure 11. The kernel associated with streamwise-aligned ridges has a distinctly different shape from those of the square posts, especially near the wall. The decay length scale, which is measured from kernel contributions further from the wall, seems to depend more upon $\lambda^{+}$ than on geometry.

We now present a direct comparison of our results from assessing scalar and momentum transport of the $\lambda^{+}=155.1$ square post case. The leading order eddy diffusivity from scalar transport and the leading order eddy viscosity from momentum transport are shown in Figure 12. The nonlocal eddy diffusivity and eddy viscosity kernels are compared in Figure 13. The eddy viscosity of momentum transport has a narrower band of nonlocality relative to the eddy diffusivity of scalar transport.

Figure 14 shows the nonlocal slip kernels of scalar and momentum transport. The slip kernel associated with scalar transport is visibly more nonlocal, as compared to the kernel associated with momentum transport. Decay length scales and nonlocality percentage are given in Tables 2 and 3, and confirm this observation: the decay length scale and nonlocality percentage of the scalar result are both significantly larger in magnitude than the momentum result.

This may, at first glance, be an unintuitive result. The scalar transport equation is completely local, while the momentum transport equation involves a nonlocal pressure term. However, this finding is consistent with previous MFM analysis of homogeneous isotropic turbulence [42]. Recall that in our homogenized setting we have averaged in the streamwise and spanwise directions, as well as in time. Thus, advection has also become a nonlocal effect in this framework. Consistent with Shirian and Mani (2022) [42], we conclude that the nonlocal advection effects are more prominent in scalar transport, or that nonlocal advection and pressure may partially cancel in momentum transport.

We now use the local and nonlocal eddy viscosity and slip length operators to calculate the solution to the RANS equations (Equation 6) for pattern size $\lambda^{+}=155.1$. Solutions were calculated with leading order eddy viscosity $D_{v,2121}^{0}$ and local slip length $\beta_{v}^{0}$, leading order eddy viscosity $D_{v,2121}^{0}$ and nonlocal slip length $\beta_{v}(y_{2})$, and nonlocal eddy viscosity $D_{v,2121}(x_{2},y_{2})$ and nonlocal slip length $\beta_{v}(y_{2})$.

Using leading order eddy viscosity and local slip length results in the largest error between DNS and RANS solutions. The maximum error is a difference of 4.5 at the boundaries. The incorporation of nonlocal slip length kernel, while still using leading order eddy viscosity, greatly improves performance in the near wall region. The maximum error is 1.1 and the location is shifted to the buffer layer ($x_{2}^{+}\approx 20$). The nonlocal eddy viscosity operator further increases the accuracy of the RANS calculation; the maximum error is 0.6. The RMS errors are 3.5 for leading order eddy viscosity and local slip length, 0.9 for leading order eddy viscosity and nonlocal slip length, and 0.2 for nonlocal eddy viscosity and nonlocal slip length.

The improvement in the RANS solutions due to the inclusion of nonlocality in the eddy viscosity and slip length operators is particularly evident for the $\lambda^{+}=155.1$ case, as we have previously measured the decay length scale and nonlocality percentage to be significant. For the smaller pattern sizes of $\lambda^{+}=77.6$ and $\lambda^{+}=103.4$, the measured influence of using the nonlocal operators is on the order of numerical and statistical error.