Table of Contents
Journal of Fluids
Volume 2016 (2016), Article ID 3587974, 15 pages
Research Article

A Volumetric Approach to Wake Reduction: Design, Optimization, and Experimental Verification

1Duke University, P.O. Box 90300, Hudson Hall, Durham, NC 27705, USA
2Duke University, P.O. Box 90291, Hudson Hall, Durham, NC 27708, USA
3Office of Naval Research, Naval Undersea Warfare Center, Newport, RI 02841, USA

Received 22 October 2015; Revised 23 February 2016; Accepted 20 March 2016

Academic Editor: Jose M. Montanero

Copyright © 2016 Dean R. Culver et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Wake reduction is a crucial link in the chain leading to undetectable watercraft. Here, we explore a volumetric approach to controlling the wake in a stationary flow past cylindrical and spherical objects. In this approach, these objects are coupled with rigid, fluid-permeable structures prescribed by a macroscopic design approach where all solid boundaries are parameterized and modeled explicitly. Local, gradient-based optimization is employed which permits topological changes in the manifold describing the composite solid component(s) while still allowing the use of adjoint optimization methods. This formalism works below small Reynolds number (Re) turbulent flow (–10,000) when simulated using small Reynolds-averaged Navier-Stokes (RANS) models. The output of this topology optimization yields geometries that can be fabricated immediately using fused deposition modeling (FDM). Our prototypes have been verified in an experimental water tunnel facility, where the use of Particle Image Velocimetry (PIV) described the velocity profile. Comparisons with our computational models show excellent agreement for the spherical shapes and reasonable match for cylindrical shapes, with well-understood sources of error. Two important figures of merit are considered: domain-wide wake (DWW) and maximum local wake (MLW), metrics of the flow field disturbance whose definitions are described.

1. Introduction

Hydrodynamic drag and the closely related phenomenon of wake are two of the most important characteristics of a vessel. While drag is directly responsible for the bulk of the energy expenditure of any propulsion system, wake properties describe the disturbance left behind in the flow field and, therefore, particularly in the cases of surface-penetrating vehicles, how easy it is to detect the vessel itself. Much of the prior research in vessel hydrodynamics has been focused on the control and reduction of drag. The methods for achieving this are numerous, including active and passive structural and passive chemical approaches. Existing work in hydrodynamic flow control via structures includes tangential surface motion control on a cylinder [1] and active and passive structures in sea creatures [2]. Chemistry and fluid heterogeneity offer another avenue into drag reduction, such as studies considering hydrophobic surfaces [3], surfactant injection [4], bubble injection [5], and supercavitation [6]. Attention to the wake in hydrodynamics was given mostly to the extent that it is useful to control the drag coefficient. In aerodynamics, wake control has found substantial relevance in aircraft maneuverability [7], efficiency of morphing wings [8], and noise mitigation [9]. Another active area of research is passive and active turbulence control; efforts were made to leverage the critical Reynolds number and postpone the onset of turbulence or flow separation through, for example, synthetic jets [10] and oscillating walls [11]. Active control of fully turbulent flows in order to lower the drag is also currently investigated; see, for example, [12] for further references on that topic.

While these approaches to drag reduction are somewhat successful, they focus almost exclusively on the boundary, attempting to control or modify the flow through the modifications of the boundary shape, or active interactions with the boundary layer, as in the case of synthetic jets and artificial cilia [13]. The possibility of controlling the flow with a multilayered, three-dimensional, volumetric structure has not received much attention until recently, due to the difficulty of modeling, designing, and fabricating any complex, three-dimensional structures. This is especially true for active scenarios. Urzhumov and Smith [14, 15] have investigated the possibilities in volumetric flow control using a simplified, two-scale description based on the local permeability of a porous medium. Motivated by recent advances in transformation electromagnetics [16] and acoustics [17], they considered both passive and active scenarios for the control of the Stokes flow around a spherical object and found that complete cancellation of both the drag term and the domain-wide wake requires an active system. These studies [14, 15] did not address the question as to whether or not any wake reduction can be achieved in hydrodynamics by passive means and also left out any consideration of flow nonlinearity and turbulence. With ever-improving computational hardware and numerical methods, it has recently become possible to perform complete physical modeling of a highly complex fluid-solid composite without the need for local homogenization. Moreover, with solution times in the range of seconds, it is feasible to perform optimization of the shapes and topologies of the solid phase, provided that smart choices in the shape parameterization and optimization algorithms are made. While global optimization of fully three-dimensional fluid-solid composites is still beyond reach for a single-CPU computer, two-dimensional and axisymmetric stationary flows through passive structures are now amenable to exhaustive investigations. The task remains feasible even with inclusion of developed-turbulence models based on RANS equations, such as, for example, the Spalart-Allmaras model [18] or a low-Re - model, both available in a commercial Computational Fluid Dynamics (CFD) simulator, COMSOL Multiphysics. This paper reports our initial effort in that direction, supplemented by experimental results on a rapid prototype sample measured in a particle velocimetry tunnel.

2. Materials and Methods

2.1. Wake and Drag Manipulation with Homogenizable Porous Media

Porous media saturated with a single fluid have been the subject of many theoretical and experimental studies since the works of Darcy [20] and Brinkman [21]. In the creeping flow limit, homogenization of porous media has been rigorously achieved by applying asymptotic theories to the Navier-Stokes model [22]. However, Brinkman-Stokes equation [23] is generally believed to be applicable even at nonvanishing Reynolds numbers (), with or without the nonlinear Forchheimer correction [24]. And yet, homogenization of turbulent porous flows remains a highly complicated subject [24, 25]. The Brinkman equation readswhere is the true (microscopic) fluid viscosity. Here, we assume that the porosity is close to 100%, implying that a small solid filling fraction in the porous media is considered. The last term on the right-hand side of (1) is the Brinkman term corresponding to the linear (Darcy) filtration law; it depends on the effective medium parameter , which is the permeability.

Permeability and porosity are the constitutive parameters of composite media that can be controlled by engineering the microstructure of the medium. Recently, composite media with engineered effective medium properties have become known as metamaterials. The main distinction between naturally occurring composite media and metamaterials is the precisely controlled and often periodic microstructure of the latter.

The effect of radial permeability distributions on the flow past spherical structures has been a subject of several studies [14, 26, 27]. Recently, the effects of these distributions on the turbulent wake of spherically symmetric structures [28] and the transition Reynolds number of cylindrical structures [15] were considered. However, the effect of more spatially complex permeability distributions remains unknown, due to both the computational cost of full three-dimensional solutions and the much larger functional space of permeability maps to be explored. Performing such a study is our starting point for evaluating the potential of passive permeable structures for wake reduction.

Here, we consider an impermeable spherical body enclosed in a permeable, spherical shell (the envelope). Ideally, we wish to know the effect of every physically possible permeability distribution on the main properties of the stationary flow, such as wake and drag. In order to sample the infinitely dimensional functional space of all possible permeability distributions, we consider the following parametrization, which introduces a finite, selectable number of parameters. One set () of these parameters describes the radial variation of permeability, and two other sets ( and ) describe the angular variation. Specifically,create an 8-dimensional parametric space for any in the layer, that is, forwhere is the number of layers. In the above, we keep only the first four coefficients on the Fourier series in the angle. The generalization to an arbitrary number of angular degrees of freedom is intuitive. In our numerical calculations, we further restrict ourselves to four radial layers. We then sample this eight-dimensional parameter space by creating a uniform grid within its boundaries and solving the Brinkman-Stokes equation (1) for each spatial distribution. The specific (exponential) form in (2) was chosen to ensure that the permeability distributions remain positive everywhere.

The Monte Carlo sweep described above was carried out for a structure having an impermeable core and an aspect ratio of 5 : 1 (Figure 1). The manifold describing the relationship between the wake and the drag of nonspherically symmetric structures indicates that passive permeable media with controlled permeability distributions can offer significant trade-offs between the drag and the wake of the structure. The manifold describing the full range of achievable wake and drag combinations is illustrated in Figure 2.

Figure 1: (a) The system schematic illustrates the principal cross section of a spherical body with radius enclosed in an envelope of prescribed permeability of radius . (b) Permeability distribution related to one of the lowest wake values in the wake-drag manifold at the maximum drag of a solid sphere; envelope aspect ratio ( : ) is 5 : 1. (c) Flow field resulting from the permeability distribution in (b).
Figure 2: Monte Carlo sampling of the 8-dimensional functional space of permeability distributions (described by the coefficients in (2)) and its projection onto the wake-drag plane. Evaluation points are collected in scatter plots illustrating the boundaries of the wake-drag manifold for (a) Reynolds number (Re) = 20 and (b) Re = 100. For the wake norm, the dimensionless domain-wide wake (DWW) figure is used; see (8) for the definition of DWW. The drag force indicated along the abscissa is normalized by the Stokes drag of a solid sphere of radius . The shaded rectangle on each plot illustrates the maximum wake and drag possible from a solid sphere of radius . The solid black curve shows DWW and drag for various solid spheres of smaller radii. Note that the manifold quickly grows in vertical and horizontal dimensions with an increasing Reynolds number.

Figure 2 offers several pieces of insight into the nature of the wake-drag manifold for various permeability distributions. First, the solid (impermeable) sphere of radius provides neither the maximum drag nor the maximum wake; it is therefore possible to have a permeable structure with substantially more drag and/or wake than what is expected of a solid object of the same outer diameter. This finding could have interesting applications in creating highly inertial buoys whose effective drag coefficient is unnaturally high. Second, there are structures whose drag coefficients are somewhat smaller than the drag of a solid sphere of radius ; these can be observed as dots on the top bar of the shaded rectangle, near its upper-right corner. Third, and most important for the purpose of this study, there is a wide distribution of points along the vertical right side of the shaded rectangle. These points correspond to permeable structures whose wake norm is less than the wake norm of a solid sphere with the exact same drag coefficient. In other words, these structures achieve wake reduction without an increase in the drag. This comprises theoretical proof that wake reduction can be achieved with permeable structures both completely passively and without an additional energy expenditure. In short, the wake of a given body may be reduced by introducing optimally designed permeable structures on its surface without changing drag.

Finally, for larger Reynolds numbers, the size of this manifold increases rapidly, implying that the efficiency of wake and drag mitigation using permeable media is substantially better for larger bodies or at higher speeds. As Re equals 100, wake reduction factor relative to a sphere with equivalent drag already exceeds a factor of two up from a factor of 25% at , indicating the potential for even larger reduction at .

2.2. Computational Methods and Design Methodology

While simple homogenization theories that yield the Brinkman equation are considered reliable at small Reynolds numbers, their usefulness in the turbulent regime remains controversial [25]. Inspired by the theoretical results of Section 2, we then pursue a design approach that avoids using homogenization in the domain of the permeable media. The entire structure, with the geometrically and topologically intricate shape of the solid and fluid phases, is modeled here explicitly using the full Navier-Stokes equations, with the standard no-slip boundary condition on the fluid-solid interfaces. The results of the grid, boundary, and domain size sensitivity studies can be found in Appendix A.

Since we are interested in evaluating a potentially very large class of complicated geometric structures, we develop a geometry parametrization scheme along the lines of the level-set method popular in two-phase flow modeling.

In this section, the Brinkman term is used again, but this time only to model the impermeable, solid phase of the fluid-solid composite in a way that is compatible with a level-set parameterization of the solid boundaries. Such compatibility enables subsequent use of efficient optimization techniques. Performing optimization on a level-set parameterized binary composite has been known as topology optimization in the design of solid composites since the works of Bendsoe and Kikuchi [29] and Bendsoe and Sigmund [30]; the technique was applied in Stokes hydrodynamics with the goal of minimized energy dissipation by Borrvall and Petersson [31] in a different context. When the permeability in the Brinkman term is small, the Brinkman term becomes large, effectively expelling the fluid out of the small-permeability region. On the other hand, where permeability is very large, the Brinkman term is negligible, and fluid flow is free.

Specifically, our parametrization of the inverse permeability iswhere is the level-set function assuming values between −1 and 1 and is a smoothed Heaviside function with a smoothing width . With this parametrization, the regions with represent the solids and regions with are solid-free. The fluid-solid boundary, given by , ends up having an effective no-slip boundary condition, due to the continuity of velocity between the solid regions (where velocity is essentially zero) and the free regions. Any continuous, bounded level-set function thus defines a particular composite structure; of interest to us are only structures that allow flow through themselves, which we term permeable volumetric composites (PVCs). PVC is a generalization of gradient-index metamaterials considered by Smith et al. [32]. Unlike the latter, PVCs are not required to have even approximate periodicity in space. These composites are wrapped around an impermeable object; the shape of this object can be arbitrary but is considered cylindrical or spherical in the examples below. The ultimate goal of this effort is to manipulate or, more specifically, reduce the wake of the impervious object to a degree that would be impossible to achieve using only the boundary-layer-control techniques cited in the introduction.

Below, we consider both spherical (as in Section 2) and cylindrical bodies, modeled using two-dimensional axisymmetric and two-dimensional CFD solvers, as the simulations will be operating at a Reynolds number before the onset of 3D effects [33]. Correspondingly, the domain occupied by the PVC was chosen to be either a cylindrical annulus (the two-dimensional case) or a spherical shell with an outside radius (as in Figure 1). In all cases, the radius of the envelope is chosen to be twice the impermeable core radius.

To implement this approach, we use a commercial software package, COMSOL Multiphysics version 4.4 [19]. The forward solutions are obtained using stationary Navier-Stokes equations in the laminar flow regime or using RANS models such as Spalart-Allmaras [18] in the turbulent flow regime, whose parameters are listed in Table 1. Several other turbulence models in the COMSOL package were tested, including -, -, and Shear-Stress Transport (SST), and the solutions were not strongly sensitive to the chosen method. The validity of the S-A model to bluff bodies and nonstreamlined structures is verified in several works, such as Schmidt and Thiele for flow over wall-mounted cubes [34]. Permeability as a coordinate-dependent parameter is readily available in the Subsurface Flow module of COMSOL Multiphysics, at least for laminar flow regimes. In turbulent (RANS) solvers, this parameter may be unavailable. In those cases, we use an equivalent description where the Brinkman term is introduced as an external volumetric force, which is then made dependent on the local flow velocity. Using the definition of inverse permeability from (4), this volume force may be appropriately prescribed with a given magnitude in the direction opposite the local flow as in

Table 1: Unitless coefficients used in the incompressible Spalart-Allmaras turbulence model, whose descriptions are available in the COMSOL CFD User’s Manual [19].

With this relationship, (1) becomes

This explicit form of the volume force can be suitably specified in CFD solvers that do not explicitly support porous effective media. With this approach, we were able to achieve consistent convergence while using the Spalart-Allmaras turbulence model [18] and the corresponding solver in COMSOL.

To reduce computational time, we consider two reduced-complexity scenarios. In the cylindrical or 2D case, the PVC-coated object is a cylinder, and the PVC shapes are obtained by extrusion of a 2D cross section out of its plane. In other words, the level-set function in (4) is independent of the -coordinate. In the spherical or axisymmetric case, the coated object is a sphere, and the PVC structure is a shape of revolution. In the approach, the PVC is either an extruded or revolved set of geometries with a small number of minimally invasive connectors added. The connectors are there to simply hold the samples together using the rapid prototyping techniques.

We selected the physical size of the impermeable cores and the free-stream flow parameters based on the limitations of the water tunnel used in the experiments described below. The choices for these parameters were based on ensuring that the finite-channel (wall proximity) and shallow-water effects of the tunnel were minimal, while maximizing the size of the PVC region to keep the minimum feature size of the structure as large as possible, so that fused deposition modeling (FDM) prototyping would be seamless. Our parameter selections also ensured operation within the specifications of the tunnel and avoided simulation parameters that could lead to computational nonconvergence. These constraints led us to use a free-stream velocity of 0.25 meters per second (m/s), a spherical (impervious) object radius of 27.5 millimeters (mm), and a cylindrical object with an 8 mm radius. The outer radii for the PVCs were twice the impervious core radii, with a spherical obstruction envelope radius of 55 mm and a cylindrical obstruction envelope radius of 16 mm. These parameters correspond to Reynolds numbers of 8000 (for the cylinders) and 27500 (for the spheres). Although different Reynolds numbers were tested (ranging from creep flow to turbulent flow), we settled on these Reynolds numbers for the experiment in order to meet the constraints of the resources available to us, such as the test chamber size, water tunnel velocity, and 3D printer supplier capabilities. These moderately high Re numbers well beyond the stable laminar flow regimes for cylinders and spheres necessitated the use of low-Re turbulence models in the forward solver.

The fluid (water) was approximated to have a density of 1000 kilograms per cubic meter () and a dynamic viscosity of pascal-seconds (). To perform optimization on a level-set parameterized PVC, we used local, gradient-based optimization techniques, such as those included in COMSOL Multiphysics. Similar optimization convergence rates and comparable results were achieved by using SNOPT and MMA optimization algorithms. Both of these algorithms are capable of using the highly efficient adjoint method for calculating the full sensitivity vector, which is needed for gradient-based optimization. The advantage of gradient-based techniques is their second order of convergence and therefore faster arrival at local minima in comparison with gradient-free techniques. On the other hand, these techniques are inherently local: they need to start with a prescribed initial guess, and they typically converge to a local minimum nearby. This is in contrast with quasiglobal techniques (genetic algorithm, Monte Carlo sampling, simulated annealing, etc.), which allow wider deviations from the initial conditions and sample a broader domain in parameter space. However, with the number of optimization degrees of freedom in the range of many thousands and forward solution times of order 1–10 seconds, it is not feasible to perform such elaborate optimization searches; we are thus focusing on gradient-based methods here.

For the initial guesses, we use a class of level-set functions constructed from periodic (such as trigonometric) functions as follows:where and are predetermined spatial frequencies in the radial and zenith angle directions, respectively. is the core cross section radius, is the PVC envelope radius, and is another parameter that controls the initial volume and size of the solid phase inclusions. As the optimization solver seeks minima for the specified wake metric, this original geometry morphs as advantageously permeable or impermeable regions are chosen. The final PVC volumes for the spherical and cylindrical cases are shown in Figure 3.

Figure 3: Three-dimensional CAD models of one-half of the PVC-coated cylinder (a) and the whole PVC-coated sphere with a mounting strut (b).

To achieve the results illustrated in Figure 3, the optimization goal chosen was the domain-wide wake norm, defined as a volume integral:where is the free-stream velocity, is the integration domain volume, is the velocity field, and is the unit vector in the direction of the free-stream velocity. We also studied the maximum of the local wake field (MLW), where the wake norm is defined asand the maximum local wake metric is defined as

Other wake metrics are possible, such as the minimum local pressure, domain-wide pressure deviation norm, or the volume of the region where pressure deviation is negative; these metrics are more relevant to the cavitation effects; however, they were not used in the optimization of structures reported here.

In performing the optimization of the PVC, we found that the downstream hemisphere (or semiannulus) is more important in the wake formation at nonnegligible Reynolds numbers, whereas in the Stokes limit the equatorial plane is essentially a symmetry plane of the flow, and the up- and downstream domains are equally important. This justifies the choices of initial guesses, which, after optimization, led to the shapes in Figure 3.

2.3. Fabrication and Experimental Methods

Four structures were examined over the course of the experiment, each offering less than a 10% channel area obstruction to minimize the wall effects on the near-field flow. All four structures, designed to be contiguous rigid bodies, were fabricated from ABS, a thermoplastic commonly used in fused deposition modeling 3D printers. The first two bodies were control samples: an impervious sphere and a cylinder spanning the height of the test section. The sphere was offset to be centered in the channel by a threaded rod. Both this threaded rod and a machine screw at the bottom of the cylinder interfaced with a mounting fixture in the bottom of the channel. The second two geometries were 3D-printed permeable volumetric composites, one of the cylindrical type and one of the spherical type as described in Section 2. Their shape and structure (again, see Figure 3) resulted from numerical optimization as described above, and these three-dimensional bodies are simply an extrusion of the cylindrical result and a revolution of the spherical result.

The prototype samples were generated from two-dimensional plots of their cross sections. In order to realize these composites in three dimensions, the cross section plots were first vectorized and converted to the DXF format and then imported into 3D CAD software such as Solidworks. Vectorization was accomplished using the trace bitmap function in the Inkscape freeware and by removing extraneous elements of the two-dimensional plot output manually. Once imported into three-dimensional Computer Aided Design (CAD), the cross sections were either extruded (for cylindrical geometries) or revolved around the symmetry axis (for spherical samples). Finally, connection elements were added manually in the CAD environment to ensure the structural integrity of the samples. Although these viscoelastic spacers alter the flow, they are designed to be as distant from the Particle Image Velocimetry (PIV) observation plane as possible. For example, the spacer cylinder connection point is a full quarter of the sample height away from the plane examined by the PIV method, and the spherical sample’s standoffs were not in-plane with the PIV sheet at all. These completed rigid solids were then three-dimensionally printed in Stratasys uPrint SE Plus with an ABS thermoplastic material.

Measurements of the flow velocity profiles were performed in the water tunnel housed at the Naval Undersea Warfare Center in Newport, Rhode Island. The water tunnel (Figure 4) is a two-story tall circuit capable of running at velocities in the 0.1–10 m/s range. The 305 × 305 mm test section is pictured in Figure 4 along with most of the data acquisition set-up for this experiment.

Figure 4: Photograph of the experimental set-up, showing a water PIV tunnel with a sample, the optical camera, and the scanning laser emitter.

The prototype samples were mounted in the same manner as the control samples. In the cylindrical cases, in order to avoid oscillation at the samples top end due to vortex shedding, a viscoelastic spacer was placed between the top of the samples and the ceiling of the channel. Furthermore, because the chamber is 305 mm in height and the three-dimensional printer was only capable of generating shapes that could fit inside of a 154 mm sided cube, two PVC-coated cylinders were assembled end-to-end to span this height.

Two-dimensional flow field behavior can be extracted from such a system via Particle Image Velocimetry (PIV). The PIV method involves projecting a laser sheet into a transparent flow field populated by reflective seeds, which are particles in the flow that will reflect the laser light so that it can be detected by a camera oriented perpendicular to the laser sheet. In Figure 5, the top view illustrates how the PIV lens sees a sheet of illuminated seed particles.

Figure 5: Schematic of our PIV data acquisition method.

When gathering data, the camera will capture two images in quick succession, and the PIV software can identify how certain seed elements in each image have moved. By calibrating the system to recognize how many pixels in the image constitute a known measurement and knowing the time between capturing the two images, local velocities may be identified. For the purposes of this experiment, thirty-five image pairs were averaged to illustrate the effective mean flow field.

This experiment was conducted with Litron Nano L-200-15 class 4 lasers and a Litron LR1550 generator. The laser was redirected and focused on the samples with an LAVision 1108453 mirror stabilizer arm, 1108405 lens, a 1101342 PIV camera with a Nikon AF 50 mm lens attached, and a cylindrical beam splitter. The beam splitter was mounted on top of a Slik AMT Tripod, and the data was processed using DaVis 7.2 PIV data acquisition software, which sampled the recording with a rectangular grid size of ninety-three points along a side. In order to appropriately compare, the computational data was interpolated to these same grid points.

3. Results and Discussion

Experimental and computational results for select figures of merit for both control and prototype cases of the sphere and cylinder are summarized in Table 2, including the difference between theory and experiment as well as percent effect of the PVC elements, which are called “cloaked” samples in this table.

Table 2: Results summary and comparison. Note that the physical location of the PVC cylinder cap caused masking errors in the experimental data, resulting in anomalous or unnecessarily absent readings near the geometry.

Figure 6 illustrates the computational behavior of the local wake near the uncloaked spherical subjects.

Figure 6: (a) Experimental local wake near the uncloaked sphere with Re = 27500. The range where reaches 8 cm. The region where spans a nearly elliptical region with a major axis of 3 cm and a minor axis of nearly 2 cm. (b) Computational local wake near the uncloaked sphere with Re = 27500. The range where reaches 11 cm. The region where spans a nearly elliptical region with a major axis of 4 cm and a minor axis of 2 cm.

Figure 7 illustrates the results concerning the behavior of the local wake near the cloaked spherical samples.

Figure 7: (a) Experimental local wake near the cloaked sphere with Re = 27500. The range where reaches 5 cm. The region where is eliminated from this test. (b) Computational local wake near the cloaked sphere with Re = 27500. The range where reaches 10 cm. The region where spans a nearly elliptical region with a major axis of 2 cm and a minor axis of less than 1 cm.

Figure 8 illustrates the behavior of the local wake near the uncloaked cylindrical subjects.

Figure 8: (a) Experimental local wake near the uncloaked cylinder with Re = 8000. The range where reaches 3 cm. A region where does not appear and (b) computational local wake near the uncloaked cylinder with Re = 8000. The range where reaches 5 cm. The region where spans a nearly elliptical region with a major axis of 2 cm and a minor axis of less than 1 cm.

Figure 9 illustrates the results concerning the behavior of the local wake near the cloaked cylindrical samples.

Figure 9: (a) Experimental local wake near the cloaked cylinder with Re = 8000. The range where reaches 6 cm. A region where does not appear and (b) computational local wake near the cloaked cylinder with Re = 8000. The range where reaches 4 cm. The region where is eliminated.

The optimization algorithms converged to certain local minima in the wake metrics depending on the prescribed initial conditions of the inverse permeability coefficient. Of all of the initial configurations tried in simulations (one or two quadrants, staggered or even obstructions, varying sizes, and spatial frequencies), three staggered downstream obstructions offered the smallest local minimum in both the spherical case and the cylindrical case. Although the DWW of the PVC cases did not drop below that of the exposed cores (control samples), which is particularly evident in Figure 9 showing a larger disturbed region of the flow, the MLW (illustrated in Figures 69) did show notable reductions relative to the control cases. Furthermore, at the tested Reynolds numbers, the size of the time-averaged vortices (the areas constituting the wake region) was dramatically reduced, which is also evident in the local wake field plots.

With the agreement of the velocity and local wake results between simulation and experimentation, other characteristics of the computational results can be examined. The greatest of these are the pressure fields, illustrated in Figures 10 and 11.

Figure 10: The pressure fields between the uncloaked sphere (a) and the cloaked sphere (b) illustrate a localization of the pressure drop in the axisymmetric cross section from the control case to the PVC case. Note the location of the PVC boundaries in (b) and how the pressure immediately downstream region has a notably reduced magnitude of negative pressure.
Figure 11: The pressure fields between the core cylinder (a) and the PVC cylinder (b) illustrate a dramatic reduction in the magnitude of the negative pressure immediately downstream the core. The outline of the PVC boundaries is illustrated in (b).

4. Conclusions

The results of this investigation illustrate that passive structures in the near-field of a given obstruction in constant flow can be configured to control the resulting wake. The near-field zone can be defined as the region where the flow past the structure deviates substantially from its asymptotic behavior. In the low-Re limit, it can be defined as simply twice as large as the structure itself, as suggested by the competition of the far-field (1/) and near-field (1/) terms in the Stokesian wake. The relative magnitude of these near- and far-field contributions can be manipulated with passive obstructions situated in the near-field zone, which is demonstrated by this study.

The results show that the maximum value of the local downstream disturbances can be substantially reduced by employing permeable volumetric composites whose geometry is optimized to meet that specific goal. This is partially due to the reduction or elimination of the large vortex in the flow immediately downstream of the large impervious core, which is a consequence of the near-field array of obstructions of optimized sizes and shapes and positioned at strategic locations. In terms of hydrodynamic “cloaking,” this is a valuable metric as it reduces the possibility of cavitation in the downstream flow. Figure 11 clearly shows that the pressure gradients are localized to the inclusions. Also, examining Figures 69 shows that other elements of the nature of the wake are improved. Although the DWW metric is not improved, the size of the areas of disturbance is reduced in every case except for the experimental cylinder. Areas of significant perturbation from the free-stream flow () are almost eliminated entirely. Finally, this study indicates that significant control over the DWW can best be accomplished through active cloaking methods, as negative permeability is necessary to affect that level of influence on the flow at high Reynolds numbers.

The experiment allowed us to quantify the differences between our CFD models and the actual time-averaged velocity profiles. At 0.27 and 0.24, respectively, for the DWW of the spherical case, a 16% model-measurement discrepancy is observed. Similarly, for the cylindrical case, DWW values of 0.0721 and 0.054 yield a 34% difference. Potential sources of these differences lie in the model set-up, experimental methods, and the data processing. From the computational side, in the case of low-Re RANS models such as Spalart-Allmaras, one may require additional tuning when a Brinkman term is added as an extra volumetric force. The virtual no-slip boundaries, created by regions of low permeability, generate additional viscous boundary layers that are not properly accounted for in the RANS code.

In the experimental set-up, further discrepancy is introduced by having to use additional structural elements whose sole purpose is to maintain the structural rigidity of the sample and to immobilize it relatively to the flow. Both of these sources of errors could be, in principle, quantified by performing full three-dimensional CFD modeling on the CAD geometries of the actual experimental samples, with all connectors and mounting fixtures included, a procedure we propose for future experiments. The experimental results were subject to many other variables, including the effects of the tunnel walls, fabrication accuracy issues, inconsistent seed illumination, and human error from visual calibration of the view sizes.

One final source of error is associated with the nature of the cross section illustrated by the computational models versus the projection viewed by the camera in the experiment. Parts of the out-of-plane physical samples masked the behavior of the flow in some regions very near to the geometric core in the PIV plane of study, making it impossible to gather flow data in those regions. As such, in these regions of question, the computational results will have a known value of two-dimensional fluid velocity whereas the experimental results provide no reliable value. This obstacle makes the DWW measure less reliable from the measurement perspective than the MLW metric, as certain regions of the flow field that were measured and accounted for in the area integral in the computational case were null in the experimental case. This may explain why the MLW discrepancy is consistently smaller than the DWW discrepancy.

In conclusion, we have shown the possibility of wake control and reduction using permeable volumetric composites (PVCs), a generalization of the hydrodynamic metamaterial concept. Our experimental tests demonstrate reasonable agreement with CFD models, with the understood sources of discrepancy. Notably, our demonstrations go beyond the laminar flow regime, reaching Reynolds numbers in the Navy-relevant range of and above. This study paves the way to practically relevant approaches to volumetric wake reduction. Global optimization studies based on our methodology are expected to reveal the maximum extent to which the wake metrics can be manipulated with a PVC of a given volume mounted on a given structure. One possible extension to our methodology is the inclusion of elastic elements capable of generating reaction forces that interact nontrivially with the oscillatory vortex motions. Another extension is the introduction of volumetrically distributed thrusters such as micropump arrays [35] as alluded to in the theoretical studies of Urzhumov and Smith [14, 15].


A. CFD Verification

A.1. Boundary and Boundary Sensitivity

Figure 12 illustrates the agreement of the flow and pressure results around the sphere independent of changes to boundary conditions and domain size. The cylindrical case yields the same insensitivity.

Figure 12: (a) Flow field around a sphere with a symmetric boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (b) Flow field around a sphere with a wall boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (c) Pressure around a sphere with a modified domain size. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (d) Tight view of the flow field around a sphere with a symmetric boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (e) Tight view of the flow field around a sphere with a wall boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (f) Tight view of the flow field around a sphere with a modified domain size. Horizontal and vertical axes are in meters, and the color scale is in meters per second. (g) Tight view of the pressure around a sphere with a symmetric boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in pascals. (h) Tight view of the pressure around a sphere with a wall boundary condition at the far side of the domain. Horizontal and vertical axes are in meters, and the color scale is in pascals. (i) Tight view of the pressure around a sphere with a modified domain size. Horizontal and vertical axes are in meters, and the color scale is in pascals.
A.2. Grid Sensitivity

Figures 13 and 14 illustrate the agreement of the flow and pressure results around the sphere independent of changes to grid. The cylindrical case yields the same insensitivity.

Figure 13: (a) View of the mesh used in the study. The cloak domain was mapped with 15 elements in the larger, inner annular thickness, 8 elements in the smaller, thinner annulus, and 100 elements along the half-circumference of the sphere cross section. This creates a total of 2300 elements. (b) View of a coarser mesh for comparison. The cloak domain was mapped with 12 elements in the larger, inner annular thickness, 5 elements in the smaller, thinner annulus, and 70 elements along the half-circumference of the sphere cross section. This creates a total of 1330 elements, cutting 570 from the model.
Figure 14: (a) View of the pressure around a sphere from the study mesh. Horizontal and vertical axes are in meters, and the color scale is in pascals. (b) View of the pressure around a sphere from the coarser mesh for comparison. Horizontal and vertical axes are in meters, and the color scale is in pascals.

By reducing the number of elements in the cloak domain by nearly half, the results do not perceptibly change.

B. Solver Sensitivity

Figures 15 and 16 illustrate almost nonexistent differences between Spalart-Allmaras and SST for the Reynolds number regime in question and only minor differences between Spalart-Allmaras and -.

Figure 15: Velocity fields around a sphere at according to various turbulent solvers.
Figure 16: Comparison of Spalart-Allmaras to SST (a) and - (b) for flow around a sphere at . Field values are calculated by , where refers to either SST or -.

Competing Interests

The authors verify that no conflict of interests is present with regard to this investigation and their work.


This work was funded by the Office of Naval Research (ONR) through the Naval Undersea Research Program (formerly a University Laboratory Initiative), award no. N00014-13-1-0743. Dean R. Culver acknowledges financial support of the Naval Research Enterprise Internship Program (NREIP), Department of the Navy. The authors are thankful to Maria Medeiros (ONR 333), Charles Henoch, James Hrubes, William Wilkinson, and the rest of the hydrodynamics group at the Naval Undersea Warfare Center Division Newport for their technical support of the experiment and to the personnel of the Center for Metamaterials and Integrated Plasmonics (Duke University) for assisting with the additive manufacturing process.


  1. R. K. Shukla and J. H. Arakeri, “Minimum power consumption for drag reduction on a circular cylinder by tangential surface motion,” Journal of Fluid Mechanics, vol. 715, pp. 597–641, 2013. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  2. F. Fish, “Passive and active flow control by swimming fishes and mammals,” Annual Review of Fluid Mechanics, vol. 38, pp. 193–224, 2006. View at Google Scholar
  3. H. Park, G. Sun, and C.-J. Kim, “Superhydrophobic turbulent drag reduction as a function of surface grating parameters,” Journal of Fluid Mechanics, vol. 747, pp. 722–734, 2014. View at Publisher · View at Google Scholar · View at Scopus
  4. S. Tamano, T. Kitao, and Y. Morinishi, “Turbulent drag reduction of boundary layer flow with non-ionic surfactant injection,” Journal of Fluid Mechanics, vol. 749, no. 1, pp. 367–403, 2014. View at Publisher · View at Google Scholar · View at Scopus
  5. S. L. Ceccio, “Friction drag reduction of external flows with bubble and gas injection,” Annual Review of Fluid Mechanics, vol. 42, pp. 183–203, 2010. View at Publisher · View at Google Scholar · View at Scopus
  6. D. Zhang, W. An, Z. Fan, and L. Liu, “Optimization of underwater body's shape,” Advanced Materials Research, vol. 291–294, pp. 1856–1859, 2011. View at Publisher · View at Google Scholar · View at Scopus
  7. I. Timor, E. Ben-Hamou, Y. Guy, and A. Seifert, “Maneuvering aspects and 3D effects of active airfoil flow control,” Flow, Turbulence and Combustion, vol. 78, no. 3-4, pp. 429–443, 2007. View at Publisher · View at Google Scholar · View at Scopus
  8. O. Bilgen and M. I. Friswell, “Piezoceramic composite actuators for a solid-state variable-camber wing,” Journal of Intelligent Material Systems and Structures, vol. 25, no. 7, pp. 806–817, 2014. View at Publisher · View at Google Scholar · View at Scopus
  9. M. Stadler, M. B. Schmitz, W. Laufer, and P. Ragg, “Inverse aeroacoustic design of axial fans using genetic optimization and the lattice-boltzmann method,” Journal of Turbomachinery, vol. 136, no. 4, Article ID 041011, 10 pages, 2013. View at Publisher · View at Google Scholar · View at Scopus
  10. A. Glezer and M. Amitay, “Synthetic jets,” Annual Review of Fluid Mechanics, vol. 34, pp. 503–529, 2002. View at Google Scholar
  11. P. S. Negi, M. Mishra, and M. Skote, “DNS of a single low-speed streak subject to spanwise wall oscillations,” Flow, Turbulence and Combustion, vol. 94, no. 4, pp. 795–816, 2015. View at Publisher · View at Google Scholar · View at Scopus
  12. M. Skote, “Scaling of the velocity profile in strongly drag reduced turbulent flows over an oscillating wall,” International Journal of Heat and Fluid Flow, vol. 50, pp. 352–358, 2014. View at Publisher · View at Google Scholar · View at Scopus
  13. P. O'Brien and R. Nuzzo, Artificial Cilia. RSC Nanoscience and Technology, The Royal Society of Chemistry, 2013.
  14. Y. Urzhumov and D. Smith, “Fluid flow control with transformation media,” Physical Review Letters, vol. 107, Article ID 074501, 4 pages, 2011. View at Google Scholar
  15. Y. Urzhumov and D. Smith, “Flow stabilization with active hydrodynamic cloaks,” Physical Review Letters, vol. 86, Article ID 056313, pp. 1–5, 2012. View at Google Scholar
  16. J. Pendry, D. Schurig, and D. Smith, “Controlling electromagnetic fields,” Science, vol. 312, no. 5781, pp. 1780–1782, 2006. View at Google Scholar
  17. S. A. Cummer and D. Schurig, “One path to acoustic cloaking,” New Journal of Physics, vol. 9, article 45, 2007. View at Publisher · View at Google Scholar · View at Scopus
  18. P. Spalart and S. Allmaras, “A one-equation turbulence model for aerodynamic flows,” in Proceedings of the 30th Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings, AIAA-Paper, Reno, Nev, USA, 1992. View at Publisher · View at Google Scholar
  19. COMSOL Inc, COMSOL Multiphysics User's Guide, Version 4.4, COM-SOL AB, 2003.
  20. H. Darcy, Public Fountains in the City of Dijon, Victor Dalmont, Paris, France, 1856.
  21. H. C. Brinkman, “Calculations on the flow of heterogeneous mixtures through porous media,” Applied Scientific Research Section A—Mechanics Heat Chemical Engineering Mathematical Methods, vol. 1, pp. 333–346, 1949. View at Google Scholar
  22. J. Rubinstein and S. Torquato, “Flow in random porous media: mathematical formulation, variational principles, and rigorous bounds,” Journal of Fluid Mechanics, vol. 206, pp. 25–46, 1989. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  23. J. Bear, “A one-equation turbulence model for aerodynamic flows,” AIAA-Paper 92-0439, 1992. View at Google Scholar
  24. C. C. Mei and J.-L. Auriault, “The effect of weak inertia on flow through a porous medium,” Journal of Fluid Mechanics, vol. 222, pp. 647–663, 1991. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  25. E. Marušić-Paloka and A. Mikelić, “The derivation of a nonlinear filtration law including the inertia effects via homogenization,” Nonlinear Analysis: Theory, Methods & Applications, vol. 42, no. 1, pp. 97–137, 2000. View at Publisher · View at Google Scholar
  26. Y. Qin and P. N. Kaloni, “Creeping flow past a porous spherical shell,” Zeitschrift fur Angewandte Mathematik und Mechanik, vol. 73, no. 2, pp. 77–84, 1993. View at Publisher · View at Google Scholar · View at MathSciNet
  27. J. H. Masliyah, G. Neale, K. Malysa, and T. G. M. Van De Ven, “Creeping flow over a composite sphere: solid core with porous shell,” Chemical Engineering Science, vol. 42, no. 2, pp. 245–253, 1987. View at Publisher · View at Google Scholar
  28. P. T. Bowen, D. R. Smith, and Y. A. Urzhumov, “Wake control with permeable multilayer structures: the spherical symmetry case,” Physical Review E, vol. 92, no. 6, Article ID 063030, 13 pages, 2015. View at Publisher · View at Google Scholar
  29. M. P. Bendsoe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Computer Methods in Applied Mechanics and Engineering, vol. 71, no. 2, pp. 197–224, 1988. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  30. M. P. Bendsoe and O. Sigmund, Topology Optimization: Theory Methods, and Applications, Springer, Berlin, Germany, 2003. View at MathSciNet
  31. T. Borrvall and J. Petersson, “Topology optimization of fluids in stokes flow,” International Journal for Numerical Methods in Fluids, vol. 41, no. 1, pp. 77–107, 2003. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  32. D. R. Smith, J. J. Mock, A. F. Starr, and D. Schurig, “Gradient index metamaterials,” Physical Review E-Statistical, Nonlinear, and Soft Matter Physics, vol. 71, no. 3, Article ID 036609, 2005. View at Publisher · View at Google Scholar · View at Scopus
  33. C. H. K. Williamson, “Vortex dynamics in the cylinder wake,” Annual Review of Fluid Mechanics, vol. 28, pp. 477–539, 1996. View at Publisher · View at Google Scholar
  34. S. Schmidt and F. Thiele, “Comparison of numerical methods applied to the flow over wall-mounted cubes,” International Journal of Heat and Fluid Flow, vol. 23, no. 3, pp. 330–339, 2002. View at Publisher · View at Google Scholar · View at Scopus
  35. H. T. G. van Lintel, F. C. M. van de Pol, and S. Bouwstra, “A piezoelectric micropump based on micromachining of silicon,” Sensors and Actuators, vol. 15, no. 2, pp. 153–167, 1988. View at Publisher · View at Google Scholar · View at Scopus