## Abstract

Understanding force generation in nonequilibrium systems is a notable challenge in statistical physics. We uncover a fluctuation-induced force between two plates immersed in homogeneous isotropic turbulence using direct numerical simulations. The force is a nonmonotonic function of plate separation. The mechanism of force generation reveals an intriguing analogy with fluctuation-induced forces: In a fluid, energy and vorticity are localized in regions of defined length scales. When varying the distance between the plates, we exclude energy structures modifying the overall pressure on the plates. At intermediate plate distances, the intense vorticity structures (worms) are forced to interact in close vicinity between the plates. This interaction affects the pressure in the slit and the force between the plates. The combination of these two effects causes a nonmonotonic attractive force with a complex Reynolds number dependence. Our study sheds light on how length scale–dependent distributions of energy and high-intensity vortex structures determine Casimir forces.

## INTRODUCTION

A fluctuating medium can exert a force on boundaries that confine the fluctuation (*1*). The most celebrated of these fluctuation-induced forces is the quantum Casimir effect: Metallic plates in a vacuum experience an attractive force because they confine vacuum fluctuations of the electromagnetic field. Recent studies have shown that the attractive force is also seen in systems containing active media that continuously consume energy on confining boundaries (*2*–*4*). Systems where large-scale structure and spatiotemporal chaos emerge through energy injection at small scales are referred to as active turbulence. The magnitude of the attractive Casimir-like force observed in these systems is determined by the partitioning of energy across various length scales. On the other hand, in classical hydrodynamic turbulence, energy is injected at large scales and allowed to transfer to smaller scales. Hydrodynamic turbulence is a paradigmatic nonequilibrium system with coherent spatial structures and well-studied mechanisms of energy transfer and dissipation. The phenomenon of force generation observed in other systems such as active media has not yet been studied in classical hydrodynamic turbulence. Understanding the interactions between objects in turbulence is crucial in explaining a plethora of phenomena, ranging from collective dynamics of plankton to volcano eruptions and multiphase flows in industrial processes.

In this work, we report the first observation of a force between objects immersed in a turbulent Newtonian fluid driven in a homogeneous and isotropic manner. We achieve this by performing direct numerical simulation of the governing Navier-Stokes (NS) equations where all the spatiotemporal scales are fully resolved without any background modeling. The nonlinearities in the NS equations cause the energy injected into the fluid cascade across length scales, energizing smaller and smaller structures. This process occurs up to length scales where viscosity dominates (the Kolmogorov scale), whereby energy is dissipated. If the forcing is continuous in time, then the system reaches a statistically stationary state with a continuous energy spectrum. Intermediate scales are energetic and play an active role in the cascade. They are known as “inertial” length scales and are neither artificially forced nor heavily dissipative. In this inertial range, there is a robust relationship between energy and wave number known as the “−5/3 law” (*5**)*, with the longest wavelengths containing the most energy.

In such a system, which is rich with nonlinearities, energy cascade, and spatiotemporal chaos, we show the existence of an attractive fluctuation-induced “Casimir” force and find that the force is a nonmonotonic function of plate separation. By analyzing the pressure distribution between the immersed plates as well as the location of the force maximum, we show that the mechanism of force generation is linked to the ability of the plates to pack, or exclude, specific flow features in between them and thus causing a Casimir-like force. Previous numerical studies of rigid spheres (*6*, *7*) and deformable bubbles (*8*) in turbulence focused on the wake instabilities of a single particle or the collective long-range effects of swarms of particles on the underlying turbulence. These results provide a mechanistic characterization of how hydrodynamic turbulence generates force between objects.

## RESULTS

Homogeneous isotropic turbulence (HIT) (*9*) is an idealized state of turbulence that can be approximated numerically by using a triply periodic computational domain, which is randomly forced at the largest wavelengths that fit in the computational box. We perform direct numerical simulations of HIT using the incompressible NS equations to study the nature of force between two closely placed objects immersed in a forced fluid. The NS equations relate fluid velocity **u** and pressure *p**t* is time, ν is the fluid kinematic viscosity, and ρ_{F} is the fluid density. **F** = **f**^{tur} + **f**^{ibm}is a force vector that is a sum of two terms; i.e., **f**^{tur} is a contribution from the forcing needed to generate homogeneous turbulence in the domain, while **f**^{ibm} is the force needed to enforce the influence of the immersed plates on the flow through the immersed boundary method, which will be described later.

We use a computational periodic cube of side ℒ, with two square parallel plates of size of *l _{p}*/ℒ = 0.25, and vary the plate distance

*d*between

*d*/ℒ = 0.05 and 0.25. Taylor-Reynolds numbers of

*Re*

_{λ}=

*u*′λ/ν = 65, 100, and 140 are simulated, where λ is the Taylor microscale

*Re*

_{λ}is a measure for the range of scales covered by the inertial range. The values simulated are enough to provide a well-developed inertial range while limiting computational costs.

Figure 1 is the central result of this work: Turbulent flow generates an attractive force between rigid plates, which increases with *Re*_{λ} and is a nonmonotonic function of plate separation. We quantify the force with the dimensionless coefficient *F* is the temporally averaged normal force on the plates and *A* is the plate area. To ensure that the force is not a numerical artifact caused by the sharp edges of the plates, we also show the results of a simulation with two spheres of radius *r*/ℒ = 0.05. In this case, the attractive force is weaker but still present and persistent over long time scales.

## DISCUSSION

We now discuss the nature and origin of the computed attractive force. Any surface (plate here) immersed in a flow field experiences two kinds of stresses: viscous stresses and pressure, which act in the tangential and normal direction, respectively. A difference in local pressure on both sides of a rigid plate leads to a force along the normal direction. The pressure in a fluid is directly related to the velocity gradients in the flow, i.e., the curl of velocity, also known as vorticity (ω), and the rate of strain in the fluid σ by the relation ∇^{2}(*p*/ρ_{F}) = ω^{2}/2 − σ^{2} = 2*Q* (*5*), where *Q* acts like a source in the pressure evolution equation. A series of classical studies in the past have looked at the consequences of the above equation and the nature of distribution of pressure in the flow (*10*–*12*). In a homogeneous and isotropic turbulent flow field, large-scale structures and strain-dominated regions are markers for positive pressure fluctuations, while thin and tubular high-intensity vorticity structures lead to negative pressure fluctuations. Another interesting observation is that the pressure distribution is highly skewed toward negative pressure with exponential tails as opposed to more Gaussian-like behavior toward positive pressure fluctuations. This is primarily a consequence of the interaction between strain-dominated and vorticity-dominated regions in the flow where it is clearly observed that high-strain regions are typically associated with high vorticity, while the converse does not necessarily hold (*10*). This leads to a bias in the source term (*Q*) for the pressure Poisson equation and consequently extremely low-pressure fluctuations within regions of high-intensity vorticity structure. The interaction between structures with different characteristics and length scales leads to a nontrivial pressure distribution in the flow.

The forces acting on the plates are directly governed by the local pressure, which, in turn, is governed by the organization of the underlying energy field and vorticity. While energy is contained in the largest wavelengths, regions of intense vorticity organize themselves into tubular-like structures, referred to as worms in literature (*13*). The radius of these tubular structures scales as the Kolmogorov length η_{K} = (ν^{3}/ε)^{1/4} (a measure of the viscous cutoff length), where ε is the energy dissipation in the system. Their length scales as the integral, or decorrelation, length scale. An estimate for this decorrelation length scale is provided by the large-eddy length scale, defined as *L* = *k*^{3/2}/ε, with *k* = 3/2ρ_{F}*u*′^{2} the kinetic energy of the flow and *u*′ the root mean squared velocity fluctuation in one direction (in our simulations, *L* is approximately constant). The complex organization of energy and vorticity into spatial structures in hydrodynamic turbulence, which affects the nature of positive and negative pressure fluctuations combined with the geometrical confinement induced by the plates, leads to the Casimir-like force.

A single plate in a randomly forced fluid will feel, on average, symmetric forces on both sides. However, two rigid plates placed next to each other restrict the structures that “fit” between the plates. Such a geometrical confinement completely excludes high-energy structures except when the plates are far apart. The geometrical confinement also affects the magnitude and direction of intense vorticity structures that fit between the plates (see the Supplementary Materials for a quantification of anisotropy). Both effects modify the pressure fluctuations and generate a net force. The mechanism of force generation is analogous to thermal or quantum Casimir forces (*1*), where a restriction in the fluctuations imposed by the boundaries leads to force generation. Here, as the energy and vorticity fluctuations organize themselves into defined and distinct length scales, we would expect the Casimir force to be nonmonotonic, with the peak force achieved when the plate separation is comparable to the length scale of the excluded structures (*4*). As *Re*_{λ} increases, the magnitude of the force increases, possibly because of increase in the intensity of thin vortex structures, which are markers for negative pressures. The plate distance at which optimal force is achieved varies between *Re*_{λ} = 65 and 140. To fully understand how the plate distance for optimal force generation varies with Reynolds number, a denser sweep of the parameter space is required, and this is anticipated in future studies.

To quantitatively understand the above arguments, we examine how the fluctuation modes distribute themselves inside the slit. Figure 2 shows an instantaneous snapshot of the flow for two plate distances. The left panel clearly shows that the energetic structures have lengths comparable to or larger than the plates. They can only penetrate between the plates for the large plate distances. Meanwhile, the right panels clearly show the tubular vorticity structures. By comparing both figures, we can see how for close distances, the plates organize the tubular structures within, while for large distances, they are more disordered.

The time-averaged kinetic energy and the marker for intense vortex structures *Q* in the planes normal to the plates (Fig. 2, C and D) further corroborate the relationship between energy, vorticity structures, and force generation. A narrower plate separation leads to lower mean energy inside the walls, agreeing with our physical picture that the plates filter the fluctuations in between and outside the plates. However, no qualitative change is seen in the profiles of kinetic energy close to the maximum, and the nonmonotonicity of the force cannot be explained from energy statistics alone. For the *Q* profiles, there is a qualitative change of behavior around the plate separation at which the maximum force is attained (green and yellow curves). The mean value of *Q* increases within the gap for plate separations at which the maximum force is found, indicating the significance of structures where vorticity dominates strain rate (*14*). At very close distances between the plates, viscous effects damps the evolution of the vortical structures in the gap. The whole flow becomes quiescent between the plates, and no structures fit. As the plate separation increases, the region between the plates becomes populated, first by intense vorticity structures and later by energy-containing structures.

By immersing the plates into the flow, we selectively filter the flow fields. In Fig. 3, we plot the pressure within and outside the immersed plates. We note that while the sharp edges of the plates are clearly visible in this figure, this is not the source of the force. To ensure this, in Fig. 1, we reported that an attractive force was also found between two spheres. Instead, we can understand the force from the fact that when we place two plates in the flow, we exclude the energy-containing structures from between the plates, which generates high pressures at the plate centers. At the same time, we force the intense vorticity structures to interact in close quarters within the confinement, which enhances vortex stretching and produces even more intense vortex structures, which decreases the pressure (*13*). The pressure on the outside section of the plate can be seen to increase greatly between the left and middle pictures, as structures are excluded from a larger volume. The pressure in the inside remains relatively low in comparison to the outside, generating the attractive force. As the plates separate further, the prominent peak in the pressure fluctuations at the midgap fades, and the pressure fluctuations profile between the plates becomes horizontal; the interaction between intense vortical structures is now less localized. This coincides with a sharp decrease in the force magnitude between the plates.

Because of the intimate link between vorticity and pressure, one can expect this effect to show up in the pressure fluctuations, shown in Fig. 3B. Around *d*/ℒ = 0.075, a sharp peak of pressure fluctuations develops at the midgap. The peak is present for the plate separations that cause the highest attractive forces. The profiles of *Q* (in Fig. 2B) and the pressure fluctuations indicate the enhancement of intense vortex structures (source of low-pressure regions), which also contributes to optimal force generation. To reassure the reader that the force in not a numerical artifact, we further characterize the nature of force generation in the Supplementary Materials.

In summary, we demonstrate the existence of an attractive turbulent Casimir-like force by direct numerical simulation of two plates immersed in a homogeneous isotropic turbulent flow. Our study sheds light on the interaction between objects in a turbulent flow at very close distances. We anticipate that this force can be experimentally realized in setups previously used for studying single-phase and multiphase HIT (*15*). Further research to understand the *R*e_{λ} dependence is needed. Moreover, we speculate that this mechanism of generating a nonmonotonic Casimir-like force due to nontrivial spatial redistribution of energy containing structures and intense vortex filaments might be a general phenomenology for active and nonequilibrium systems and could appear in large classes of biological fluids comprising microbial suspensions that exhibit notable analogies with turbulent flows (*16*–*19*).

## MATERIALS AND METHODS

The NS equations are solved using an energy-conserving second-order centered finite difference scheme in a Cartesian domain with fractional time stepping. An explicit Adams-Bashforth scheme is used to discretize the nonlinear terms, while an implicit Crank-Nicholson scheme is used for the viscous terms (*20*). Time integration is performed via a self-starting fractional-step third-order Runge-Kutta scheme, and the time step is dynamically chosen so that the maximum Courant-Friedrich-Lewy condition number is 1.2. The domain is periodic in all three directions with a periodicity length *L*. To spatially discretize the equations, a cubic grid is used, with 240^{3} points for *Re*_{λ} = 65, 360^{3} points for *Re*_{λ} = 100, and 480^{3} points for *Re*_{λ} = 140. The formulation of the force vector **f**^{tur} is based on random processes driving the time evolution of a selected number of large scales (or small–wave number modes). Additional details on the forcing scheme and its corresponding parameters can be found in the study by Eswaran and Pope (*21*) and the study by Chouippe and Uhlmann (*22*). Further details on how this forcing injects energy and how the forces on the objects are affected by it can be found in the Supplementary Materials.

The influence of rigid plates on the surrounding fluid is simulated using an immersed boundary method based on the moving least squares approximation. This involves discretizing the immersed plates using several (∼10^{4}) triangular computational elements. Further details about the method are described in the Supplementary Materials (*23*, *24*). The forces acting on the immersed rigid plates can be computed from the pressure interpolated on both sides individual triangular elements. Temporal convergence of the forces was assured by running the simulations until the force from the hydrodynamic pressure on both plates were equal (but oppositely signed) to within 3%. This defined the magnitude of the error bars. In practice, this meant running up to more than 100 large-eddy turnover times defined by *T _{e}* =

*u*′

^{2}/ε.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/14/eaba0461/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank R. Verzicco, A. Pumir, Y. B. Sinai, and M. P. Brenner for valuable discussions.

**Funding:**This work was completed, in part, with resources provided by the University of Houston Center for Advanced Computing and Data Science. A.A.L. acknowledges the support of the Winton Programme for the Physics of Sustainability.

**Author contributions:**Conceptualization, A.A.L.; methodology, R.O.-M.; software and formal analysis, V.S.; running simulations, R.O.-M., V.S., and D.P; writing, all authors.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).