With the advance of additive manufacturing, many researchers are increasingly interested in planar acoustic lenses that are not only easier to fabricate than typical convex/concave lenses, but also have excellent imaging performance. However, the planar acoustic lenses reported so far cannot work for a short-duration pulse used in conventional imaging systems due to their inherent dispersive characteristics. This study addresses the challenge by devising a transient topology optimization formulation to design a planar acoustic lens that works effectively for a short-duration pulse. A planar lens consists of two materials where optimal combination and distribution are obtained with a crisp interface via the level-set method. Design is based on the transient acoustic responses, which are calculated from a time-dependent acoustic model solved by the Newmark method. The proposed method uses the area-fraction approach to compute the acoustic properties of a cut element by the interface. A localizing time-window function is introduced so that acoustic energy can be focused within the desired time range as much as possible. We obtain optimum design solutions designed with the proposed method and verify its effectiveness through the numerical investigations.

## I. INTRODUCTION

Imaging systems such as a scanning acoustic microscope (SAM) image an object by measuring the reflection (or echo) waves caused by the impedance mismatch with the background medium (Yu and Boseck, 1995). The clear separation of incident and echo waves as well as the generation of a large echo wave play an important role in realizing a satisfactory imaging performance. To achieve this goal, conventional imaging systems use acoustic energy focusing devices such as acoustic lenses (Che *et al.*, 2015; Khuri-Yakub, 1993; Lemons and Quate, 1974). Acoustic lenses are generally classified into two types according to their shapes: (1) convex/concave and (2) planar. A convex/concave acoustic lens focuses acoustic energy by adjusting the phase delay of the transmitted waves through the curvature of its surface (Labréche *et al.*, 1985; Yang *et al.*, 2017; Yokosawa, 1996). Due to its high transmission efficiency, conventional imaging systems have primarily utilized this type of lens. However, the dimensions are typically required to be large compared to an operating wavelength, and the focusing performance is very sensitive to the polishing condition of the lens surface. In contrast, a planar acoustic lens, such as a Fresnel zone plate (FZP) (Calvo *et al.*, 2015) or a super-oscillatory acoustic lens (SOAL) (Hyun *et al.*, 2018), utilizes the constructive interference by diffracting the transmitted waves through its micro-slits. This lens type has several favorable features in that it is compact, is easy to fabricate, and has a good lateral imaging resolution. In addition, with the recent development of metamaterials/metasurfaces and additive manufacturing (AM), planar acoustic lenses with tailored characteristics are becoming more viable (Xie *et al.*, 2014).

The planar acoustic lenses reported so far have dispersive characteristics that make them difficult to apply to actual imaging systems, in which the focal position depends on the operating frequency (Chen *et al.*, 2019; Hyun *et al.*, 2018; Shen *et al.*, 2019; Tarrazó-Serrano *et al.*, 2019). The dispersive characteristics cause a narrow operating frequency range, making it challenging to use the planar acoustics lenses in many applications. Meanwhile, in efforts to use such lenses in imaging systems, either an *N*-cycle tone burst input pulse or an additional deconvolution method has been employed for generating the images of target objects through signal processing (Shen *et al.*, 2019). The *N*-cycle tone burst pulse approach can be easily utilized, but the number of cycles of the pulse should be large enough so that its center frequency approximates the operating frequency of the lens, making it difficult to clearly separate incident and echo waves. The deconvolution method can extract the information of the imaging object from the measured echo waves but has the disadvantage of requiring knowledge of the point spread function in advance, making it difficult to apply real-time imaging due to the additional signal processing. Thus, to achieve clear separation between the incident and echo waves in the time domain, a short-duration pulse is suggested to be an ideal approach (Parker *et al.*, 2010). However, if this type of pulse is utilized with existing planar acoustic lenses, most of the acoustic energy within the operating frequency range cannot be focused on the imaging object and is wasted due to the dispersive characteristics. Therefore, the design of a planar acoustic lens working with a short-duration pulse remains a research challenge.

A short-duration pulse has the characteristics of a broad frequency band. Therefore, one way to design a device (e.g., lens) working with a short-duration pulse is to make the focal position constant for the frequency range of the target short pulse, achieving broadband focusing performance (Hyun *et al.*, 2020). Since this approach considers the pulse in the frequency domain, not in the time domain, it is referred to as an indirect method. However, this approach is computationally expensive because a large number of frequencies must be selected for the design to achieve a satisfactory broadband performance. Also, it is difficult to consider a generalized time-dependent pulse, even though this type of pulse can be handled in the frequency domain through the fast Fourier transform (FFT). Hence, transient responses must be taken into account “*directly*” to design an acoustic lens working with the generalized time-dependent pulse including a short-duration pulse.

With the development of AM technology, topology optimization (TO) has attracted significant attention in recent years and has been applied to the design of acoustic energy focusing devices including lenses. The multilayered convex/concave acoustic lens was topologically optimized to minimize the comatic aberration problem that the focal position is changed according to the angles of incident waves. This was achieved by making the focal position the same for the multiple incidence angles, which is beneficial because it allows the imaging object to be made on the same imaging plane (Tran *et al.*, 2017). The layout of an ultrathin planar acoustic lens, the SOAL, was topologically optimized and experimentally verified to achieve ultrasonic sub-diffraction focusing performance by maximizing the acoustic intensity at a target focal position, being subject to a constraint of the specified beam width (Hyun *et al.*, 2018). In this study, the sub-diffraction focusing performance that overcomes the diffraction limit was achieved for the first time through a planar acoustic lens. Yoon *et al.* (2018) recently reported that a set of annular piezoelectric rings could be designed using TO to focus the acoustic energy in the target focal position. These studies designed focusing devices (lenses) working at a single operating frequency through the time-harmonic analyses. However, to the best of our knowledge, TO has not been developed for planar acoustic lenses considering the transient responses.

The objective of this study was to develop a new TO formulation and methodology for transient acoustic wave propagation and apply it to the inverse design of a planar acoustic lens working with a short-duration pulse. A boundary-based TO method called the level-set topology optimization (LSTO) method is applied to design a lens with crisp boundaries such that the resulting design can be realized in practice. The design sensitivity for updating a structural boundary is calculated via the adjoint variable method (AVM). The effectiveness of the proposed method is demonstrated by optimizing planar acoustic lenses for various design conditions, including incident pulse type, lens size, target focal length, and materials used, and the design solutions are analyzed and discussed.

## II. TRANSIENT ACOUSTIC PHYSICAL MODEL

### A. Time-dependent acoustic wave equation and boundary conditions

As mentioned earlier, a recent study clearly stated that the transient response should be considered for achieving optimal focusing performance of an acoustic lens; otherwise, the focusing performance, including focal intensity and lateral and axial resolutions, is degraded (Pérez-López *et al.*, 2020). Therefore, for the sake of demonstrating this importance more clearly and simply, this study assumes that only compressional (or longitudinal) waves are considered. In this case, in the linear structural elastic equation, the effect of shear modulus is neglected, and only the effect on bulk modulus is considered, so that the purely scalar wave equation can be adopted as the governing equation. A couple of studies on acoustic device designs considered the same assumption for polymeric materials, such as rubber, Stycast 1090SI, and acrylic plastic [poly(methyl methacrylate) (PMMA)] materials (Lu *et al.*, 2013; Tran *et al.*, 2017; Xia *et al.*, 2013). Assuming inviscid and irrotational fluid, the time-dependent acoustic wave equation is expressed as Eq. (1) with the given initial conditions,

where $p$ is the acoustic pressure (Pa), $\rho $ is the mass density ($kg/m3$) of medium, $\kappa $ is the bulk modulus (Pa) of medium, and $Qpt$ is the point (monopole) source, $\u2212\u2202qpt\u2202t\delta x\u2212xp$, with $\delta $ being the Dirac-delta function, $x$ the global coordinates, and $xp$ the specified source position. Here, $qpt$ is the volume velocity defined on a unit volume. Four types of boundary conditions are usually considered in an acoustic problem:

*Prescribed pressure*(2)$p=p0t\u2003on\Gamma p,$*Normal velocity*(3)$\u2212n\u22c51\rho \u2207p=\u2212\u2202vint\u2202t\u2003on\Gamma v,$*Surface impedance*(4)$\u2212n\u22c51\rho \u2207p=1Zst\u2202p\u2202t\u2003on\Gamma z,$*Radiation*(5)$\u2212n\u22c51\rho \u2207p=1Zo\u2202p\u2202t\u22122\u2202pint\u2202t\u2003on\Gamma rad,$

where $n$ is the unit normal vector, $p0$ is the prescribed pressure, $vin$ is the input normal velocity, $Zs$ is the surface impedance, $Zo$ is the characteristic impedance of medium, which is the water in the present study, and $pin$ is the input pressure. Applying the weighted residual method, the following weak formulation of the time-dependent acoustic wave equation can be derived (Bathe, 1996):

### B. Finite element discretization

In this study, acoustic pressure field $p$ and its gradient are approximated via the isoparametric quadrilateral finite elements (FEs) with bi-linear shape function, $N$, as expressed in Eqs. (7) and (8),

where $N$ is the shape function, $B$ is the gradient of the shape function, and $p$ is the nodal acoustic pressure vector of an element. The weak formulation is discretized through the standard Galerkin approach, in turn, leading to the transient FE equation in a matrix form as

and

where *NE* is the total number of finite elements, and *e* is each element index.

### C. Transient response analysis

The transient FE equation of (9) is solved by the implicit Newmark method, which is one of the most stable direct time integration methods. The transient FE equation is discretized for the time step $tn+1$,

To solve Eq. (15), the updates of the first and second derivatives of pressure fields, $p\u0307n+1$ and $p\u0308n+1,$ at the time step $tn+1$ should be defined, respectively, as

where $\beta 1$ and $\beta 2$ are the control parameters, set to 0.25 and 0.50, respectively, in this study (Bathe, 1996), and $\Delta t$ is the time step. Substituting Eqs. (16) and (17) into Eq. (15), the motion of equation is obtained, Eq. (18), and the time-dependent pressure fields $p(t)$ can be calculated by solving this equation. For further details, readers are referred to Bathe (1996),

where

### D. Gaussian-modulated pulse: Simulating short-duration pulse

Many acoustic imaging systems use a short-duration pulse as the incident pulse to clearly separate the incident and echo waves without their overlapping. This short-duration pulse is typically modeled as a Gaussian-modulated pulse, expressed as Eq. (19) (Dahl *et al.*, 2008),

where $f0$ is the center frequency, and $t0$ is the center time of the pulse. The $bw$ defines the fractional bandwidth defined as the frequency range from the center frequency $f0$ to the frequency that intensity drops below a threshold of −60 dB compared to peak level. To determine an appropriate $bw$ value, we first obtain the frequency response of a pulse with an arbitrary $bw$ value through the FFT and subsequently select the $bw$ value by comparing it with the frequency range of interest. The Gaussian-modulated pulse is sampled according to the time step $\Delta t$ and then applied to the radiation boundary condition, thereby constructing $Frt$ and $Crt$ in Eq. (9).

## III. LSTO

In the present study, we adopt the LSTO method (Dunning and Kim, 2015; Kambampati *et al.*, 2020), as the systematic design approach for a planar acoustic lens. Compared to a typical density-based TO method such as SIMP (Bruns, 2005), one of the advantages of LSTO is a crisp boundary (or interface) of an optimized structure without gray regions or intermediately mixed materials. This section will explain both how the LSTO method updates a structural boundary and the problem formulation.

### A. Structural boundary update via the level-set method

The level-set method employs an implicit representation of the structural boundary,

where $\Omega s$ and $\Omega v$ are the structural and void domains, respectively, $\Gamma $ is the structural boundary, $\varphi \chi $ is typically an implicit signed-distance function, and $\chi $ is the point inside the design domain. The structural boundary is moved under a given design velocity field $Vn\chi $ using the following level-set advection equation, the so-called Hamilton–Jacobi equation (Sethian and Vladimirsky, 2000),

This equation is solved numerically using the following discretized equation:

where $k$ is the design iteration number, $i$ is the index of a grid point, $\Delta tp$ is the pseudo-time step, $Vn,i$ is the design velocity at the grid point $i$, and $\u2207\varphi ik$ is computed using the Hamilton–Jacobi weighted essentially non-oscillatory (HJ-WENO) scheme (Sethian, 1999). Here, a generic optimization problem can be stated mathematically as

where $J$ is the objective function and $gm$ is the *m*th inequality constraint. For a given level-set topology, these objective function and constraints can be linearized with the help of the shape sensitivities, thereby constructing a sub-optimization problem.

### B. Linearized sub-optimization problem

It can be seen from the discretized Hamilton–Jacobi equation, Eq. (22), that to update the implicit level-set function $\varphi i$, the boundary velocity is required. This is obtained by solving the following problem linearized by the first-order Taylor expansion:

where $J$ is the objective function, $gm$ is the *m*th constraint with constrained value $g\xafm$, $Vn,j\Delta tp=\Delta \Omega j$, and $dcfl$ is the maximum distance the boundary point can move along the normal direction inside the design domain, i.e., the Courant–Friedrichs–Lewy (CFL) stability condition. The integrals in Eq. (24) are estimated by

where $lj$ is the discrete length of the boundary segment around the boundary point $j$, $SJ$ and $Sgm$ are vectors containing the product of boundary lengths and boundary sensitivities, and $Vn$ is the vector of normal velocities. With Eq. (25), $Vn$ can be obtained to minimize the objective function, given any constraints. This provides a numerical implementation of the steepest-descent direction of the objective function with the constraints for the set of $Vn$ by introducing the Lagrange multipliers,

where $\lambda J$ and $\lambda gm\u2009(m=1,\u20092,\u2009\u2026,Ng)$ are the design variables to be determined and the Lagrange multipliers mathematically. The problem Eq. (24) is then reformulated as a resulting linearized sub-optimization problem,

where $\Delta g\xafm$ is the allowed change of the *m*th constraint function at every design iteration. The boundary velocities are calculated by solving this linearized sub-optimization problem, and subsequently the velocities are extended through the fast velocity extension algorithm to those at grid points within the narrow band close to the structural boundary. Then the level-set function is updated using the extended design velocities through the discretized Hamilton–Jacobi equation, Eq. (22). Since our LSTO method uses the narrow band update scheme, the frequent re-initialization of the level-set function to maintain the signed-distance function property is not required. Our LSTO method has the following features: (a) reduced re-initialization frequency through the fast velocity extension and (b) reduced number of design variables through the sub-optimization problem linearized by shape sensitivities. The sub-optimization problem is solved at every design iteration via the simplex algorithm (Oki, 2012), yielding the optimum boundary point velocities. Further details can be found in Dunning and Kim (2015).

### C. Objective function

A generic objective function $J$ of a transient optimization problem is to either minimize or maximize the integral of a specific function of interest $f$ in the time interval from $Tinit$ to $Tend$. Here, setting $Tinit$ to 0 in general, $J$ is expressed as

To optimize a planar acoustic lens working with a short-duration pulse, the objective is to maximize the integral of the squared pressure $pobj2$ at the target focal position in a specified time interval, [0, $Tend$]. In the present study, $Tend$ is set to a sufficiently long time to realize the satisfactory focusing performance of the lens. The logarithmic function is employed for better numerical scaling. The localizing time-window function $ht\u2212tc$ is introduced, where $tc$ is the center time of the function,

### D. Localizing time-window function

This study introduces a localizing time-window function, $ht\u2212tc$, which maximizes the focusing performance in the desired time range and minimizes it in other time ranges. This prevents a possible unwanted focusing outside the target time range. The significance of the time-window function is twofold: (a) controlling a specific time delay and (b) localizing acoustic energy at the target focal point, making it easily applicable to acoustic lens design. In this function, the center time $tc$ is introduced to specify a specific time range in which focusing is achieved, and it is different from $Tend$, which defines the total analysis time in the objective function of Eq. (29). Figure 1 shows the behavior of this function expressed by the sum of the two relaxed Heaviside functions, $H\u0303t\u2212tlower$ and $H\u0303t\u2212tupper,$

where

Here, $\Delta tc$ is set equal to the time-bandwidth of the incident pulse, and the center time $tc$ is determined as the average between the fastest time $\Delta tfastest$ for the wave paths to reach the target focal position from the input source and the slowest time $\Delta tslowest$, as illustrated in Fig. 2. The $\Delta tfastest$ and $\Delta tslowest$ are calculated using the simple relationship between speed and propagation time, i.e., speed = frequency × wavelength. Also, since acoustic waves propagate in the order of water, lens, and water, they are calculated by considering the approximate average speed of the lens domain. The effect of the localizing time-window function on the optimized lenses will be discussed in more detail in Sec. VI C,

## IV. SENSITIVITY ANALYSIS

### A. Material interpolation functions

The proposed LSTO method employs the fixed Eulerian grid, and an interface can cut through the grid cells, as shown in Fig. 3. For the finite element method, the material property of a bi-material finite element is approximated using the rational approximation of material properties (RAMP)-like interpolation functions (Lee and Kim, 2009),

where $\rho 1$ and $\rho 2$ denote the mass density of the design materials 1 and 2, respectively, and $\kappa 1$ and $\kappa 2$ denote the bulk modulus. $\alpha e$ is the area-fraction of cut element *e* and defined as

where $\Omega 1,e$ is the area of cut element *e* that lies inside the material-1 phase, and $\Omega e$ is the total area of the element. As illustrated in Fig. 3, the LSTO method finds the distribution and combination of design materials 1 and 2 in the lens domain based on the area-fraction $\alpha e$.

### B. Sensitivity computation via the AVM

The shape sensitivity of an objective function, $J$, with respect to the boundary $\Omega $ is required to update the level-set functions. We first derive the element density sensitivity $dJ/d\alpha e$ by applying the adjoint sensitivity method. The corresponding augmented Lagrangian function $J\u0303$ is formed by introducing an adjoint variable vector, $\lambda t,$ as

Hereafter, $t$ is omitted for convenience. Differentiating the augmented Lagrangian function defined above with respect to the area-fraction $\alpha e$ using the chain rule, it gives

Applying the integration by parts for the last term in Eq. (36), both $\u2202p\u0308/\u2202\alpha e$ and $\u2202p\u0307/\u2202\alpha e$ are converted to the functions of $\u2202p/\u2202\alpha e$, and considering that $p0$ and $p\u03070$, i.e., initial conditions, have no dependence on design, Eq. (37) is obtained,

The adjoint equation of Eq. (38) with the condition of $\lambda Tend=\lambda Tend=0$ is obtained from the second term in the right-hand side in Eq. (37) to compute the adjoint variable $\lambda $,

We convert the terminal value problem of Eq. (38) to the initial value problem, Eq. (39), through the time-mapping $\tau =Tend\u2212t$. The Newmark time integration scheme used to solve the primal equation can also be applied to the adjoint equation. Here, the notation $Tend\u2212\tau $ means that the adjoint load vector $\u2202f/\u2202p$ should be applied as opposed to the time integration direction of acoustic pressure $p$ obtained by the primal equation,

$\lambda \u0303$ obtained by solving the adjoint equation is converted back to $\lambda $ through the time-mapping relationship, and by substituting it into Eq. (37), the resultant sensitivity formulation is expressed as

Now applying the material interpolation functions and the element matrices, the sensitivity of the generic objective function $J$ can be finalized as

where $M0$ and $K0$ are the unit element matrices, respectively. For our objective function of Eq. (29) containing the scaling logarithmic function, the element density sensitivities are expressed as Eq. (42) by the chain rule,

In this equation, $f$ is the term $pobj2t\u22c5ht\u2212tc$ in the time integral, and thus the adjoint load vector $\u2202f/\u2202p$ is expressed by

where $L=0,\u2009\u2026,\u20090,\u20091,\u20090,\u2009\u2026,\u20090\u2009T\u2208R1\xd7n$, with 1 at the target focal position entry and 0 elsewhere.

As mentioned in Sec. III, the shape sensitivity should be calculated to update the level-set function. In our LSTO method, the structural boundary is discretized into a number of boundary points by zero level-set, and the shape sensitivity measures the amount of change in the function of interest $J$ when these boundary points change in their normal direction. As a result, the shape sensitivity $SJ$ is expressed as a vector set containing the shape sensitivities evaluated at the discretized boundary point $j$ as expressed in Eq. (44). Here, $\Omega j$ means the *j*th discretized boundary point itself,

where *j* is the index of boundary points, and NB is the number of boundary points. After calculating the element density sensitivities $dJ/d\alpha e$ at Gauss points in each element *e*, the weighted least-squares interpolation method is used to obtain the shape sensitivities $dJ/d\Omega j$ (Picelli *et al.*, 2018) by Eq. (45),

## V. OPTIMIZATION METHODOLOGY

The above procedure has been implemented in the object-oriented LSTO framework using c++ (Kambampati *et al.*, 2018). The overall methodology is summarized as follows:

Initialize the FE analysis and level-set function objects, including material properties and optimization parameters.

Initialize the level-set function $\varphi $ in the design domain with initial holes.

Discretize the structural boundary into

*j*boundary points via the current implicit level-set surface.Compute the area-fractions $\alpha e$ of cut elements with Eq. (34).

Update the material properties of each element using the material interpolation functions, Eqs. (32) and (33).

Assemble and solve the transient primal FE equation using the implicit Newmark method to obtain the acoustic pressure field $pt$ using Eq. (18).

Calculate the objective function, Eq. (29), via the trapezoidal rule.

Assemble the adjoint load vector using Eq. (43) in time domain.

Solve the transient adjoint FE equations, Eq. (39), with the adjoint load vector using the implicit Newmark method to obtain the adjoint variable vector, $\lambda t$.

Calculate shape sensitivities at the boundary points via the weighted least-squares interpolation scheme.

Compute the design velocities at the boundary points.

Update level-set function $\varphi $ at the grid points using Eq. (22).

Check the convergence by checking the termination criterion, e.g., the maximum iteration or the change in the objective function during the last five iterations.

Iterate from step (3).

## VI. NUMERICAL INVESTIGATION

### A. Problem definition

Figure 4(a) shows the configuration of a numerical model for optimizing planar acoustic lens, which consists of a rectangular domain with width = 39 mm and height = 30 mm. The radiation boundary condition is imposed at the upper boundary to model the input source that is expressed by the Gaussian-modulated pulse. This study investigates ultrasonic acoustic lenses, so the center frequency $f0$ of an incident pulse is set to $1\u2009MHz$. Figure 4(b) shows the waveform of the Gaussian-modulated pulse when $t0=7.5\u2009\mu s$ and $bw=4\xd71012$, which exhibits the behavior of a short-duration pulse. A mesh is created with 130 × 200 = 26 000 isoparametric quadrilateral elements. The maximum element size is determined as $\lambda 1\u2009MHz/10$, which is small enough to ensure the accuracy of the model. Additionally, the surface impedance boundary condition with $Zs=\rho wcw$ is applied at all other edges to suppress the reflected wave, where $\rho w$ and $cw$ are the mass density and the sound speed of water, respectively.

The design domain of a planar acoustic lens is 1.5 mm away from the input source. Also, the target position at which acoustic energy can be focused is $dfocus$ away from the front end of the lens. The time range of interest is from 0 to $200\u2009\Delta t=20\u2009\mu s$, where $\Delta t=0.1\u2009\mu s$. The *x*-symmetry model is used. The material properties used in this problem are presented in Table I (Tran *et al.*, 2017).

. | Mass density, ρ (kg/m^{3})
. | Bulk modulus, κ (Pa)
. |
---|---|---|

Water (background) | 1000 | 2.25 × 10^{9} |

Stycast 1090SI | 560 | 3.67 × 10^{9} |

F3301 | 1730 | 1.80 × 10^{9} |

. | Mass density, ρ (kg/m^{3})
. | Bulk modulus, κ (Pa)
. |
---|---|---|

Water (background) | 1000 | 2.25 × 10^{9} |

Stycast 1090SI | 560 | 3.67 × 10^{9} |

F3301 | 1730 | 1.80 × 10^{9} |

Based on the problem definition described above, we state the TO formulation as Eq. (46),

where $Tend$ is the $20\u2009\mu s$, $\nu e$ is the volume of element *e*, $V0$ is the volume of the entire design domain, and the target global volume fraction factor (*VFF*) is fixed at 0.6. The number of initial holes is set to $7\xd77=49$.

### B. Sensitivity validation

To validate the calculated sensitivities, a test example is considered in Fig. 5. Here, the objective function is to maximize the squared pressure at a specific focal point over the entire time range from 0 to $Tend$, as expressed in Eq. (46). The focal point where the objective function is calculated is indicated by a black cross symbol in Fig. 5, and the design domain has two holes as shown. The sensitivities are evaluated at the boundary points discretized by zero level-set represented implicitly along the boundary of these holes.

We first performed adjoint sensitivity analysis for the objective function and compared with the finite difference method (FDM)-based sensitivities calculated by perturbating the boundary points by $\Delta \Omega $ in the normal direction. For this comparison, one of the discretized boundary points is randomly selected. The sensitivity comparison is performed for two parameters, the number of time steps used in the transient primal finite element analysis (FEA) and the magnitude ($\Delta \Omega $) of the perturbation in the normal direction of the boundary point. To investigate the effect of the number of time steps on the calculated sensitivities, the magnitude of time step ($\Delta t$) is fixed to 0.1 μs, the value used in FEA, and only the number of time steps is adjusted.

Figure 6(a) shows the sensitivity comparison result according to the number of time steps. It is observed that the difference between the FDM-based sensitivity and the adjoint sensitivity decreases linearly as the number of time steps increases. Then Fig. 6(b) shows that as the perturbation magnitude ($\Delta \Omega $) decreases, the difference between the two decreases as well and then starts to increase again from a certain point, i.e., $\Delta \Omega \u224310\u22128$. This is due to the inherent truncation error of the sensitivity computed based on FDM. Thus, the AVM is used to efficiently calculate accurate sensitivity without such truncation error.

On the other hand, the adjoint sensitivity analysis method we applied in this study is also called the semi-discrete sensitivity analysis (Dilgen and Aage, 2020), because it only considers the spatial discretization, not both the spatial and time discretization. In recent years, it was shown that for transient optimization, fully discrete sensitivity analysis in a general matrix form can give more consistent sensitivities rather than the commonly used semi-discrete sensitivity analysis (Dilgen and Aage, 2020). The fully discrete sensitivity analysis method could be adopted in the future to obtain consistent sensitivities in the transient optimization of acoustic lens.

### C. Effect of localizing time-window function

This optimization study investigates the effect of the localizing time-window function by performing the four different optimizations with varying center times, $tc$. The lens size is specified to be $hlens=7.5\u2009mm$, and the target focal length $dfocus=3\u2009mm$.

Figure 7 presents the converged optimization solutions, and design material 1, 1090SI, and design material 2, F3301, are indicated with green and blue colors, respectively. Here, the convergence property of these optimization solutions can be confirmed by investigating the first-order necessary condition for optimality, as explained in Appendix A. Figure 7(a) is the optimum design when the time-window function is not considered (i.e., $\Delta tc$ is infinity), Fig. 7(b) when $tc=13.5\u2009\mu s$ using Eq. (31), and in Figs. 7(c) and 7(d) the $tc$ values are arbitrarily selected, $tc=9.5\u2009\mu s$ and $tc=17.5\u2009\mu s$, respectively.

It can be observed from the normalized intensity fields in the bottom row that the optimization focuses the acoustic energy at the target focal position for all solutions. The middle row figures depict the normalized pressure fields, and many scattering waves can be observed around the focal position for Figs. 7(a), 7(c), and 7(d). This demonstrates that calculating the center time $tc$ using Eq. (31) is an effective choice for the initial conceptual design of acoustic lenses. Since the optimal lens layout in Fig. 7(c) is obtained with the smallest center time of $tc=9.5\u2009\mu s$, the interface angle between the design materials placed in the side is the smallest to achieve focusing within the shortest time.

The transient responses of the pressure measured at the target focal position make the significance of the localizing time-window function clearer. Figure 8 shows the optimal lens layouts and their corresponding transient pressure responses. When the time-window function is not considered, Fig. 8(a) shows that the unwanted waves are additionally generated after the main focusing time. These waves make a clear separation between the incident and echo waves difficult, degrading the overall imaging performance. In contrast, Figs. 8(b)–8(d) show that acoustic energy is mainly focused in the time ranges specified through the time-window function, and the unwanted waves are suppressed. Moreover, it is possible to easily adjust the time range in which focusing performance occurs, by specifying the center time of the time-window function, $tc$. Figure 8 thus confirms that the localizing time-window function is effective in designing the acoustic lens.

### D. Effect of volume constraint

This section investigates the effect of volume constraint on the optimization results. Three optimizations are performed with and without the volume constraint. The center time $tc$ of the localizing time-window function is set to $tc=13.5\u2009\mu s$ for all three optimizations. The lens size is specified to be $hlens=7.5\u2009mm$, and the target focal length $dfocus=3\u2009mm$.

The optimization solutions are shown in Fig. 9. As illustrated in Fig. 9(c), the lens optimized without the volume constraint uses a greater amount of material 1, 1090SI, to focus more acoustic energy. In addition, in the case of *VFF* $=0.2$, Fig. 9(a), it can be seen that some of the incident acoustic energy cannot be focused and is discarded into the surrounding medium. This characteristic can be seen more clearly through the transient responses of pressure in Fig. 10.

On the other hand, for dynamic problems including acoustic lens optimization, the volume constraint may not play an important role in the optimization. This is because the optimal direction does not need to fully fill the design domain with material, as confirmed from the optimization results here. For following numerical examples, we will set $VFF$ to 0.6 to quantitatively compare optimization results.

### E. Effect of incident pulse type

This section elucidates how a short pulse-based acoustic lens differs from the more common lenses designed for a single operating frequency in terms of the focusing performance. The single operating frequency can be modeled as a single-tone continuous pulse in the time domain. The lenses are optimized for a short-duration pulse and continuous pulse using the proposed LSTO method, with the design conditions of the lens size $hlens=7.5\u2009mm$ and the target focal length $dfocus=3\u2009mm$. To simulate a continuous input pulse as close as possible to the actual one in the time domain, the sinusoidal pulse is incident on the lens for 20 $\mu s,$ as shown in Fig. 12(b). This is about 20 times larger than one period at 1 MHz and is a long enough time range to have a near single-tone. Please note that this is not an ideal single-tone but has a sufficiently narrow frequency bandwidth, as shown in the FFT results in Appendix B.

Figure 11 shows the optimized solutions, and they are completely different for the different types of the incident pulses. For the short-duration pulse in Fig. 11(a), the interface between the two design materials on the sides of the lens domain is inclined more toward the target focal position, and this is seen to focus the acoustic energy effectively. In contrast, the continuous pulse lens is more effective in maximizing the constructive interference of waves over the entire time range of interest, as observed from the normalized pressure fields in the middle row. Moreover, the acoustic energy accumulated by the continuous pulse lens is much greater than that accumulated by the short pulse-based lens.

To demonstrate the need for the consideration of transient responses for a short-duration pulse lenses, we present a cross-check analysis by applying the different types of incident pulses to the two optimized designs. The following four cross-check models and their corresponding results are shown in Fig. 12:

Model 1: The optimized lens for short-duration pulse with short-duration pulse.

Model 2: The optimized lens for continuous pulse with short-duration pulse.

Model 3: The optimized lens for continuous pulse with continuous pulse.

Model 4: The optimized lens for short-duration pulse with continuous pulse.

Figure 12(a) presents the squared pressure responses when the short-duration pulse is applied to the two optimized designs. It is observed from model 2 that the acoustic energy is focused at an undesired time of $\u223c10\u2009\mu s$ even though the center time is set at $13.5\u2009\mu s$. The peak amplitude is also reduced by about 4 times compared to that of model 1. Figure 12(b) shows the continuous pulse results for both lenses. Interestingly, for model 4, the lens is designed based on the short pulse, but some acoustic energy is focused when applying the continuous pulse. In other words, its squared pressure response is greater than that of the baseline without the lens. This is because the lens of model 4 is optimized for the short-duration pulse, where center frequency is the same as that of the continuous pulse. However, most of the acoustic energy is not focused and is wasted due to the dispersive characteristics of the lens. In other words, this lens has little focusing capability for frequencies other than the center frequency, 1 MHz, of the continuous pulse. To quantify the focusing performances, $\u222b0Tendptobj2dt$ are calculated for the four cross-check models and presented in Table. II. It can be clearly seen from Table II that the transient responses should be considered in the optimization formulation to design the planar acoustic lens.

$\u222b0Tendp(t)obj2dt$ . | Incident pulse type: Short-duration pulse . | Incident pulse type: Continuous pulse . |
---|---|---|

Lens optimized for short-duration pulse | 5.1 | 95.1 |

Lens optimized for continuous pulse | 2.6 | 193.0 |

$\u222b0Tendp(t)obj2dt$ . | Incident pulse type: Short-duration pulse . | Incident pulse type: Continuous pulse . |
---|---|---|

Lens optimized for short-duration pulse | 5.1 | 95.1 |

Lens optimized for continuous pulse | 2.6 | 193.0 |

### F. Effect of the lens size, $hlens$

We investigate the effects of the different lens sizes by optimizing for the lenses with $hlens=1.5$, $4.5$, and $7.5\u2009mm$. Figure 13 illustrates the design solutions obtained with the focal length of $dfocus=12.0\u2009mm$. All optimal lens layouts have a structure that allows acoustic waves to be direct toward the target focal position so that focusing performance is well accomplished. However, comparing the acoustic energy accumulated by the optimized lenses, Fig. 13(c), with the largest lens size shows the most improved results, which means that most of the acoustic energy can be focused. This is not surprising, as the $hlens$ has the greatest design space available to manipulate the waves.

The effect of lens size becomes clearer from the transient pressure responses presented in Fig. 14. As the lens size increases, the time-bandwidth in which acoustic energy is focused is gradually narrowed, thereby improving the focusing performance. From a physical point of view, the key principle of focusing acoustic waves at the target focal point is summarized in two ways: (a) the phase delay (or time delay) of acoustic waves through the lens domain and (b) the constructive/destructive interference of the transmitted acoustic waves. Therefore, as the lens size increases, the degree of experiencing the phase delay of acoustic waves further increases, thereby improving the focusing performance due to the more enhanced interferences and narrow time-bandwidth. On the contrary, if the lens size is very small, the focusing performance is remarkably reduced because the degree of phase delay is very small as well. This is more evident by looking at the squared pressure response in Fig. 14(b). As the lens size increases, more acoustic energy is focused in the time range of interest. It can be concluded from these results that the lens size should be large enough compared to the operating wavelength to achieve a satisfactory focusing performance of both the large amplitude and the narrower time-bandwidth.

### G. Effect of the size of focal length, $dfocus$

This section investigates how the focal length affects the optimized lens designs. Figure 15 depicts the optimized solutions for the three different target focal lengths, $dfocus=3.0$, $7.5$, and $12.0\u2009mm$. The lens size $hlens$ is set to 7.5 mm. Note that the three different optimization results are obtained using the different center times $tc$ of localizing time-window function. In other words, as the target focal length gets longer, $tc$ gets larger.

Interestingly, in the case of Fig. 15(a) with the smallest focal length, the acoustic energy accumulated by the lens is the greatest. This can be explained with the numerical aperture (NA), which is the geometrical property of the lens, defined as $NA=sintan\u22121rin/dfocus$, where the half width of the input source $rin=19.5\u2009mm$ and the focal length $dfocus$. Since the achievable focusing beam width is defined by the Rayleigh focusing limit of $0.61\lambda 1\u2009MHz/NA$ theoretically, the lens with a greater NA exhibits a narrower beam width (Hyun *et al.*, 2018; Maznev and Wright, 2017). As a result, the smaller the $dfocus$, the narrower the beam width that is achieved from the Rayleigh focusing limit. Moreover, if there is no energy loss, the total amount of focused energy is the same for all optimized solutions, which is defined as the product of intensity and the beam width at the target focal position. Therefore, when the narrower beam width is achieved, the intensity at the focal position becomes greater, resulting in a clearer image.

Figure 16 presents the transient responses of pressure measured at the target focal position for the optimized solutions. As the focal length $dfocus$ increases, the peak amplitude decreases and the time range in which acoustic energy is mainly focused is delayed from $13.5$ to $17.3\u2009\mu s$.

### H. Effect of the materials used

In the previous numerical tests, two design materials of F3301 and 1090SI were used, and it was shown that a satisfactory focusing performance can be achieved through their optimal distribution and combination. In this section, we investigate the focusing performance when only a single material is used. The following three material combinations are considered: (a) 1090SI-water (*single design materials*); (b) F3301-water (*single design materials*); and (c) F3301–1090SI (*two design materials*). For all combinations, the lens size $hlens$, the target focal length $dfocus$, and the center time of the localizing time function $tc$ are set to 7.5, 12, and 13.5 $\mu s$, respectively.

Figure 17 shows the optimization results, and for all material combinations, acoustic energy is focused well around the target focal position. On the other hand, comparing the values of accumulated acoustic energy $\u222b0Tendptobj2dt$, it seems that the 1090SI-water solution in Fig. 17(a) exhibits the best performance as it has the greatest value. However, the normalized pressure fields in the middle row suggest that the 1090SI-water solution has more scattering waves around the focal position relative to the other solutions. This indicates that some of the acoustic energy is not focused in the time domain specified by $tc$.

To further investigate these observations, we measure the transient responses of pressure at the target focal position for the optimized designs. Figure 18(a) shows that for the 1090SI-water solution, the time-bandwidth in which the acoustic energy is mainly focused is distributed more widely than in the other two cases. Figure 18(b) shows that although the maximum peak amplitude of the squared pressure is the greatest for the 1090SI-water solution, the F3301-1090SI solution exhibits the narrowest time-bandwidth of the solutions. Since the main objective of this study is to design a planar acoustic lens used in actual imaging systems, it is beneficial to clearly separate the incident and echo waves; hence, a narrower time-bandwidth yields a better imaging performance.

### I. Analysis of optimized lens with coupled acoustic-structure model

We optimized the acoustic lenses in Fig. 17(c) by considering only the pure acoustic model, excluding the elastic characteristics of the materials. We now examine the optimum lens performances when considering the elastic material behavior also by applying a fully coupled acousto-mechanics model via COMSOL Multiphysics (COMSOL, 2014). The bulk modulus $\kappa $ is given as $3.67\xd7109$ and $1.80\xd7109\u2009Pa$ for Stycast 1090SI and F3301, respectively, in Table 1 (Tran *et al.*, 2017), and Poisson's ratio of a typical solid has a small range of 0.33–0.45. For the plain-strain condition, Young's modulus, *E*, is given by $E=2\kappa 1+\nu 1\u22122\nu $ (Reddy, 2007). We therefore conduct three analyses with $\nu $ = 0.33, $\nu \u2009$= 0.39, and $\nu \u2009$= 0.45 with the acoustic properties in Table 1.

Figure 19 compares the acoustics performance results of the fully coupled acousto-mechanics analyses and the pure acoustic analysis for the design in Fig. 17(c). All acoustic intensity fields are normalized by their maximum values to show the focal point clearly. Figures 19(a)–19(d) show that the smaller the Poisson's ratio, the farther away the focal point is from the lens. The analysis result for Poisson's ratio of 0.33 shows that the focal length is shifted by about 16% compared to the target focal length, i.e., 12 mm, of the pure acoustic model. However, for the coupled models with Poisson's ratio above 0.39, the difference reduces to 2%, approximately approaching the ideal target focal length. Moreover, it can be seen from Figs. 19(e) and 19(f) that as Poisson's ratio decreases, the acoustic pressure and intensity at the focal point are decreased. This can be attributed to the destructive interference of the longitudinal and transverse waves within the lens. The discrepancies between the pure acoustic model and the coupled acousto-mechanics model are primarily due to the shear effect of design materials. This causes a mode conversion of longitudinal and transverse waves at the interfaces between the lens and the surrounding medium, water (Lee *et al.*, 2017).

From the coupled acousto-mechanics numerical analyses above, the actual focal length may differ by about 2%–16% from that obtained by the acoustics only model depending on the Poisson's ratio of the design materials. For increased accuracy, a higher fidelity analysis and design may be employed in consideration of the coupled acousto-mechanics effect (Dilgen and Aage, 2020).

## VII. CONCLUSIONS

This paper introduced a transient LSTO method that is demonstrated to be effective in the design of a planar acoustic lens working with a short-duration pulse. It is shown that the proposed method is a direct way to consider a generalized time-dependent pulse in designing an acoustic lens. The proposed design method is verified through a series of numerical investigations for various design parameters. Particularly, to effectively design a lens working with a short-duration pulse, it is beneficial to consider the localizing time-window function so that unwanted waves can be suppressed. Also, to achieve a satisfactory focusing performance, the lens size in the wave propagation direction needs to be sufficiently large compared to the operating wavelength. Otherwise, most of the acoustic energy cannot be focused and is wasted. In addition, the use of two or more design materials can provide improved solutions in terms of the focusing performance due to the possibility of exploration of a larger design space, opening up the applicability of a multi-material TO approach.

## ACKNOWLEDGMENTS

We acknowledge the support of National Science Foundation (NSF) Civil, Mechanical and Manufacturing Innovation (CMMI) Grant No. 1762530.

### APPENDIX A: CHECK OF THE FIRST-ORDER NECESSARY CONDITION FOR OPTIMALITY

We examined the first-order necessary condition for optimality, the Karush–Kuhn–Tucker (KKT) condition, over the design iteration. According to the KKT condition, if $\Omega *$ is the local minimum for the objective function $J\Omega $ subject to the constraints $gm\Omega \u2009m=1,\u20092,\u2009\u2026,\u2009Ng$ and they are continuously differentiable in the neighborhood of $\Omega *$, then $\u2207L\Omega *=0$ is satisfied. In our LSTO method, $L$ is the Lagrangian function,

where $\lambda J$ and $\lambda g1$ are the Lagrange multipliers, which are determined by solving the linearized sub-optimization problem at every design iteration. From the KKT condition,

Therefore, we solved the optimization problem for $tc=13.5\u2009\mu s$ in Sec. VI C and presented the histories of $\u2225\u2207J\u2225$, $\u2225\u2207g1\u2225$, and $\u2225\u2207L\u2225$ during the design iterations, as shown in Fig. 20(b), where $\u2225\u2218\u2225$ is the Euclidean norm of $\u2218$. It can be seen from this evolution history plot that $\u2207J$ and $\u2225\u2207g1\u2225$ are gradually converged, and $\u2225\u2207L\u2225$ keeps a very small value close to 0 throughout the entire design iterations, securing the convergence of optimal solutions. For the sake of simplicity, in practice, we checked the convergence through the predefined maximum iteration number or the change in the objective function during the last five iterations [Fig. 20(a)].

### APPENDIX B: FREQUENCY RESPONSES OF THE PULSES OBTAINED USING THE FFT

The frequency responses for the two pulses used in Sec. VI E, the short-duration Gaussian-modulated pulse and the continuous pulse, are obtained using the FFT. As shown in Fig. 21, the short-duration pulse has a wide frequency bandwidth from the center frequency of 1 MHz, and the continuous pulse has a very narrow frequency bandwidth, i.e., near single-tone.