Robust design of an asymmetrically absorbing Willis acoustic metasurface subject to manufacturing-induced dimensional variations^a)

The realization of acoustic metamaterials (AMMs)¹ and acoustic metasurfaces (AMSs)² has been enabled, in part, by advanced manufacturing techniques such as additive manufacturing (AM).³ Although the scope of AMM and AMS research is very broad, these advanced acoustic materials can be divided into two classes: active and passive.⁴ Whereas active AMMs derive their unique response from external control systems^5–7 or modulation^8–10 to affect the acoustic propagation behavior, the response of passive AMMs relies solely on material composition and geometry of the structure.^1,3,4 The successful design of the subwavelength structure of an AMM leads to macroscopically observable behavior, which is unachievable in the conventional materials.^1,3 Traditionally, exceptional control of the acoustic wave propagation has been achieved by exploiting a few different underlying physical phenomena, such as Bragg scattering, to construct the phononic crystals,^11,12 AMMs with local resonances,¹³ the generation of anisotropic dynamic density,^14,15 and pentamodal effective properties.^16,17 However, passive AMMs designed to exhibit an asymmetric scattering response^18–21 with respect to the direction of propagation and active AMMs displaying non-reciprocal behavior,^22,23 are particularly interesting for applications in which one desires direction-dependent tunability. Of particular interest to the present study are a class of passive AMMs known as the Willis materials.^18,24 Willis materials are described with nonclassical coupled constitutive relations for the local momentum and pressure-strain behavior, as shown in Eq. (1), due to the presence of the asymmetric subwavelength structure. In the frequency domain, the constitutive relations are expressed as

where $μ$⁠, the momentum density, and $u$⁠, the acoustic particle velocity, are vector fields and the volume strain, $ϵ$⁠, and acoustic pressure, $p$⁠, are scalar fields. In Willis materials, the momentum and strain fields are related by the anisotropic mass density tensor, $ρ_$⁠, adiabatic compressibility, $β$⁠, and Willis-coupling vectors, $η$ and $γ$⁠.²⁵ The unconventional pressure-momentum and strain-velocity coupling provides an additional degree of freedom for the design of materials and interfaces to control wave propagation and scattering behavior, which has been demonstrated in numerous previous publications.^{18–21,24–29} It is important to note that although the constitutive behavior described by Eq. (1) can include non-reciprocal wave phenomena,^7,23,30,31 the present work only considers passive behavior, specifically direction-dependent reflection and absorption of acoustic waves. This is achieved by designing a subwavelength scatterer arrayed in a two-dimensional (2D) surface to create an AMS with Willis properties and asymmetric, frequency-dependent absorption. This interesting behavior is the result of scatterers constructed from lossy multi-material inclusions with an identical orientation in space within the AMS.

Although there has been considerable theoretical work on the topic of Willis materials and some proof-of-concept experimental demonstrations, pushing the boundaries of AMM behavior requires manufacturing methods that can produce the novel material arrangements conceived by engineers. In practice, manufacturing capabilities often restrict the way that material arrangements can be realized. To this end, the ongoing development of AM is a paradigm shift, which enables the realization of metamaterials that were previously impossible to create.³² AM enables the material to be patterned with greater control and for complex features to be built at small scales. To achieve predictable behavior of complex metamaterial designs, the manufacturing method must be able to reliably fabricate them. For a passive asymmetrically absorbing metamaterial, the behavior is sensitive to geometric variation because the physical response relies on resonating inclusions to affect the scattered fields.^33–35 

This work seeks to prescribe the limits on the fabrication accuracy for the creation of Willis-coupled AMSs, displaying asymmetric absorption. Due to the combination of small inclusion sizes and multiple interfaces between different materials, a highly accurate multi-material manufacturing approach is necessary. Multi-material AM is an obvious candidate for use in manufacturing these metamaterials. The existing literature on multi-material AM offers a few potential approaches. Direct ink writing is capable of patterning materials at the submillimeter scale,³⁶ and microfluidic devices are able to deposit inclusions as small as 200 $μm$ in a direct ink writing process,³⁷ but depositing multiple materials within each inclusion is an outstanding challenge. Ultrasonic³⁸ and magnetic³⁹ fields can be used to align and arrange inclusions but they have not been applied to multi-material inclusions. Aerosol deposition methods can create the type of material interfaces required⁴⁰ and have achieved deposition of silver in lines as narrow as 10 $μm$ wide⁴¹ and with very thin layers (on the order of 1 $μm$⁠),⁴⁰ therefore, layering enough material for the inclusions in question may be prohibitively time-consuming. Because none of these methods have been used to create an AMM or AMS with architected features like the asymmetric absorber we seek to design, a robust design approach is sought to help specify the requirements of a manufacturing process that could feasibly fabricate this novel AMS structure to demonstrate a reliable asymmetric absorption performance using existing or future AM technology.

All of the expected variation in the geometry or material properties induced by manufacturing process must be considered during the design process. By doing so, the design task becomes an exercise in manufacturing-aware robust design. Robust design processes seek to identify the designs that exhibit high levels of performance despite underlying variation in the design itself and/or the environment in which it functions.^42,43,45 These approaches account for randomness in subsystems or manufacturing processes as an integral part of the design process when deterministic design optimization is insufficient because it does not account for variable production and operating conditions. Using robust design is especially important for design tasks in which the effect of variations among fabricated components is significant with respect to the expected performance. Robust topology optimization techniques have also been developed to incorporate stochastic variation in material properties and geometry into a topology optimization process,^44,46,47 but they require the underlying analysis to be amenable to gradient-based topology optimization formulations. The approach proposed in this paper is more general and broadly applicable than topology optimization approaches. First applied by the authors to a negative stiffness metamaterial in Ref. 48, it makes use of a unique double-classification procedure to identify the high performance designs and then the subset of those designs that meet the performance targets reliably, even with manufacturing-induced variability, is taken into account. The approach is extended here to investigate the robust design of asymmetrically absorbing AMSs.

This paper presents a class of Willis AMS that consists of an ellipsoidal multi-material inclusion embedded in a background material, which exhibits asymmetric absorption together with a robust design methodology. By designing for robust as-built performance, the expected absorptive response of an additively manufactured AMS can be evaluated against inherent process variability, which is of fundamental importance in realizing metamaterial concepts for real-world applications. Furthermore, from the results of the robust design process, a manufacturer gains the information required to appropriately tune—or even develop—an AM machine to meet the manufacturing requirements for future AMM and AMS designs. In this work, the robust design process is used to ensure the as-manufactured success of a class of asymmetrically absorbing AMSs but is practically extensible to any other AMM or AMS when considering manufacturing variability.

Figure 1 provides a primitive of the AMS design considered in this work. The AMS consists of identically oriented ellipsoidal inclusions embedded in an elastomeric matrix. The inclusions are assumed to be fabricated using a mixture of two materials, one material is a lossless, isotropic elastic solid, such as a metal, and the second material is a lossy isotropic material such as a viscoelastic polymer. These heterogeneous inclusions are embedded in a background material (the matrix) that has a good characteristic acoustic impedance match to the background fluid in which it operates, which is water in this work. The inclusions are assumed to have identical orientation with respect to the global coordinate system and are arranged with square periodicity in the y-z plane such that the AMS is orthogonal to incident waves propagating parallel to the x axis. Figure 1 shows a general representation of the heterogenous inclusions where the distribution of the two materials is graded in the x-direction.

The design of an AMS displaying strongly asymmetric absorption is achieved by simulating the scattering response of heterogeneous inclusions in a matrix material using the finite element analysis (FEA) in the commercial software package Comsol Multiphysics, as detailed in Sec. IV A, and measuring the back- and forward-scattered fields for plane waves propagating in the +x- and −x-directions. The degree of asymmetry in the absorption is determined by performing two numerical experiments—one experiment with excitation via equal amplitude plane waves traveling in the positive x-direction and the second experiment when the excitation is from the negative x-direction—and comparing the backscattered fields from both cases for a wide range of frequencies of interest. The numerical experiments are denoted as cases 1 and 2 for incident wave propagation directions along the unit vectors $n→p,1=1,0,0$ and $n→p,2=−1,0,0$⁠, respectively. The scattered pressure field is then probed on planes that are parallel to the AMS on both sides of the AMS for cases 1 and 2 to determine the absorption from the scattered field as expressed in Eqs. (2)–(4), which are summarized in the following paragraph.^33–35,49 Finally, the geometry of the primitive is varied in ways that are typical of the AM fabrication error until a satisfactory level of asymmetry at the targeted frequency is achieved, which is robust to manufacturing variability.

The plane-waves scattered from the AMS can be represented using a scattering matrix formalism in which the scattering matrix, $S$⁠, relates the incoming (incident) waves to the outgoing (scattered) plane waves.^33–35 When considering the AMS geometry and plane waves propagating along the x axis, this system can be presented by a two-port system using the matrix relationship,

where $pi+$ and $pi−$ are the incident waves traveling in the +x-direction and −x-direction, respectively. Likewise, $po+$ and $po−$ represent the outgoing (scattered) plane waves in the +x-direction and −x-direction, respectively. Given this definition, we can write cases 1 and 2, defined above, as

where we have identified the incoming waves for cases 1 and 2 as $pi,1$ and $pi,2$⁠, respectively, and the outgoing waves are identified for each case with the subscripts “r” and “t” to denote the reflected and transmitted plane waves, respectively, in agreement with the notation provided in Fig. 1. It is now possible to provide an expression for the absorption coefficient of case n as defined by

Note that Eq. 3(b) identifies the coefficients of the scattering matrix as the reflection and transmission coefficients for each case, which are defined as $R1=pr,1/pi,1$⁠, $R2=pr,2/pi,2$⁠, $T1=pt,1/pi,1$⁠, and $T2=pt,2/pi,2$⁠.

Expressing absorption in this form in three dimensions assumes that the incident and scattered wave fields are planar. This assumption is valid in the metamaterial paradigm in which the characteristic inclusion sizes $a$ and inter-inclusion spacing are subwavelength, i.e., $a≪λ$ and $dp≪λ$⁠, and the scattered field is probed multiple wavelengths from the scattering surface such that evanescent fields are extinguished. The absorption asymmetry can then be characterized with the ratio

which is a measure of the dependence of the acoustic absorption on the direction of incidence and will serve as the performance metric for this study. Whereas the performance goal is application specific, one obvious goal would be to maximize the asymmetry as a function of the design features, namely, the excitation frequency, the composition of the heterogeneous inclusions, and the distribution of the constituent materials within the inclusion domain.

The arrangement of dissimilar materials in the inclusions that make up the AMS are essential to generating asymmetric acoustic absorption from the Willis scattering primitive design even when the inclusions of the AMS are deeply subwavelength.^18,33–35 By arranging the inclusions as shown in Fig. 1, the input impedance of the AMS is very different for incident waves from either direction.¹⁸ When one inclusion material is assumed to be lossy, the acoustic absorption may be strongly dependent on the direction of incidence.^33–35 Here, we define the material loss of the inclusion materials by assuming that they are constructed from materials displaying a complex plane wave modulus, $M=M′1−iη$⁠, where $M′$ is the storage modulus, $η$ is the isotropic viscoelastic loss factor, and we have assumed the $e−iωt$ time convention. It is important to note that asymmetry in reflection and absorption does not violate the reciprocity and, thus, the transmission is symmetric, i.e., $T1=T2$⁠. The asymmetry in reflection is, therefore, solely responsible for asymmetric absorption as one would anticipate from previous work.^23,33,34 For this reason, asymmetry of the reflection coefficient $R2/R1$ is also a valid performance metric.

The accuracy of the manufacturing method and the sensitivity of AMS performance must be considered during the design process to be certain that an as-built AMS will meet the as-designed performance criteria. The multi-material AM is an obvious choice for the creation of this type of AMS because of the small size scale as well as the unique material and geometric challenges of building the heterogeneous inclusions to specification.⁵⁰ The multi-material AM processes enable a wide variety of materials and geometries to be explored to create an AMS, demonstrating asymmetric absorption. This work considers combinations of material choices and manufacturing accuracies that are not currently achievable to demonstrate how robust design can be used to help guide process development toward an AM process that is capable of fabricating similar AMMs and AMSs.

A robust design method for AM must account for manufacturing-induced stochastic variations in geometry and the resulting deviations from the target performance of the AMS. By studying the sensitivity of AMS performance to manufacturing variation, the development of a prospective process can directly target the accuracy metrics required to produce an AMS with an as-built performance that meets design requirements. Quantifying those metrics, in turn, reduces the need for experimental iterations in process development. The true performance of an as-built AMM is expected to deviate from the nominal performance, depending on how sensitive the performance is to error (manufacturing variability) associated with the as-built inclusion geometry and the magnitude of that error. To quantify AMS performance, we express the acoustic response as a stochastic function, $g(x,f)$⁠, which has an error component represented as a random variable, $ϵ(z,T)$⁠,

where x represents the design features, $f$ is the frequency, L is the deterministic response, $ϵ$ is the error, $T∈Ω is the AMprocess random variable set$⁠, and $z∈ℝ3$ is the three-dimensional (3D) position in the build chamber. As defined, the error component depends on the interaction between a process-specific set of random variables $T$ as well as the position $z$ in the physical space of the build chamber. The spatial error manifests itself as an inaccuracy of the size, shape, and location of each material and the interfaces between them. There is an abundance of literature attempting to quantify error for specific AM processes,^51–54 as well as translating process error to the as-built part inaccuracy.^44,45,48,55 Regardless of the source of the build accuracy information, a proper robust design workflow must be able to use it.

Once the as-built part variability is quantified, the next step is to translate it into the expected variation in performance and the resulting manufacturability requirements. For the present study, the variation in performance is quantified in terms of the degree of asymmetry of acoustic absorption, $α¯$⁠. This type of analysis is typically performed using a Monte Carlo (MC) approach, but MC methods can be prohibitively expensive when the underlying FEA is computationally expensive. This is especially true when the MC analysis must be repeated iteratively as part of an optimization process. In this work, the computational expense is reduced by replacing the FEA with a classifier-based metamodel of the FEA. The modeling and design for AM are described in detail in Sec. IV.

Identifying a robust design and specifying the manufacturing requirements for this class of asymmetric absorber requires a deliberate design approach to be computationally feasible without requiring supercomputing capacity. Feasibility is a concern because FEA is necessary to evaluate candidate designs, and the potential design space is high-dimensional and vast. Furthermore, the consideration of the robustness of any given design to the manufacturing variability requires additional sampling near optimally performing designs to quantify the sensitivity of performance to manufacturing variation. In lieu of brute force methods, a progressively refining, coarse-to-fine design approach is employed here to address this challenge. A visualization of the design process is provided as a flow chart in Fig. 2. The focus in the present study is on an AMS displaying asymmetric acoustic absorption at a specified design frequency, but the robust design framework presented is applicable for studying the manufacturability of any AMM or AMS that has a quantifiable performance and can be functionally related to the geometric features.

Although the design method seeks to minimize the computational expense, simulating the behavior of the asymmetric absorber using finite elements is a critical element of the design because the spatial variability of the geometry precludes analytical solutions. A thorough description of the FEA specific to this AMS is provided in Sec. IV A. The design method uses FEA to gather data to train a classifier-based metamodel, which closely replicates the input/output behavior of the FEA. For fabrication processes like AM, which have variable process conditions, estimating the success rate for manufacturing a complex component whose performance is sensitive to geometric variability would require many calls to the FEA, making an uncertainty analysis computationally expensive. Alternatively, a metamodel of the simulation input/response relationship can be queried quickly and for a much lower computational expense. Metamodels are surrogate models, which are trained on data obtained from an underlying physics-based model, and approximate the predictions of the underlying model.^56–58 There are many types of metamodels and selecting the appropriate one requires information from studies, benchmark techniques based on the characteristics of the task,^58–60 or directly comparing a suite of techniques for the task at hand. The metamodel used for the asymmetric absorber design is detailed in Sec. IV B. An accurate metamodel then replaces the FEA for MC simulations to evaluate the manufacturability of the asymmetric absorber on a particular machine.

The efficacy of the design methodology for the use of AM for applications in acoustics is demonstrated by designing an AMS with asymmetric absorbing capability. Details of the design process are presented here with results organized together in Sec. V. Some initial performance and design constraints are specified a priori to eliminate superfluous variables and reduce the size of the search space. An analysis is performed on a unit cell consisting of ellipsoidal inclusions inside a matrix material with a cross-sectional length and width of 1.2 mm. The cross-sectional dimensions of the unit cell have been chosen to be much smaller than the acoustic wavelengths in the matrix material at the frequencies of interest (<160 kHz) to avoid Bragg scattering. The characteristic size of the inclusions is necessarily constrained to be subwavelength, $ka≪1$⁠, to ensure that the scattered field is planar in the far-field and we remain in the long-wavelength scattering limit. For satisfactory performance from a design perspective, the asymmetric absorber must exhibit at least a 6 dB difference of the reflection with respect to the direction of the incident wave propagation: $20log10(|R2/R1|)≤−6 dB$⁠. The exploration of the asymmetry with respect to $ka$ is performed by holding $a$ constant (for the manufacturing-agnostic design) with $d=1 mm$ and varying the frequency in the range of 1–160 kHz. Finally, only geometric variations due to manufacturing are included in this design. Material property variability is excluded but will be considered in future work.

The finite element modeling approach and analysis are presented for the general case first because it is called multiple times during the design process with the specific inclusion materials introduced as part of Sec. IV A 3. Figure 3 illustrates the baseline FEA geometry of a unit cell for the case of an incident plane wave propagating in the $+x$-direction.

The geometry for the model consists of a matrix of the background material with a semi-infinite layer of ellipsoidal inclusions having the radii $ra, rb,$ and $rc$⁠, aligned with the x, y, and z axes, respectively, and centered at $xa, xb, xc=(0, 0, 0)$⁠, halfway through the matrix domain with respect to the direction of incident wave propagation. The inclusions are arranged on the y-z plane in a square-packed lattice with center-to-center side-lengths of $dp$ as shown in Fig. 4(a). Each inclusion includes two materials separated by—and bonded along—an interface defined by a $y$-$z$ plane, which intersects the ellipsoid at a distance $h$ from its center [Fig. 4(b)].

The coupling between the pressure acoustics in the matrix and solid mechanics of the inclusion is implemented with the acoustic-structure boundary condition in Comsol Multiphysics. This boundary condition couples the stress of the structure to the acceleration experienced by the fluid matrix. To efficiently model an infinite layer of inclusions, a single unit cell, as shown in Fig. 3, is treated as infinitely periodic in the $x=0$ plane by applying Floquet periodic boundary conditions⁶¹ along both the $y$- and $z$-directions. Elimination of internal reflections in the $x$-direction is performed with perfectly matched layers (PMLs) at either end of the matrix. The PMLs are domains that emulate the Sommerfeld radiation condition by progressively damping acoustic energy over a finite distance.^62–64 With these boundary and domain conditions, the simulated scattering behavior is the result of a square lattice of inhomogeneous inclusions, which constitutes the asymmetric absorbing AMS. The incident wave is a planar background field, which originates between the PML and the inclusions, and propagates toward the inclusions with the matrix serving as a waveguide. The scattered field is measured as a steady state single-frequency result in each separate finite element solution on both sides of the inclusion layer plane. The measurements are taken sufficiently distant from the inclusion layer to ensure that all of the evanescent fields near the inclusions have decayed such that the total scattered field is planar. Because the matrix is lossless, there is no loss of information by extending the distance between the inclusions and probe plane, but for these analyses, this distance is set at $3λ$ from the $x=0$ plane. The field outputs are averaged across each probe plane to account for numerical noise and then post-processed to yield $R$ and $T$ for each simulation.

We note that the physical implementation of this AMS would take the form of a layer of matrix material containing aligned heterogeneous inclusions embedded in another material or fluid. The scattering from the external boundaries of the matrix is application specific and outside the scope of this work. However, we note that the matrix material selected for the case considered here is the silicone-based polymer polydimethylsiloxane (PDMS), which was selected for its impedance match to water. The range of realizable impedance for the unaltered PDMS is dependent on the manufacturing factors but falls between 1 and 1.1 MRayl.^65,66 This leads to the reflection coefficients of approximately 0.15 at a PDMS-water boundary. Additives, such as ZnO and TiO₂, can be used to adjust the density and, therefore, the impedance of the PDMS while minimally affecting the loss such that the water–PDMS interfaces reflect almost no energy.^65,66 The PDMS is treated as a lossless fluid to investigate the absorptive performance of the inclusion layer alone. The shear deformation can be ignored at these frequencies because the shear waves attenuate very rapidly at the frequencies of interest (over 100 kHz).^67,68

Simulating these structures with a subwavelength inclusion layer and sufficient waveguide requires a domain that is multiple orders of magnitude larger than the smallest feature. Consequently, it is both unnecessary and computationally infeasible for the mesh to be made of uniformly sized elements throughout. Rather, mesh density must be higher in and around the inclusion to capture near-field effects. Mesh density is lower where the scattered field is approximately planar with wavelengths on the order of the incident wave. This mesh refinement is implemented by densely meshing the near-field region of the domain within ±2.5 $ra$ and then extruding the mesh to the extents of the remaining domain with a much coarser resolution in the $x$-direction. The characteristic element sizes, $l$⁠, within the near-field domain are bounded as $λ/1000<l<λ/10$ with a maximum element growth rate of 1.3 and curvature factor of 0.2. Comsol automatically controls the resolution of the mesh while using a Delaunay tessellation.⁶⁹ Figure 5 shows the near-field mesh and a portion of the extruded coarse mesh. The extruded mesh scales with the length of the matrix, hence, the element count is nearly identical for all of the simulations at approximately $1.9×105$ elements with approximately $7.6×105$ degrees of freedom.

Another requirement of the mesh is to maintain accuracy along the Floquet periodic boundary conditions. To avoid interpolation and the associated error, the meshes on each periodic boundary are made identical by defining the mesh on one face and copying it exactly to the opposite face. Finding the resolution that efficiently yields accurate results for both of these specific mesh regions and the mesh overall is confirmed with a mesh refinement study before accepting the results from design evaluations.

The design instances are evaluated by parameterizing all of the pertinent geometric variables so that they may be changed programmatically. The simulations are executed with the PARDISO solver.^70,71 The scattered fields are averaged across the probe planes and post-processed to yield the reflection (⁠$R$⁠), transmission (⁠$T$⁠), and absorption (⁠$α$⁠) coefficients. Given all of the constraints, the manufacturing-agnostic ideal design is chosen from a coarse exploratory sample of feature combinations with the asymmetry of reflection, $R2/R1$⁠, as the performance indicator. All of the combinations of the features presented in Table I are included in the initial exploration for a total of 240 unique geometric designs.

The aggregate results of this initial exploration are used to identify the ideal design. From this sample, the most promising ideal design is made of silver and polystyrene with a material interface at $h=200 μm$ for excitation frequencies, $f≅125 kHz$⁠. The strong narrowband asymmetry in the acoustic absorption of this design is shown in Fig. 6. The narrow peaks of the asymmetric response indicate that this design is sensitive to the excitation frequency or, conversely, the performance at the design frequency may be very sensitive to the geometric parameters. The strong asymmetric response of this instance can also be visualized in the 3D plots of the total pressure fields in the matrix and 2D cross sections of the stress in the inclusions. The stress inside the inclusion shows a strong dependence on the direction of the incident waves. When propagating in the +x-direction, much more acoustic energy is converted to energy of deformation in the heterogeneous solid inclusion and, therefore, more energy is absorbed by the viscoelastic material. Recall that the AMS is an interface of finite but subwavelength thickness, separating two half-spaces, and asymmetric absorption occurs in a narrow band of frequencies around a resonance frequency of the inclusions. Due to the asymmetric material distribution within the inclusions, the mode shape of the resonance is also asymmetric and, thus, there are different local impedances for the waves incident from each side. Specifically, on one side of the AMS, the resonance leads to very little motion, whereas on the other side, the local particle velocity and pressure fields better match the characteristic acoustic impedance in the background medium. The asymmetric mode, consequently, leads to poor impedance matching and strong reflections for the plane waves traveling in the +x-direction (case 1), whereas the waves traveling along the −x-direction (case 2) are associated with a good impedance match and strong acoustic absorption (see Fig. 6). This behavior is analogous to the resonances in the detuned Helmholtz resonators reported by Long et al.³³ and Merkel et al.,³⁴ but the present work exploits the material property asymmetry in the subwavelength structures to create an asymmetric mode shape and, thus, asymmetric acoustic absorption with reciprocal transmission. The visualizations of the resulting fields are presented in Figs. 7(a)–7(d).

Like any part whose performance is sensitive to variability, manufacturing an AMS to achieve a strong asymmetry in absorption requires tolerance specifications on the accuracy. It is impractical to design parts like this one to exact feature sizes without accounting for manufacturing variation, especially in the case of AM, which is well-known to suffer from high part-to-part variation.⁷² To account for manufacturing variation, each source of variability can be modeled as a random variable with a specified distribution based on the known fabrication techniques. For AM processes, the distribution that describes the variation of any given feature is often determined empirically via a metrology study. However, the variation in AM processes is inconsistent across machines, raw material types, process parameter settings, etc.⁷² As a result, the success rate of building a satisfactory asymmetric AMS absorber will be different as each process variable changes. For effective and efficient production of these parts, an estimate of the build success rate given a specific set of variable process conditions is invaluable such that the build success rate quantifies the likelihood that a fabricated AMM meets a threshold for satisfactory performance.

For application to the asymmetric absorber design, the ability to specify a minimum magnitude of asymmetry is required. To that end, a classifier is a suitable type of metamodel because the classifiers separate the inputs into classes.⁵⁹ In this work, the classes are separated in a binary sense based on whether the performance of a given design is satisfactory or unsatisfactory. We specify a performance threshold for the acoustic absorption of 6 dB of directional asymmetry of reflection. A shallow neural network (NN) architecture, called a MLP classifier, is used as the metamodel for this application because it is a general purpose nonlinear approximator, which has many options for tuning the performance.⁵⁹ The same classifier would be applicable to any case for which a performance specification requires a minimum (threshold) value to render the part satisfactory.

Generating the data necessary to train the MLP classifier requires querying the FEA sufficiently throughout the space of possible design feature settings. For a classifier that maps the performance to the manufacturing variability, this means starting with the ideal design, estimating the maximum range of manufacturing variability for each feature, and generating a space-filling sample of the variability space. The space-filling sample is generated using a Halton sequence, which is a quasi-random sequence that distributes samples quasi-uniformly throughout the space of the possible feature settings. The Halton sequence is one of many low-discrepancy sample generation methods, which include the Sobol and Faure sequences, as well as Latin hypercube sampling. Although Latin hypercube sampling is more common, the Halton sequence exhibits a lower discrepancy (greater uniformity) over a broad range of dimensions.^73,74 Training the classifier with a space-filling sample improves its accuracy by minimizing the distance between any prospective design instance and a ground-truth point from the FEA.

The ideal silver and polystyrene asymmetric absorber identified in Sec. IV A 3 is being designed for robustness to the spatial variability of the build process. Using the best reported aerosol deposition raster widths as a reference,^40,41 the range of each spatial variable is estimated to be $x±10 μm$ where $x=[xa,xb,xc,ra,rb,rc,h]$⁠. A seven-dimensional Halton sequence is used to generate 500 unique structural designs in this space, which are then evaluated at 20 discrete frequencies in a 2 kHz frequency range spanning 124–126 kHz. This frequency range bounds the strong asymmetric response at 125 kHz in the ideal design. As a result of the spatial variability in the sample, the peak asymmetry of the structures subject to manufacturing variability varies in the frequency and magnitude, as shown in Fig. 8. The training data generation details are in Table II. The data are generated by calling the FEA model for each design at each frequency and post-processing the data to calculate the asymmetric response. The computational expense of generating 10 000 data points is quite high as expected. Approximately 100 computing hours are required on a Linux machine with a 12-core Intel Xeon Gold 5118 central processing unit (CPU) and 250 GB of random access memory (RAM) for an evaluation rate of 100 samples/h via the FEA.

The training data are used to fit the MLP classifier algorithm to accurately represent the FEA outputs to reliably classify the designs that were not explicitly modeled. The MLP classifier is a network (graph) of connected nodes, organized into three types of layers: (i) an input layer, (ii) an output layer, and (iii) one or more hidden layers as shown in Fig. 9.^75,76 The edges and nodes of the network are represented by weights and biases that are multiplied and summed from layer-to-layer such that each output node is a nonlinear function of the preceding layers.

Finally, each output node is scaled by an activation function.⁷⁷ The MLP is trained using backpropagation,⁷⁸ and the mapping between the input and output is optimized using stochastic gradient descent (SGD). These options for the MLP architecture are called hyperparameters. A more detailed description of the MLP classifier algorithm is outside the scope of this work, and the interested reader is referred to a wealth of available literature.^59,75–79

The data from the FEA are used to train the MLP classifier, and the optimization of its hyperparameters is performed using a threefold cross-validation (CV)⁷⁹ scheme, where the prediction accuracy, denoted by ACC, is maximized

where TP and TN represent the true positive and true negative predictions, respectively, and N is the total number of predictions. As part of the CV, the predictive quality of a classifier trained on a subset of data is evaluated for comparison to other classifier training iterations. Training and testing the MLP takes less than a minute on a Windows machine with an Intel i7–8700 CPU and 32 GB of RAM. CV is the most straightforward way to optimize the hyperparameters, and the CV process also helps avoid overfitting by training and testing the MLP on multiple unique slices of the training data. The specified hyperparameter settings are tested for all of the unique combinations. For this design, the optimized hyperparameters are the hidden layer quantity and size, activation function, the weight optimizer, and regularization parameter (⁠$λw$⁠). Because training and testing the model requires very little computational expense, a large hyperparameter space was explored. The hyperparameter search space and optimal settings are displayed in Table III. With optimized parameters, the training accuracy of the classifier is repeatably 98%.

Once the classifier is trained to a high accuracy, it is used to evaluate the success of manufacturing the ideal design on a machine with a nonideal size and position resolution (i.e., spatial variability). The trained classifier is called 1 × 10⁶ times by the MC simulation to evaluate the expected success rate of the parts subject to the manufacturing variability. 1 × 10⁶ calls to the trained classifier take approximately 2 s on the same machine used to train the classifier. The success rate of this design instance is evaluated for multiple different idealized AM machines, where the geometric variables from Table II are distributed random variables. These random geometric variables are characterized by multivariate truncated Gaussian distributions and multivariate beta distributions. The truncated and beta distributions are chosen to represent the centered and biased distributions while being bounded. As a matter of practicality, the distributions must be bounded to prevent the inclusion geometry from being built outside of the matrix. The biased beta distributions represent a manufacturing process that commonly skews one way (i.e., shrinkage or over-sintering).

The truncated Gaussian distributions are truncated at ±1$σ$⁠, $±2σ$⁠, and ±3$σ$ for $μx=0 mm$⁠, $μr=0.5 mm$⁠, $μh=0.2 mm$⁠, and $σx,r,h∈1,3,5,8,10,15,20 μm$⁠. The range of standard deviations, $σ$⁠, was selected such that the lower bound matches the thinnest known layer thickness acheived from the aerosol deposition techniques,⁴⁰ whereas the upper bound represents double the raster width.⁴¹ The beta distributions are defined as $Beta(2,6)$ and $Beta(6,2)$ with the lower bound (LB) and upper bound (UB) defined as $LB,UB=μ, μ+10 μm$⁠. The probability density functions (PDFs) of these distributions are shown in Figs. 10(a)–10(d) and 11 for a single variable. For a real AM machine, the variability of every spatial feature would likely be described by a different probability distribution. For the sake of this demonstration, the distributions and standard deviations for each variable are assumed to be defined identically [e.g., $xa,b,c ∼ Nμx,σ$⁠, $ra,b,c ∼ Nμr,σ$⁠, $h ∼ Nμh,σ$], and the frequency is sampled uniformly, $f ∼ U(124 kHz, 126 kHz)$ for all of the MC simulations.

Finally, to validate the metamodel, a subset of 500 designs, generated by each MC simulation, are evaluated with the FEA and compared with the classifier metamodel's predicted performance. The validation step will ensure sufficient accuracy of the classifier for a specific set of process conditions.

The results of the MC simulation using the MLP classifier metamodel provide important information for the prospective manufacturing of this asymmetric absorber. Most clearly, it shows the proportion of the as-built parts with dimensions specified by the idealized design, which would be expected to achieve a satisfactory magnitude of asymmetry. Conversely, the required manufacturing variability for an AM process can be specified to reach a satisfactory as-built success rate. As a result, the manufacturer gains information to balance the risk of part failure with the cost of accuracy. For the example case considered here, the predicted success rates are higher for lower variability, as one would expect, with a success rate of approximately 40% for a machine with a $σ=1 μm$ accuracy for all of the specified geometric variables. Conversely, running a machine with $σ=20 μm$ would lead to only about a 14% success rate. The success rate results for all of the ±1$σ$ truncated Gaussian distributions used in the MC simulations are shown in Fig. 12(a). Figure 12(b) shows the results for ±1 $σ$⁠, $±2σ$⁠, and ±3$σ$ truncated Gaussian distributions. The results of a validation subset, generated by evaluating the FEA simulation, are also shown in Fig. 12 and labeled as “true.”

The analysis of the validation samples shows that the classifier is ≥95% accurate when the samples are generated from the MC simulations with ±1 $σ$ truncated Gaussian distributions, as documented in Fig. 13(a). As shown in Fig. 13(b), the classifier accuracy for the ±3$σ$ truncated Gaussian distribution with $σ=10 μm$ is very high at 98.2% despite no designs with a feature size of $μ±30 μm$ being included in the training set. Although very few samples are taken from the tail probability near $σ=±30 μm$⁠, it is notable that the metamodel is effectively extrapolating in these cases.

In contrast to the non-biased Gaussian distributions that concentrate the as-built instances near the ideal design, the beta distributions are used in the MC simulations to show the effect of bias in a manufacturing process. This case represents a machine with a process parameter maladjusted, leading to the mean spatial result being skewed away from the expected mean. When the geometric variables are skewed toward the design setting with $Beta2,6,$ the manufacturability performance is very similar to the Gaussian distribution results. However, when they are skewed away from the design setting and closer to the upper bound of $10 μm$ with $Beta6,2$⁠, there are no satisfactory as-built instances identified in the validation set. The classifier accuracy is still good for both cases; it even successfully identifies all instances as unsatisfactory when the variables are distributed as $Beta6,2$⁠. The results of the skewed MC simulations appear in Fig. 14.

Ideally, the accuracy of a metamodel is the same for any spatial feature distributions tested. The demonstration problem shows very strong accuracies of 95%–99% as long as the metamodel is fit to a range of data such that the studies require interpolation only. Whereas some extrapolation may be possible, the metamodel will lose accuracy very rapidly with respect to extension beyond its training data range. For this reason, the selection of an appropriate training data range is important. The range ±10 $μm$ was selected for training in this demonstration because it became clear during the initial design exploration that very few samples were satisfactory for the larger geometric feature ranges at the frequencies of interest.

Clearly, the performance of this structural asymmetric absorber is very sensitive to the spatial variation during manufacturing. This is evident in that the manufacturing success rate is only 40% with a modest standard deviation of $σ=±1 μm$ on all of the geometric features. The success rate declines nonlinearly as the distributions widen, whether due to more tail probabilities included (⁠$±1σ$ to ±3$σ$ truncations) or a larger $σ$⁠. The success rate decline flattens as the truncated Gaussian distributions get wider because they are all centered on the mean value of the ideal design. Comparatively, the metamodel predicts zero instances of successful builds for the beta distribution, which is skewed away from the mean. This shows the value of the unbiased error in manufacturing these metamaterials. Although it may be unlikely that all of the spatial features skew that much from the ideal design dimensions, it demonstrates a worst-case scenario and shows that the metamodel's predictive capability is equally strong in the unsatisfactory regions of the design space.

This work has demonstrated that the fabrication of an AMS with a Willis constitutive behavior type with asymmetric absorbing capabilities and well-matched to water is possible using subwavelength asymmetric inclusions in an elastomeric matrix material. However, for the design to function properly at a design frequency of 125 kHz requires strict levels of accuracy in the fabrication process. The information provided here lends itself to a few options. Most obviously, the size of the structure could be increased. However, increasing the size will necessitate designing such that the lower frequency incident waves stay in the metamaterial regime (⁠$ka≪1$⁠) and maintain approximately planar scattered fields. Increasing the size of the metamaterial may also limit its applicability as its size and weight will increase accordingly. Another option is to develop a rigorous multi-material AM process that is suitable for building the AMS so it can work at higher frequencies. The methodology developed in this work will provide the specifications of manufacturing process accuracy required to build this type of asymmetric absorber for any desired operating frequency and size scale. At the time of writing, aerosol deposition development has achieved a $1 μm$ size resolution⁴⁰ but only in the direction normal to the deposition surface. Aerosol deposition looks promising for building inclusions, although future work will need to more fully characterize the spatial variability of the process as well as address the co-deposition of the viscoelastic polymers. Additionally, the layer-wise deposition of many layers is difficult due to the thin layer sizes and layer deposition abrading the previous layer.

In addition to the spatial parameters, constitutive material properties are a source of manufacturing variability. This study only considered the spatial parameters to reduce the complexity and improve the clarity of the method and demonstration. In practice, the material properties are also random variables, which are dependent on AM process parameters. For example, the material properties of PDMS are functionally related to resin viscosity, curing conditions, additives, material aging, etc. Metals produced with AM processes often have material properties that differ from reference values and are usually anisotropic based on build directions. For these reasons, the variability of the as-built material properties cannot be ignored when building sensitive AMSs such as the asymmetric absorber considered here. Fortunately, as the AM industry grows and looks toward the future, many organizations are characterizing their machines and performing exactly the type of variability studies needed to execute a classifier-based robust design with the MC simulation.⁸⁰ However, there will continue to be a need for design methods that help designers predict and manage the effect of process-induced variability on acoustic material performance.

Although this work focuses on a particular AMS with a single layer of elliptical inclusions, the design methodology presented is applicable to the design of any structure whose performance is sensitive to spatial or material variability, which is characteristic of manufacturing with AM machines. The reusability of the metamodel is significant to reducing the computational expense of evaluating the robustness of manufacturing the part on an AM machine, but it is crucial that the designer understand the limits on reusability. A metamodel can be reused when the variability of the manufacturing process stays within or very near the range of the data used to train the metamodel. Any significant changes in the geometry will require training of a separate metamodel. No metamodel is universal, but the method presented here satisfies the need to evaluate many machines given a target design and performance specification.

The authors acknowledge support from the National Science Foundation EFRI under Award No. 1641078 and Office of Naval Research through YIP Grant No. N00014-18-1-2335.