#### Abstract

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.

**(a)**

**(b)**

**(c)**

**(a) Re = 20**

**(b) Re = 100**

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

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.

**(a)**

**(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.

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.

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.

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

**(a) Exp**

**(b) Comp**

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

**(a) Exp**

**(b) Comp**

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

**(a) Exp**

**(b) Comp**

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

**(a) Exp**

**(b) Comp**

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 6–9) 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.

**(a) Uncloaked**

**(b) Cloaked**

**(a) Uncloaked**

**(b) Cloaked**

#### 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 6–9 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].

#### Appendix

#### 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

##### 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.

**(a)**

**(b)**

**(a)**

**(b)**

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 -.

**(a)**

**(b)**

#### Competing Interests

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

#### Acknowledgments

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.