Search

Superhydrophobic surfaces (SHS) are textured hydrophobic surfaces which have the ability to trap air pockets when immersed in water. This can result in significant drag reduction, due to substantially lower viscosity of air resulting in substantial effective slip velocity at the interface. Past studies of both laminar and turbulent flows model this slip velocity in terms of a homogenized Navier slip boundary condition with a slip length relating the wall slip velocity to the wall-normal velocity gradient. In this work, we seek to understand the effects of superhydrophobic surfaces in the context of mean scalar and momentum mixing. We use the macroscopic forcing method (Mani and Park, 2021) to compute the generalized eddy viscosity and slip length operators of a turbulent channel over SHS, implemented as both a pattern-resolved boundary condition and homogenized slip length boundary condition, for several pattern sizes and geometries. We present key differences in the mixing behavior of both boundary conditions through quantification of their near-wall eddy viscosity. Analysis of transport in turbulent flow over pattern-resolved surfaces reveals substantial nonlocality in teh measured homogenized slip length for both scalar and momentum mixing when the Reynolds and Peclet numbers based on pattern size are finite. We present several metrics to quantify this nonlocality and observe possible trends relating to Reynolds number, texture size, and pattern geometry. The importance of nonlocality in the slip length operator and in the eddy diffusivity operator is demonstrated by examining impact on Reynolds-averaged solutions for the mean scalar and velocity fields.

Models for solving the Reynolds-averaged Navier-Stokes equations are popular tools for predicting complex turbulent flows due to their computational affordability and ability to provide or estimate quantities of engineering interest. However, results depend on a proper treatment of unclosed terms, which require progress in the development and assessment of model forms. In this study, we consider the Reynolds stress transport equations as a framework for second-moment turbulence closure modeling. We specifically focus on the terms responsible for decay of the Reynolds stresses, which can be isolated and evaluated separately from other terms in a canonical setup of homogeneous turbulence. We show that by using anisotropic forcing of the momentum equation, we can access states of turbulence traditionally not probed in a triply-periodic domain. The resulting data span a wide range of anisotropic turbulent behavior in a more comprehensive manner than extant literature. We then consider a variety of model forms for which these data allow us to perform a robust selection of model coefficients and select an optimal model that extends to cubic terms when expressed in terms of the principal coordinate Reynolds stresses. Performance of the selected decay model is then examined relative to the simulation data and popular models from the literature, demonstrating the superior accuracy of the developed model and, in turn, the efficacy of this framework for model selection and tuning.

Low Stokes number particles at dilute concentrations in turbulent flows can reasonably be approximated as passive scalars. The added presence of a drift velocity due to buoyancy or gravity when considering the transport of such passive scalars can reduce the turbulent dispersion of the scalar via a diminution of the eddy diffusivity. In this work, we propose a model to describe this decay and use a recently developed technique to accurately and efficiently measure the eddy diffusivity using Eulerian fields and quantities. We then show a correspondence between this method and standard Lagrangian definitions of diffusivity and collect data across a range of drift velocities and Reynolds numbers. The proposed model agrees with data from these direct numerical simulations, offers some improvement to previous models in describing other computational and experimental data, and satisfies theoretical constraints that are independent of Reynolds number.

Reynolds-averaged Navier--Stokes (RANS) closure must be sensitive to the flow physics, including nonlocality and anisotropy of the effective eddy viscosity. Recent approaches used forced direct numerical simulations to probe these effects, including the macroscopic forcing method (MFM) of Mani and Park ($\textit{Phys. Rev. Fluids}$ $\textbf{6}$, 054607 (2021)) and the Green's function approach of Hamba ($\textit{Phys. Fluids}$ $\textbf{17}$, 115102 (2005)). The resulting nonlocal and anisotropic eddy viscosities are exact and relate Reynolds stresses to mean velocity gradients at all locations. They can be used to inform RANS models of the sensitivity to the mean velocity gradient and the suitability of local and isotropic approximations. However, these brute-force approaches are expensive. They force the mean velocity gradient at each point in the averaged space and measure the Reynolds stress response, requiring a separate simulation for each mean velocity gradient location. Thus, computing the eddy viscosity requires as many simulations as degrees of freedom in the averaged space, which can be cost-prohibitive for problems with many degrees of freedom. In this work, we develop an adjoint-based MFM to obtain the eddy viscosity at a given Reynolds stress location using a single simulation. This approach recovers the Reynolds stress dependence at a location of interest, such as a separation point or near a wall, on the mean velocity gradient at all locations. We demonstrate using adjoint MFM to compute the eddy viscosity for a specified wall-normal location in an incompressible turbulent channel flow using one simulation. In contrast, a brute-force approach for the same problem requires $N=144$ simulations (the number of grid points in the non-averaged coordinate direction). We show that a local approximation for the eddy viscosity would have been inappropriate.

The importance of nonlocality of mean scalar transport in 2D Rayleigh-Taylor Instability (RTI) is investigated. The Macroscopic Forcing Method (MFM) is utilized to measure spatio-temporal moments of the eddy diffusivity kernel representing passive scalar transport in the ensemble averaged fields. Presented in this work are several studies assessing the importance of the higher-order moments of the eddy diffusivity, which contain information about nonlocality, in models for RTI. First, it is demonstrated through a comparison of leading-order models that a purely local eddy diffusivity is insufficient in capturing the mean field evolution of the mass fraction in RTI. Therefore, higher-order moments of the eddy diffusivity operator are not negligible. Models are then constructed by utilizing the measured higher-order moments. It is demonstrated that an explicit operator based on the Kramers-Moyal expansion of the eddy diffusivity kernel is insufficient. An implicit operator construction that matches the measured moments is shown to offer improvements relative to the local model in a converging fashion.

The macroscopic forcing method (MFM) of Mani and Park and similar methods for obtaining turbulence closure operators, such as the Green's function-based approach of Hamba, recover reduced solution operators from repeated direct numerical simulations (DNS). MFM has been used to quantify RANS-like operators for homogeneous isotropic turbulence and turbulent channel flows. Standard algorithms for MFM force each coarse-scale degree of freedom (i.e., degree of freedom in the RANS space) and conduct a corresponding fine-scale simulation (i.e., DNS), which is expensive. We combine this method with an approach recently proposed by Sch\"afer and Owhadi (2023) to recover elliptic integral operators from a polylogarithmic number of matrix-vector products. The resulting Fast MFM introduced in this work applies sparse reconstruction to expose local features in the closure operator and reconstructs this coarse-grained differential operator in only a few matrix-vector products and correspondingly, a few MFM simulations. For flows with significant nonlocality, the algorithm first "peels" long-range effects with dense matrix-vector products to expose a local operator. We demonstrate the algorithm's performance for scalar transport in a laminar channel flow and momentum transport in a turbulent one. For these, we recover eddy diffusivity operators at 1% of the cost of computing the exact operator via a brute-force approach for the laminar channel flow problem and 13% for the turbulent one. We observe that we can reconstruct these operators with an increase in accuracy by about a factor of 100 over randomized low-rank methods. We glean that for problems in which the RANS space is reducible to one dimension, eddy diffusivity and eddy viscosity operators can be reconstructed with reasonable accuracy using only a few simulations, regardless of simulation resolution or degrees of freedom., Comment: 16 pages, 10 figures. S. H. Bryngelson and F. Sch\"afer contributed equally to this work

A major challenge in developing accurate and robust numerical solutions to multi-physics problems is to correctly model evolving discontinuities in field quantities, which manifest themselves as interfaces between different phases in multi-phase flows, or as shock and contact discontinuities in compressible flows. When a quick response is required to rapidly emerging challenges, the complexity of bespoke discretization schemes impedes a swift transition from problem formulation to computation, which is exacerbated by the need to compose multiple interacting physics. We introduce "inverse asymptotic treatment" (IAT) as a unified framework for capturing discontinuities in fluid flows that enables building directly computable models based on off-the-shelf numerics. By capturing discontinuities through modifications at the level of the governing equations, IAT can seamlessly handle additional physics and thus enable novice end users to quickly obtain numerical results for various multi-physics scenarios. We outline IAT in the context of phase-field modeling of two-phase incompressible flows, and then demonstrate its generality by showing how localized artificial diffusivity (LAD) methods for single-phase compressible flows can be viewed as instances of IAT. Through the real-world example of a laminar hypersonic compression corner, we illustrate IAT's ability to, within just a few months, generate a directly computable model whose wall metrics predictions for sufficiently small corner angles come close to that of NASA's VULKAN-CFD solver. Finally, we propose a novel LAD approach via "reverse-engineered" PDE modifications, inspired by total variation diminishing (TVD) flux limiters, to eliminate the problem-dependent parameter tuning that plagues traditional LAD. We demonstrate that, when combined with second-order central differencing, it can robustly and accurately model compressible flows.

Computational modeling of charged species transport has enabled the analysis, design, and optimization of a diverse array of electrochemical and electrokinetic devices. These systems are represented by the Poisson-Nernst-Planck (PNP) equations coupled with the Navier-Stokes (NS) equation. Direct numerical simulation (DNS) to accurately capture the spatio-temporal variation of ion concentration and current flux remains challenging due to the (a) small critical dimension of the diffuse charge layer (DCL), (b) stiff coupling due to fast charge relaxation times, large advective effects, and steep gradients close to boundaries, and (c) complex geometries exhibited by electrochemical devices. In the current study, we address these challenges by presenting a direct numerical simulation framework that incorporates (a) a variational multiscale (VMS) treatment, (b) a block-iterative strategy in conjunction with semi-implicit (for NS) and implicit (for PNP) time integrators, and (c) octree based adaptive mesh refinement. The VMS formulation provides numerical stabilization critical for capturing the electro-convective flows often observed in engineered devices. The block-iterative strategy decouples the difficulty of non-linear coupling between the NS and PNP equations and allows the use of tailored numerical schemes separately for NS and PNP equations. The carefully designed second-order, hybrid implicit methods circumvent the harsh timestep requirements of explicit time steppers, thus enabling simulations over longer time horizons. Finally, the octree-based meshing allows efficient and targeted spatial resolution of the DCL. These features are incorporated into a massively parallel computational framework, enabling the simulation of realistic engineering electrochemical devices. The numerical framework is illustrated using several challenging canonical examples., Comment: Accepted for publication in Journal of Computational Physics, 34 pages, 12 figures

Second-order phase field models have emerged as an attractive option for capturing the advection of interfaces in two-phase flows. Prior to these, state-of-the-art models based on the Cahn-Hilliard equation, which is a fourth-order equation, allowed for the derivation of surface tension models through thermodynamic arguments. In contrast, the second-order phase field models do not follow a known energy law, and deriving a surface tension term for these models using thermodynamic arguments is not straightforward. In this work, we justify that the energy-based surface tension model from the Cahn-Hilliard context can be adopted for second-order phase field models as well and assess its performance. We test the surface tension model on three different second-order phase field equations; the conservative diffuse interface model of Chiu and Lin [1], and two models based on the modified Allen-Cahn equation introduced by Sun and Beckermann [2]. Additionally, we draw the connection between the energy-based model with a localized variation of the continuum surface force (CSF) model. Using canonical tests, we illustrate the lower magnitude of spurious currents, better accuracy, and superior convergence properties of the energy-based surface tension model compared to the CSF model, which is a popular choice used in conjunction with second-order phase field methods, and the localized CSF model. Importantly, in terms of computational expense and parallel efficiency, the energy-based model incurs no penalty compared to the CSF models., Comment: 13 pages, 5 figures, Revision submitted to Journal of Computational Physics

A decade ago, two-dimensional microscopic flow visualization proved the theoretically predicted existence of electroconvection roles as well as their decisive role in destabilizing the concentration polarization layer at ion-selective fluid/membrane interfaces. Electroconvection induces chaotic flow vortices injecting volume having bulk concentration into the ion-depleted diffusion layer at the interface. Experimental quantification of these important flow patterns have so far only been carried out in 2D. Numerical direct simulations suggest 3D features, yet experimental proof is lacking. 3D simulations are also limited in covering extended spacial and temporal scales. This study presents a new comprehensive experimental method for the time-resolved recording of the 3D electroconvective velocity field near a cation-exchange membrane. For the first time, the spatio-temporal velocity field can be visualized in 3D at multiples of the overlimiting current density. In contrast to today's simulations, these experiments cover length and time scales typical for actual electrodialytic membrane processes. We visualize coherent vortex structures and reveal the changes in the velocity field and its statistics during the transition from vortex rolls to vortex rings with increasing current density. The transition is characterized by changes in the rotational direction, mean square velocity, and temporal energy spectrum with only little influence on the spatial spectrum. These findings indicate a more significant impact of EC's structural change on the mean square velocities and temporal spectra than on the spatial spectra. This knowledge is a prerequisite for engineering ion-selective surfaces that will enable the operation of electrically driven processes beyond the diffusion-limited Nernst regime.

Algebraic closure models with spatially nonlocal operators that are associated with both unresolved advective transport and nonlinear reaction terms in a Reynolds-averaged Navier-Stokes context are presented in this work. In particular, a system of two species subject to binary reaction and transport by advection and diffusion are examined by expanding upon analysis originally developed for binary reactions in the context of Taylor dispersion of scalars. This work extends model forms from weakly-nonlinear extensions of that dispersion theory and the role of nonlocality in the presence of reactions is studied and captured by analytic expressions. These expressions can be incorporated into an eddy diffusivity matrix that explicitly capture the influence of chemical kinetics and flow conditions on the closure operators and we demonstrate that the model form derived in a laminar context can be directly translated to an analogous setup in homogeneous isotropic turbulence, which has implications as a subgrid scale model. We show that this framework can improve prediction of mean quantities compared to previous purely local results, but does not fully close unresolved terms., Comment: arXiv admin note: text overlap with arXiv:2111.10551

This work presents algebraic closure models associated with advective transport and nonlinear reactions in a Reynolds-averaged Navier-Stokes context for a system of species subject to binary reactions and transport by advection and diffusion. Expanding upon analysis originally developed for non-reactive transport in the context of Taylor dispersion of scalars, this work extends the modified gradient diffusion model explicated by Peters (Turbulent Combustion, 2000) and based on work by Corrsin (JFM, vol. 11, p.407-416) beyond single-component transport phenomena and involving nonlinear reactions. The presented model forms, from this weakly-nonlinear extension of the original dispersion theory, lead to an analytic expression for the eddy diffusivity matrix that explicitly captures the influence of the reaction kinetics on the closure operators. Furthermore, we demonstrate that the derived model form directly translates between flow topologies through a priori and a posteriori testing of a binary species system subject to homogeneous isotropic turbulence. Using two- and three-dimensional direct numerical simulations involving laminar and turbulent flows, it is shown that this framework improves prediction of mean quantities compared to previous results. Lastly, the presented model form, collapses to the earlier gradient diffusion and its modified version derived by Corrsin in the limits of non-reactive species and linear reactions, respectively.

This study considers advective and diffusive transport of passive scalar fields by spatially-varying incompressible flows. Prior studies have shown that the eddy diffusivities governing the mean field transport in such systems can generally be nonlocal in space and time. While for many flows nonlocal eddy diffusivities are more accurate than commonly-used Boussinesq eddy diffusivities, nonlocal eddy diffusivities are often computationally cost-prohibitive to obtain and difficult to implement in practice. We develop a systematic and more cost-effective approach for modeling nonlocal eddy diffusivities using matched moment inverse (MMI) operators. These operators are constructed using only a few leading-order moments of the exact nonlocal eddy diffusivity kernel, which can be easily computed using the inverse macroscopic forcing method (IMFM) (Mani and Park (2021)). The resulting reduced-order models for the mean fields that incorporate the modeled eddy diffusivities often improve Boussinesq-limit models since they capture leading-order nonlocal effects. But more importantly, these models can be expressed as partial differential equations that are readily solvable using existing computational fluid dynamics capabilities rather than as integro-partial differential equations.

When reporting statistics from simulations of statistically stationary chaotic phenomenon, it is important to verify that the simulations are time-converged. This condition is connected with the statistical error or number of digits with which statistics can be reliably reported. In this work we consider homogeneous and isotropic turbulence as a model problem to investigate statistical convergence over finite simulation times. Specifically, we investigate the time integration requirements that allow meaningful reporting of the statistical error associated with finiteness of the temporal domain. We address two key questions: 1) What is the appropriate range of sampling frequency in large eddy time units? and 2) How long should a simulation be performed in terms of large eddy time so that the statistical error could be reliably reported. Our results indicate that proper sampling frequency is on the order of 10 large eddy time scale. More importantly, we find that reliable reporting of statistical errors requires simulation durations orders of magnitude longer than typically performed. Our observations of sampling frequency for homogeneous isotropic turbulence are also shown to hold in turbulent channel flow.

This work describes three diffuse-interface methods for the simulation of immiscible, compressible multiphase fluid flows and elastic-plastic deformation in solids. The first method is the localized-artificial-diffusivity approach of Cook (2007), Subramaniam et al., (2018), and Adler and Lele (2019), in which artificial diffusion terms are added to the conservation equations. The second method is the gradient-form approach that is based on the quasi-conservative method of Shukla et al., (2010) and Tiwari et al., (2013), in which the diffusion and sharpening terms (together called regularization terms) are added to the conservation equations. The third approach is the divergence-form approach that is based on the fully conservative method of Jain et al., (2020), in which the regularization terms are added to the conservation equations. The primary objective of this work is to compare these three methods in terms of their ability to maintain: constant interface thickness; the conservation property; and accurate interface shape for long-time integration. The second objective of this work is to consistently extend these methods to model interfaces between solid materials with strength. To assess and compare the methods, they are used to simulate a wide variety of problems, including (1) advection of an air bubble in water, (2) shock interaction with a helium bubble in air, (3) shock interaction and the collapse of an air bubble in water, and (4) Richtmyer-Meshkov instability of a copper-aluminum interface. The current work focuses on comparing these methods in the limit of relatively coarse grid resolution, which illustrates the true performance of these methods. This is because it is rarely practical to use hundreds of grid points to resolve a single bubble or drop in large-scale simulations of engineering interest., Comment: 37 pages, 12 figures, submission for journal of computational physics

This study aims to quantify how turbulence in a channel flow mixes momentum in the mean sense. We applied the macroscopic forcing method (Mani and Park, Physical Review Fluids, 2021, p.054607) to direct numerical simulation (DNS) of a turbulent channel flow at Re_tau=180 using two different forcing strategies that are designed to separately assess the anisotropy and nonlocality of momentum mixing. In the first strategy, the leading term of the Kramers-Moyal expansion of the eddy viscosity is quantified, revealing all 81 tensorial coefficients that essentially characterize the local-limit eddy viscosity. The results indicate: (1) the eddy viscosity has significant anisotropy, (2) Reynolds stresses are generated by both the mean strain rate and mean rotation rate tensors associated with the momentum field, and (3) the local-limit eddy viscosity generates asymmetric Reynolds stress tensors. In the second strategy, the eddy viscosity is quantified as an integration kernel revealing the nonlocal influence of the mean momentum gradient at each wall-normal coordinate on all nine components of the Reynolds stresses over the channel width. Our results indicate that while the shear component of the Reynolds stress is reasonably controlled by the local mean gradients, other components of the Reynolds stress are highly nonlocal. These results provide an understanding of anisotropy and nonlocality requirements for closure modeling of momentum transport in wall-bounded turbulent flows., Comment: 31 pages, 24 figures

In this article, we propose a novel scalar-transport model for the simulation of scalar quantities in two-phase flows with a phase-field method (diffuse-interface method). In a two-phase flow, the scalar quantities typically have disparate properties in two phases, which results in effective confinement of the scalar quantities in one of the phases, in the time scales of interest. This confinement of the scalars lead to the formation of sharp gradients of the scalar concentration values at the interface, presenting a serious challenge for its numerical simulations. To overcome this challenge, we propose a model for the transport of scalars. The model is discretized using a central-difference scheme, which leads to a non-dissipative implementation that is crucial for the simulation of turbulent flows. Furthermore, the provable strengths of the proposed model are: (a) the model maintains the positivity property of the scalar concentration field, a physical realizability requirement for the simulation of scalars, when the proposed criterion is satisfied, (b) the proposed model is such that the transport of the scalar concentration field is consistent with the transport of the volume fraction field, which results in the enforcement of the effective zero-flux boundary condition for the scalar at the interface; and therefore, prevents the artificial numerical diffusion of the scalar across the interface. Finally, we present numerical simulations using the proposed model in a wide range of two-phase flow regimes, spanning laminar to turbulent flows; and assess: the accuracy and robustness of the model, the validity of the positivity property of the scalar concentration field, and the enforcement of the zero-flux boundary condition for the scalar at the interface., Comment: 38 pages, 24 figures, submitted to the journal of computational physics

Shock electrodialysis is a recently developed electrochemical water treatment method which shows promise for water deionization and ionic separations. Although simple models and scaling laws have been proposed, a predictive theory has not yet emerged to fit experimental data and enable system design. Here, we extend and analyze existing "leaky membrane" models for the canonical case of a steady shock in cross flow, as in recent experimental prototypes. Two-dimensional numerical solutions are compared with analytical boundary-layer approximations and experimental data. The boundary-layer theory accurately reproduces the simulation results for desalination, and both models predict the data collapse of the desalination factor with dimensionless current, scaled to the incoming convective flux of cations. The numerical simulation also predicts the water recovery increase with current. Nevertheless, both approaches cannot quantitatively fit the transition from normal to over-limiting current, which suggests gaps in our understanding of extreme electrokinetic phenomena in porous media.

Realistic two-phase flow problems of interest often involve high $Re$ flows with high density ratios. Accurate and robust simulation of such problems requires special treatments. In this work, we present a consistent, energy-conserving momentum transport scheme in the context of a second order mass-conserving phase field method. This is achieved by (1) accounting for the mass flux associated with the right-hand-side of the phase field equation in the convective flux of the conservative form of the momentum transport equation---a correction absent in previous phase field simulations (2) utilization of non-dissipative spatial discretization. We demonstrate accuracy and robustness improvements from our proposed scheme via numerical tests, including a turbulent case of a water jet subject to air cross-flow. Our proposed modifications to the momentum transport equation can be quite readily extended to general conservative phase field methods. Additionally, in the framework of this phase field method, we present a free energy-based surface tension force calculation scheme. This scheme, which significantly reduces spurious currents, is based on a general paradigm that can also be extended to other two-phase flow methods.

Unconditionally stable implicit time-marching methods are powerful in solving stiff differential equations efficiently. In this work, a novel framework to handle stiff physical terms implicitly is proposed. Both physical and numerical stiffness originating from convection, diffusion and source terms (typically related to reaction) can be handled by a set of predefined Time-Accurate and highly-Stable Explicit (TASE) operators in a unified framework. The proposed TASE operators act as preconditioners on the stiff terms and can be deployed to any existing explicit time-marching methods straightforwardly. The resulting time integration methods remain the original explicit time-marching schemes, yet with nearly unconditional stability. The TASE operators can be designed to be arbitrarily high-order accurate with Richardson extrapolation such that the accuracy order of original explicit time-marching method is preserved. Theoretical analyses and stability diagrams show that the $s$-stages $s$th-order explicit Runge-Kutta (RK) methods are unconditionally stable when preconditioned by the TASE operators with order $p \leq s$ and $p \leq 2$. On the other hand, the $s$th-order RK methods preconditioned by the TASE operators with order of $p \leq s$ and $p > 2$ are nearly unconditionally stable. The only free parameter in TASE operators can be determined a priori based on stability arguments. Unlike classical implicit methods, the TASE methodology allows for solving non-linear problems with arbitrary order without requiring solving a nonlinear system of equations. A set of benchmark problems with strong stiffness is simulated to assess the performance of the TASE method. Numerical results suggest that the proposed framework preserves the high-order accuracy of the explicit time-marching methods with very-large time steps for all the considered cases.

In this article, we propose a novel conservative diffuse-interface method for the simulation of immiscible compressible two-phase flows. The proposed method discretely conserves the mass of each phase, momentum and total energy of the system. We use the baseline five-equation model and propose interface-regularization (diffusion--sharpening) terms in such a way that the resulting model maintains the conservative property of the underlying baseline model; and lets us use a central-difference scheme for the discretization of all the operators in the model, which leads to a non-dissipative implementation that is crucial for the simulation of turbulent flows and acoustics. Furthermore, the provable strengths of the proposed model are: (a) the model maintains the boundedness property of the volume fraction field, which is a physical realizability requirement for the simulation of two-phase flows, (b) the proposed model is such that the transport of volume fraction field inherently satisfies the total-variation-diminishing property without having to add any flux limiters that destroy the non-dissipative nature of the scheme, (c) the proposed interface-regularization terms in the model do not spuriously contribute to the kinetic energy of the system and therefore do not affect the non-linear stability of the numerical simulation, and (d) the model is consistent with the second law of thermodynamics. Finally, we present numerical simulations using the model and assess (a) the accuracy of evolution of the interface shape, (b) implementation of surface tension effects, (c) propagation of acoustics and their interaction with material interfaces, (d) the accuracy and robustness of the numerical scheme for simulation of complex high-Reynolds-number flows, and (e) performance and scalability of the method., Comment: 50 pages, 24 figures

Ion transport in perm-selective media has been extensively investigated in recent years owing to its applications in advancements of fundamental understanding of nanoscale electrokinetics as well as innovative engineering applications. A key phenomenon in this context is ion concentration polarization (ICP) that occurs near perm-selective nanoporous membranes or nanochannels under dc bias. In classical settings, involving voltage-driven transport of ions through perm-selective membranes, concentration polarization is well understood as formation of steep concentration drop to reach near zero concentration (depletion layer) at the anodic side of cation-selective membranes. In contrast to this classical description, we demonstrate here that when ICP is driven in a microchannel and coupled with electroosmotic flows a long-tailed fluorescent layer in front of diffusion layer is formed, which can be characterized as an additional plateau in the concentration profile. Using a micro/nanofluidic platform, in situ visualizations, and local concentration measurements by microelectrode arrays we quantify that this plateau layer has a concentration of 60 % of bulk electrolyte concentration, regardless of the bulk concentration itself. Using a mathematical analysis of conservation laws for flow, charge, and salt, we demonstrate that this plateau region is related to the formation of a flow stagnation point due to the competition between electro-osmotic and induced pressure driven flows. Consistent with the experiments, this analysis predicts the plateau concentration to be 60% of the bulk electrolyte. Thus, this plateau concentration region completes the picture of ICP in microchannels which has drawn significant attention from the fields of electrokinetics and micro/nanofluidics.

Overlimiting current (OLC) through electrolytes interfaced with perm-selective membranes has been extensively researched in recent years for understanding the fundamental mechanisms of transport and developing efficient applications from electrochemistry to sample analysis and separation. Predominant mechanisms responsible for OLC include surface conduction, convection by electro-osmotic flow, and electro-osmotic instability depending on input parameters such as surface charge and geometric constrictions. This work studies how a network of microchannels in a non-uniform array, which mimicks a natural pore configuration, can contribute to OLC. To this end, micro/nanofluidic devices are fabricated with arrays of parallel microchannels with non-uniform size distributions. All cases maintain the same surface and bulk conduction to allow probing the sensitivity only by the non-uniformity of the channels. Both experimental and theoretical current-voltage relations demonstrate that OLCs increase with increasing non-uniformity. Furthermore, the visualization of internal recirculating flows indicates that the non-uniform arrays induce flow loops across the network enhancing advective transport. These evidences confirm a new driving mechanism of OLC, inspired by natural micro/nanoporous materials with random geometric structure. Therefore, this result can advance not only the fundamental understanding of nanoelectrokinetics but also the design rule of engineering applications of electrochemical membrane.

We use the recently developed Macroscopic Forcing Method [Mani and Park, Physical Review Fluids, 6:054607, 2021] to compute the scale-dependent eddy diffusivity characterizing ensemble-averaged scalar and momentum transport in incompressible homogeneous isotropic turbulence. For scales larger than the energy containing eddies, eddy diffusivity is found to be constant and consistent with the Boussinesq approximation. However for small scales eddy diffusivity is found to vanish inversely proportional to the wavenumber. Behavior at all scales is reasonably captured by a non-local eddy diffusivity operator modeled as $D/\sqrt{\mathcal{I}-l^2\nabla^2}$, where $D$ is the eddy diffusivity in the Boussinesq limit, and $l$ is a constant on the order of the large-eddy length. These results present the first direct measurement of eddy diffusivity in turbulence with implications in turbulence modeling.

This study presents a numerical procedure, which we call the macroscopic forcing method (MFM), which reveals the differential operators acting on the mean fields of quantities transported by underlying fluctuating flows. Specifically, MFM can precisely determine the eddy diffusivity operator, or more broadly said, it can reveal differential operators associated with turbulence closure for scalar and momentum transport. We present this methodology by considering canonical problems with increasing complexity. Starting from the well-known problem of dispersion of passive scalars by parallel flows we elucidate the basic steps in quantitative determination of the eddy viscosity operators using MFM. Utilizing the operator representation in Fourier space, we obtain a stand-alone compact analytical operator that can accurately capture the nonlocal mixing effects. Furthermore, a cost-effective generalization of MFM for analysis of nonhomogeneous and wall-bounded flows is developed and is comprehensively discussed through a demonstrative example. Extension of MFM for analysis of momentum transport is theoretically constructed through the introduction of a generalized momentum transport equation. We show that closure operators obtained through MFM analysis of this equation provide the exact RANS solutions obtained through averaging of the Navier-Stokes equation. We introduce MFM as an effective tool for quantitative understanding of non-Boussinesq effects and assessment of model forms in turbulence closures, particularly, the effects associated with anisotropy and nonlocality of macroscopic mixing., Comment: 44 pages, 9 figures

Soft solids in fluids find wide range of applications in science and engineering, especially in the study of biological tissues and membranes. In this study, an Eulerian finite volume approach has been developed to simulate fully resolved incompressible hyperelastic solids immersed in a fluid. We have adopted the recently developed reference map technique (RMT) by Valkov et. al (J. Appl. Mech., 82, 2015) and assessed multiple improvements for this approach.These modifications maintain the numerical robustness of the solver and allow the simulations without any artificial viscosity in the solid regions (to stabilize the solver). This has also resulted in eliminating the striations ("wrinkles") of the fluid-solid interface that was seen before and hence obviates the need for any additional routines to achieve a smooth interface. An approximate projection method has been used to project the velocity field onto a divergence free field. Cost and accuracy improvements of the modifications on the method have also been discussed.

We present a numerical study of nonlinear electrokinetic transport in porous media, focusing on the role of heterogeneity in a porous microstructure on ion concentration polarization and over-limiting current. For simplicity, the porous medium is modeled as a network of long, thin charged cylindrical pores, each governed by one-dimensional effective transport equations. For weak surface conduction, when sufficiently large potential is applied, we demonstrate that electrokinetic transport in a porous network can be dominated by electroconvection via internally induced flow loops, which is not properly captured by existing homogenized models. We systematically vary the topology and "accessivity" of the pore network and compare with simulations of traditional homogenized parallel-pore (capillary-bundle) models, in order to reveal the effects of regular and hierarchical connectivity. Our computational framework sheds light on the complex physics of nonlinear electrokinetic phenomena in microstructures and may be used to design porous media for applications, such as water desalination and purification by shock electrodialysis., Comment: 41 pages, 15 page supporting Information Material

Micro-bubbles are of significant interest due to the long-living signature they leave behind naval ships. In order to numerically model and predict these bubbles in naval applications, subgrid-scale models are required because of the extreme separation of length- and time-scales between the macroscopic geometries and the physical processes leading to the formation of these bubbles. Yet, there is much that is unknown about the mechanism behind the entrainment of such bubbles. Furthermore, quantitative information regarding their size distribution and dependence on flow parameters is very limited. Impact events are hypothesized to be the main contributor to the generation of micro-bubbles. This is due to a phenomenon known as Mesler entrainment, which has been observed in the context of the drop-pool impact problem. Namely, when a water droplet with diameter of $\mathcal{O}(1mm)$ impacts a deep water pool with an impact velocity of $\mathcal{O}(1m/s)$, hundreds of air micro-bubbles are entrained into the pool. These bubbles have been found to be remnants of a very thin air film entrapped between the two liquid bodies. These films have extremely high aspect ratios and after being punctured, shed micro-bubbles while retracting on time scales much shorter than the outer flow time scales. This separation of time-scales, along with a lacking of studies on retracting thin gas films has motivated us to study this fundamental problem with numerical simulations in two and three dimensions. Using a diffuse interface method, we perform two-phase simulations of retracting thin gas films in initially static liquid backgrounds to gain understanding regarding this problem and gather statistics that can be potentially used in a subgrid-scale model to predict micro-bubble shedding from a thin gas film. Coupling a model like this to a model that predicts the thin film characteristics ..., Comment: Proceedings of 32nd Symposium on Naval Hydrodynamics, Hamburg, Germany, 5-10 Aug. 2018

A wide variety of interface capturing methods have been introduced for simulating two-phase flows throughout the years. However, there is a noticeable dearth of literature focusing on objective comparisons between these methods, especially when they are coupled to the momentum equation and applied in physically relevant regimes. In this article, we compare two techniques for simulating two-phase flows that possess attractive qualities, but belong to the two distinct classes of diffuse interface (DI) and volume of fluid (VOF) methods. Both of these methods allow for mass-conserving schemes that can naturally capture large interfacial topology changes omnipresent in realistic two phase flows. The DI solver used in this work is based on a conservative and bounded phase field method, developed recently. Similar to level set methods, this diffuse interface method takes advantage of the smoothness of the phase field in computing curvature and surface tension forces. Geometric VOF methods track the fractional tagged volume in a cell. The specific geometric VOF scheme used here is a discretely conservative and bounded implementation that uses geometric algorithms for unsplit advection and interface reconstruction, while employing height functions for normal and curvature calculation. We present a quantitative comparison of these methods on Cartesian meshes in terms of their accuracy, convergence rate, and computational cost using canonical two-dimensional (2D) two-phase test cases: a very dense drop moving through a quiescent gas, the Rayleigh-Taylor instability, an equilibrium static drop, an oscillating drop and the damped surface wave. We further compare these methods in their ability to resolve thin films by simulating the impact of a water drop on a deep water pool. Using results of these studies, we suggest qualitative guidelines for selection of schemes for two-phase flow calculations., Comment: Under review for publication in Journal of Computational Physics

Central finite difference schemes have long been avoided in the context of two-phase flows for the advection of the phase indicator function due to numerical overshoots and undershoots associated with their dispersion errors. We will show however, for an incompressible flow, in the context of a specific diffuse interface model, one can maintain the boundedness of the phase field while also taking advantage of the low cost and ease of implementation of central differences to construct a non-dissipative discretization scheme for the advective terms. This is made possible by combining the advection and reinitialization steps of a conservative level set scheme introduced by Olsson and Kreiss [J. Comput. Phys., 210, 225 (2005)] to form a phase field equation similar to that of Chiu and Lin [J. Comput. Phys., 230, 185 (2011)]. Instead of resorting to specialized upwind methods as in these articles, we prove that the boundedness of the phase field is guaranteed for certain choices of the free parameters ($\epsilon$ and $\gamma$) for a specific central difference scheme that we propose. The proposed discretely conservative and bounded phase field equation, which is free of any reinitialization or mass redistribution, possesses desirable properties that can be leveraged in the coupled finite difference discretization of the two-phase momentum equation. Additionally, as compared to the state-of-the-art conservative and bounded two-phase flow methods, the proposed method boasts competitive accuracy-vs-cost trade-off, small memory requirements, ease of implementation and excellent parallelizability, providing a viable alternative for realistic two-phase flow calculations., Comment: Under review for publication in Journal of Computational Physics

Preferential concentration of inertial particles by turbulence is a well recognized phenomenon. This study investigates how this phenomenon impacts the mean heat transfer between the fluid phase and the particle phase. Using direct numerical simulations of homogeneous and isotropic turbulent flows coupled with Lagrangian point particle tracking, we explore this phenomenon over wide range of input parameters. Among the nine independent dimensionless numbers defining this problem, we show that particle Stokes number, defined based on large eddy time, and a new identified number called heat mixing parameter have the most significant effect on particle to gas heat transfer, while variation in other non-dimensional numbers can be ignored. An investigation of regimes with significant particle mass loading, suggests that the mean heat transfer from particles to gas is hardly affected by momentum two-way coupling. Using our numerical results we propose an algebraic reduced order model for heat transfer in particle-laden turbulence.

Superhydrophobic surfaces demonstrate promising potential for skin friction reduction in naval and hydrodynamic applications. Recent developments of superhydrophobic surfaces aiming for scalable applications use random distribution of roughness, such as spray coating and etched process. However, most of previous analyses of the interaction between flows and superhydrophobic surfaces studied periodic geometries that are economically feasible only in lab-scale experiments. We conduct direct numerical simulations of turbulent flows over randomly patterned interfaces considering a range of texture widths $w^+\approx 4-26$, and solid fractions $\phi_s=11\%$ to $25\%$. Slip and no-slip boundary conditions are implemented in a pattern, modeling the presence of gas-liquid interfaces and solid elements. Our results indicate that slip of randomly distributed textures under turbulent flows are about $30\%$ less than those of surfaces with aligned features of the same size. In the small texture size limit $w^+\approx 4$, the slip length of the randomly distributed textures in turbulent flows is well described by a previously introduced Stokes flow solution of randomly distributed shear-free holes. By comparing DNS results for patterned slip and no-slip boundary against the corresponding homogenized slip length boundary conditions, we show that turbulent flows over randomly distributed posts can be represented by an isotropic slip length in streamwise and spanwise direction. The average pressure fluctuation on gas pocket is similar to that of the aligned features with the same texture size and gas fraction, but the maximum interface deformation at the leading edge of the roughness element is about twice larger when the textures are randomly distributed., Comment: 15 pages, 7 figues

Concentration gradients in a fluid along a reactive surface due to contrast in surface reactivity generate convective flows. These flows result from contributions by electro and diffusio osmotic phenomena. In this study we have analyzed reactive patterns that release and consume protons, analogous to bimetallic catalytic conversion of peroxide. Here, we present a simple analytical model that accurately predicts the induced potentials and consequent velocities in such systems over a wide range of input parameters. Our model is tested against direct numerical solutions to the coupled Poisson, Nernst-Planck, and Navier-Stokes equations. Our analysis can be used to predict enhancement of mass transport and the resulting impact on overall catalytic conversion, and is also applicable to predicting the speed of catalytic nanomotors.

The clustering of small heavy inertial particles subjected to Stokes drag in turbulence is known to be minimal at small and large Stokes number and substantial at $\rm St = \mathcal O(1)$. This non-monotonic trend, which has been shown computationally and experimentally, is yet to be explained analytically. In this study, we obtain an analytical expression for the Lyapunov exponents that quantitatively predicts this trend. The sum of the exponents, which is the normalized rate of change of the signed-volume of a small cloud of particles, is correctly predicted to be negative and positive at small and large Stokes numbers, respectively, asymptoting to $\tau Q$ as $\tau \to 0$ and $\tau^{-1/2} |Q|^{1/4}$ as $\tau \to \infty$, where $\tau$ is the particle relaxation time and $Q(\tau)$ is the difference between the norm of the rotation- and strain-rate tensors computed along the particle trajectory. Additionally, the trajectory crossing is predicted only in hyperbolic flows where $Q<0$ for sufficiently inertial particles with a $\tau$ that scales with $|Q|^{-1/2}$. Following the onset of crossovers, a transition from clustering to dispersion is predicted correctly. We show these behaviors are not unique to three-dimensional isotropic turbulence and can be reproduced closely by a one-dimensional mono-harmonic flow, which appears as a fundamental canonical problem in the study of particle clustering. Analysis of this one-dimensional canonical flow shows that the rate of clustering, quantified as the product of the Lyapunov exponent and particle relaxation time, is bounded by $-1/2$, behaving with extreme nonlinearity in the hyperbolic flows and always remaining positive in the elliptic flows. These findings, which are stemmed from our analysis, are corroborated by the direct numerical simulations., Comment: full article, 16 figures

Drag laws for particles in fluids are often expressed in terms of the undisturbed fluid velocity, defined as the fluid velocity a particle sees before the disturbance develops in the fluid. In two-way coupled point-particle simulations the information from the undisturbed state is not available and must be approximated using the disturbed velocity field. Horwitz and Mani (2016) recently developed a procedure to estimate the undisturbed velocity for particles moving at low Reynolds number and obeying the Stokes drag law. Using our correction, we demonstrated convergence of numerical simulations to expected physical behavior for a range of canonical settings. In this paper we further extend and examine that correction scheme for particles moving at finite Reynolds number, by considering the Schiller-Nauman drag law. Tests of particle settling in an otherwise quiescent fluid show reasonable predictions of settling velocity history for particle Reynolds numbers up to ten. The correction scheme is shown to lead to $\geq 70\%$ reduction in particle velocity errors, compared to standard trilinear interpolation of disturbed velocities. We also discuss the modelling of unsteady effects and their relation to Stokes number and density ratio. Finally, we propose a regime diagram to guide scheme selection for point-particle modelling., Comment: Submitted to International Journal of Multiphase Flow. 14 pages, 7 figures

We present an efficient and robust numerical model for simulation of electrokinetic phenomena in porous networks over a wide range of applications including energy conversion, desalination, and lab-on-a-chip systems. Coupling between fluid flow and ion transport in these networks is governed by the Poisson-Nernst-Planck-Stokes equations. These equations describe a wide range of transport phenomena that can interact in complex and highly nonlinear ways in networks involving multiple pores with variable properties. Capturing these phenomena by direct simulation of the governing equations in multiple dimensions is prohibitively expensive. We present here a reduced order computational model that treats a network of many pores via solutions to 1D equations. Assuming that each pore in the network is long and thin, we derive a 1D model describing the transport in pore's longitudinal direction. We take into account the non-uniformity of potential and ion concentration profiles across the pore cross-section in the form of area-averaged coefficients in different flux terms representing fluid flow, electric current, and ion fluxes. Distinct advantages of the present framework include: a fully conservative discretization, fully bounded tabulated area-averaged coefficients without any singularity in the limit of infinitely thick electric double layers (EDLs), a flux discretization that exactly preserves equilibrium conditions, and extension to general network of pores with multiple intersections. By considering a hierarchy of canonical problems with increasing complexity, we demonstrate that the developed framework can capture a wide range of phenomena. Example demonstrations include, prediction of osmotic pressure built up in thin pores subject to concentration gradient, propagation of deionization shocks and induced recirculations for intersecting pores with varying properties.

Recent investigations have revealed that ion transport from aqueous electrolytes to ion-selective surfaces is subject to electroconvective instability that stems from coupling of hydrodynamics with electrostatic forces. Electroconvection is shown to enhance ion mixing and the net rate of transport. However, systems subject to electroconvection inherently involve fluid density variation set by salinity gradient in the bulk fluid. In this study we thoroughly examine the interplay of gravitational convection and chaotic electroconvection. Our results reveal that buoyant forces can significantly influence the transport rates, otherwise set by electroconvection, when the Rayleigh number $Ra$ of the system exceeds a value $Ra \sim 1000$. We show that buoyancy forces can significantly alter the flow patterns in these systems. When the buoyancy acts in the stabilizing direction, it limits the extent of penetration of electroconvection, but without eliminating it. When the buoyancy destabilizes the flow, it alters the electroconvective patterns by introducing upward and downward fingers of respectively light and heavy fluid.

We have developed a fully parallel C++/MPI based simulation code for variable-density particle-laden turbulent flows. The fluid is represented through a uniform Eulerian staggered grid, while particles are modeled using a Lagrangian point-particle framework. Spatial discretization is second-order accurate, and time integration has a fourth-order accuracy. Two-way coupling of the particles with the background flow is considered in both momentum and energy equations. The code is fully modular and abstracted, and easily can be extended or modified. We have considered two different boundary conditions. We have also developed a novel parallel linear solver for the variable density Poisson equation that arises in the calculation., Comment: In 2015, Annual Research Briefs, Center for Turbulence Research, Stanford University

An asymptotic solution is derived for the motion of inertial particles exposed to Stokes drag in an unsteady random flow. This solution provides the finite-time Lyapunov exponents as a function of Stokes number and Lagrangian strain- and rotation-rates autocovariances. The sum of these exponents, which corresponds to a concentration-weighted divergence of particle velocity field, is considered as a measure of clustering. For inertial particles dispersed in an isotropic turbulent flow our analysis predicts maximum clustering at an intermediate Stokes number and minimal clustering at small and large Stokes numbers. Direct numerical simulations are performed for quantitative validation of our analysis, showing a reasonable agreement between the two., Comment: 11 pages and 4 figures

In this work, we propose and test a method for calculating Stokes drag applicable to particle-laden fluid flows where two-way momentum coupling is important. In the point-particle formulation, particle dynamics are coupled to fluid dynamics via a source term that appears in the respective momentum equations. When the particle Reynolds number is small and the particle diameter is smaller than the fluid scales, it is common to approximate the momentum coupling source term as the Stokes drag. The Stokes drag force depends on the difference between the undisturbed fluid velocity evaluated at the particle location, and the particle velocity. However, owing to two-way coupling, the fluid velocity is modified in the neighborhood of a particle, relative to its undisturbed value. This causes the computed Stokes drag force to be underestimated in two-way coupled point-particle simulations. We develop estimates for the drag force error as function of the particle size relative to the grid size. We then develop a correction method that estimates the undisturbed fluid velocity from the computed disturbed velocity field. The correction scheme is tested for a particle settling in an otherwise quiescent fluid and is found to reduce the error in computed settling velocity by an order of magnitude compared with common interpolation schemes., Comment: Submitted to Journal of Computational Physics. 38 pages, 16 figures

Most electrochemical processes, such as electrodialysis, are limited by diffusion, but in porous media, surface conduction and electro-osmotic flow also contribute to ionic fluxes. In this paper, we report experimental evidence for surface-driven over-limiting current (faster than diffusion) and deionization shocks (propagating salt removal) in a porous medium. The apparatus consists of a silica glass frit (1 mm thick with 500 nm mean pore size) in an aqueous electrolyte (CuSO$_4$ or AgNO$_3$) passing ionic current from a reservoir to a cation-selective membrane (Nafion). The current-voltage relation of the whole system is consistent with a proposed theory based on the electro-osmotic flow mechanism over a broad range of reservoir salt concentrations (0.1 mM - 1.0 M), after accounting for (Cu) electrode polarization and pH-regulated silica charge. Above the limiting current, deionized water ($\approx 10 \mu$ $M$) can be continuously extracted from the frit, which implies the existence of a stable shock propagating against the flow, bordering a depleted region that extends more than 0.5mm across the outlet. The results suggest the feasibility of "shock electrodialysis" as a new approach to water desalination and other electrochemical separations., Comment: 39 pages, 9 figs