Adaptive shape optimization with NURBS designs and PHT-splines for solution approximation in time-harmonic acoustics

]Javier Videla$^a$
]Ahmed Mostafa Shaaban$^b$
]Elena Atroshchenko$^{a,}$ [Corresponding author<EMAIL_ADDRESS><EMAIL_ADDRESS>
[add1]School of Civil and Environmental Engineering, University of New South Wales, Sydney, Australia
[add2]Institute of Structural Mechanics, Bauhaus-Universität Weimar, Weimar, Germany

Geometry Independent Field approximaTion (GIFT) was proposed as a generalization of Isogeometric analysis (IGA), where different types of splines are used for the parameterization of the computational domain and approximation of the unknown solution. GIFT with Non-Uniform Rational B-Splines (NUBRS) for the geometry and PHT-splines for the solution approximation were successfully applied to problems of time-harmonic acoustics, where it was shown that in some cases, adaptive PHT-spline mesh yields highly accurate solutions at lower computational cost than methods with uniform refinement. Therefore, it is of interest to investigate performance of GIFT for shape optimization problems, where NURBS are used to model the boundary with their control points being the design variables and PHT-splines are used to approximate the solution adaptively to the boundary changes during the optimization process.

In this work we demonstrate the application of GIFT for 2D acoustic shape optimization problems and, using three benchmark examples, we show that the method yields accurate solutions with significant computational savings in terms of the number of degrees of freedom and computational time.

Geometry Independent Field approximaTion Helmholtz equation shape optimization NURBS PHT-splines

§ INTRODUCTION

Due to rapid civil and transport development, noise pollution has become an important public health concern, requiring efficient noise control solutions, such as design of structures (e.g. noise barriers) with optimal acoustic performance. Numerical design is based on the solution of the wave propagation problem, modeled by the Helmholtz equation.

Numerical methods for the Helmholtz equation encounter two major challenges: the so-called “pollution error" and treatment of unbounded domains. The pollution error results from the numerical dispersion error, which is related to the discrepancy between the exact and numerical wave number $(k)$. It has been proven theoretically that for 2D and 3D problems, the pollution error cannot be fully eliminated <cit.>. It is controlled by adapting the mesh size to the wave length, i.e. the element size is held a few times smaller than the wave length. Hence, the number of elements grows proportionally to the wave number, making computations for large $k$ unfeasible. However, pollution error is inversely proportional to the order of shape functions and hence can be significantly reduced by the use of higher-order approximations. This makes isogeometric analysis (IGA), with its natural feature of polynomial degree elevation or $p-$refinement, an attractive alternative for solution of the Helmholtz equation, as demonstrated in <cit.>.

Another source of error in modeling wave propagation in unbounded domains results from introducing a boundary truncation domain (usually a circle of radius $R$ in 2D and a sphere of radius $R$ in 3D), in which the Sommerfeld radiation condition is modeled by the Absorbing Boundary Condition (ABC). The accuracy of modeling the ABC increases with $R$, nonetheless, larger computational domain leads to a larger system and higher computational cost. Another approach to reduce domain truncation error is by using higher order ABCs. This requires higher order derivatives of the shape functions, which present a problem in linear FEM, but can be easily overcome by the use of higher order NURBS in IGA. A detailed study of various ABCs in the context of Isogeometric Collocation was performed in <cit.>.

Isogeometric analysis (IGA) and its subsequent variations have been a topic of active research since IGA was originally introduced in 2005 by Hughes et al. <cit.>. The main objective of IGA is to connect the Computer-Aided Design (CAD) model directly with numerical analysis. In most cases, Non-Uniform Rational B-Spline (NURBS) are able to preserve the original geometry replacing Lagrange polynomials in the FEM discretization. It has been successfully proven that, IGA could be implemented in several engineering applications, such as structural mechanics <cit.>, fluid-structure interaction <cit.>, fracture mechanics <cit.>, electromagnetics <cit.>, Helmholtz equation <cit.>, among others. Ref. <cit.> presents a complete review for the application of IGA in different engineering aspects.

The necessity to truncate the boundary presents a serious disadvantage of domain-type methods in comparison with the boundary-type methods, such as boundary element method (BEM) and its isogeometric variation, namely IGABEM. In BEM, the Sommerfeld radiation condition at infinity is automatically satisfied by the fundamental solutions and the mesh burden is reduced to the boundary discretization <cit.>. BEM and IGABEM have been widely applied to acoustic problems in exterior domains <cit.>. Moreover, IGABEM has been efficiently paired with optimization methods for acoustic shape optimization in different research works <cit.>.

One of the fundamental problems of the numerical formulas which utilize either NURBS or B-splines as basis functions (such as IGA) is that, in 2D and 3D applications, basis functions are generated as a tensor-product of 1D structures without the ability of local refinement. This makes them expensive in terms of computational resources. Different basis functions have been proposed to overcome this disadvantage, such as T-splines <cit.>, truncated hierarchical B-splines (THB-splines) <cit.> and Polynomial splines over Hierarchical T-meshes (PHT-Splines) <cit.>. All of the aforementioned basis functions have been efficiently applied to the problems of statics, dynamics, acoustics. See <cit.>.

Geometry Independent Field approximaTion (GIFT) is a numerical scheme proposed in <cit.> as an alternative to iso-parametric methods that decouples the geometry and field spaces. It consists in employing certain basis functions (for example, NURBS) to model the computational domain, whereas employing different basis functions to approximate the field. GIFT with NURBS to model the geometry and PHT-splines to model the unknown solution were applied to problems of linear elasticity and Laplace equation <cit.>, bending of cracked Kirchhoff-Love plates <cit.>, and time-harmonic acoustics <cit.>. It was demonstrated that, the adaptive local refinement of the solution can obtain optimal convergence rates in the cases of solutions of reduced continuity. Additionally, Jansari et al. <cit.> developed the research work of <cit.> by enriching the PHT-splines solution with the partition of unity property using plane waves, showing in some cases that, the enrichment provides an enhancement of several orders of magnitude in the overall solution error.

In the field of optimization, IGA formulations have been applied to many engineering applications, such as: structural shape optimization <cit.>, composite structural optimization <cit.>, acoustic shape optimization <cit.>, piezoelectric energy harvesters <cit.>, thermal meta-materials <cit.>, heat conduction problems <cit.> and fluids <cit.> (for a comprehensive review, the reader is referred to <cit.>). In addition to the better accuracy of IGA over FEM per degree of freedom because of the NURBS higher order and higher continuity, IGA addresses an important issue in FEM shape optimization associated with the boundary representation and the evolution of mesh following the changing boundary. In IGA, the set of control variables that parameterize the boundary is naturally chosen as design variables, providing a tight link between the design, analysis and optimization models. The same approach to shape optimization is applicable in GIFT.

The optimization methods are classified in two families: gradient-free and gradient-based methods. The family of gradient-free methods includes, for example, Particle Swarm Optimization (PSO) <cit.> implemented in shape optimization problems in <cit.>, Genetic Algorithm <cit.> and its optimization applications in <cit.>, among others. The main advantage of gradient-free methods is their ability to find global minimum without any sensitivity analysis, however, at a much higher computational cost than gradient-based methods. Gradient-based methods <cit.> have much higher convergence rate, however, depending on the initial guess, they may converge to a local minimum. The performance of gradient-based optimization methods can be further improved by providing exact gradients, obtained from the shape derivatives of the objective function, the weak form of the problem and constraints to perform the sensitivity analysis <cit.>. However, in many applications shape derivatives are difficult to obtain. In such case, gradients are calculated using finite-difference approximations. This can be done efficiently for small and medium-size problems. In this work, we use Sequential Quadratic Programming (SQP) algorithm, which is considered as one of the most effective methods for solving nonlinear constraint optimization problems <cit.>. In SPQ, the sequence of quadratic sub-problems is solved at each iteration to obtain the search direction. The algorithm is efficiently implemented within Matlab  function.

Combining shape optimization with adaptive refinement can be traced back to the 80's. For example, the works of Kikuchi et al. <cit.> and Canales et al. <cit.>, are devoted to study the shape optimization problems for linear elasticity using FEM and adaptive refinements. Both works show that it is possible to achieve a process of automatic mesh generation and shape optimization, while the main drawback is that the mesh is prone to distortion, which is a usual problem with FE formulations. In a more recent work, Mohite and Upadhyay <cit.> proposed an adaptive shape optimization framework for laminated composite plates, showing the advantage of adaptive refinement over uniform refinement.

As pointed out by the review study from Upadhyay et al. <cit.>, adaptive mesh refinement is not fully integrated in practical applications, since the mesh refinement requires accesses to the exact geometry, which means, an automatic communication with CAD model is needed. The later is a drawback that can be alleviated with our GIFT formulation.

In the literature related to adaptive optimization with splines, Chen et al. <cit.> proposed an adaptive shape optimization scheme using T-splines and finite cell method for structural shape optimization problems. Topology optimization using morphable components and Hierarchical B-splines (HB-splines) <cit.>, as well as Truncated Hierarchical B-splines (THB-Splines) has been introduced by Xie et al.<cit.> showing superior numerical performance over its uniform counterpart.

It is of particular interest the work of Gupta et al. <cit.> on adaptive topology optimization using GIFT. In their work, they used NURBS to model the geometry, and PHT-splines to approximate the unknown displacement field and the density function in the framework of Solid Isotropic Material with Penalization (SIMP) method. It was shown that adaptive GIFT can achieve up to 90% CPU time reduction in comparison with uniform meshes.

In this paper, an adaptive shape optimization scheme for Helmholtz problems is proposed. GIFT scheme is employed, using PHT-splines to discretize the sound field, while NURBS are utilized to model the boundary and the interior of the computational domain. The shape optimization is performed over the control points that define the NURBS boundary, while adaptive refinement is done over the PHT-spline basis functions. This allows us to reduce both the number of the degrees of freedom, as well as the computational time of the full optimization process. The adaptive optimization in this work is based on the recovery-based error estimator, previously proposed in <cit.> for the Helmholtz equation and PHT-splines.

The remainder of the paper is organized as follows. Section <ref> is assigned for the preliminaries, where the boundary value problem for the Helmholtz equation, NURBS, PHT-splines, GIFT formulation, and a generic shape optimization problem are introduced. Section <ref> in particular is devoted to the adaptive optimization. In Section <ref>, numerical results for three benchmark examples are discussed. Conclusions are drawn in Section <ref>.

§ PRELIMINARIES

§.§ Boundary Value Problem (BVP) for the Helmholtz equation

The boundary value problem for the Helmholtz equation, as explained in Figure <ref>, in domain $\Omega\in\mathbb{R}^{2}$ with boundary $\Gamma$ consists in finding the spacial component of the acoustic pressure, $u$, such that:

Δu + k^2 u    = 0          in   Ω

u    = g         on   Γ_D

∂u∂n    = ikh         on   Γ_N

∂u∂n + αu     = f          on   Γ_R

where $\Delta$ is the scalar Laplace operator, $\boldsymbol{n}$ is a unit normal vector on $\Gamma$, outward to $\Omega$, $k$ is the wave number, and $i$ is the imaginary unit. Dirichlet, Neumann and Robin boundary conditions are prescribed on parts of the boundary, $\Gamma_{D} ,\Gamma_{N}$ and $\Gamma_{R}$ respectively ($\Gamma = \Gamma_{D}\cup\Gamma_{N}\cup\Gamma_{R}$ and $\Gamma_D\cap\Gamma_N = \Gamma_D\cap\Gamma_R = \Gamma_R\cap\Gamma_N = \emptyset$), in which $g$, $h$, $\alpha$ and $f$ are prescribed functions. For exterior problems, $u$ refers to a scattered wave and the boundary conditions on the scatterer boundary are written in terms of the incoming wave $u^{\text{inc}}$ and its normal derivative, <cit.>, i.e.:

\begin{equation}
\begin{split}
u &= -u^{\text{inc}} \quad\text{on} \quad \Gamma_{D} \\
\dfrac{\partial u}{\partial\boldsymbol{n}} &= -\dfrac{\partial u^{\text{inc}}}{\partial\boldsymbol{n}} \quad \text{on} \quad \Gamma_{N}
\end{split}
\end{equation}

Helmholtz acoustic problem.

Additionally, the solution satisfies the Sommerfeld radiation condition prescribed at infinity to truncate all possible reflections of spurious acoustic waves from the far-field. This condition is formulated for 2D problems as follows <cit.>:

\begin{equation}
\lim\limits_{r\rightarrow\infty}  \sqrt{r} \left(\dfrac{\partial u}{\partial r} - i k u\right) = 0. \label{eq:Sommerfeld_condition}
\end{equation}
where $r$ is the distance from the origin. In domain-type methods, the Sommerfeld radiation condition is replaced by the Absorbing Boundary Condition (ABC) on the truncation boundary $\Sigma$ (usually given by a circle of radius $R$ in 2D) <cit.>. In this work we use Bayliss–Gunzburger–Turkel condition of order 1, referred as BGT1, formulated in polar coordinates $(r, \theta)$ as:

\begin{equation}
\dfrac{\partial u}{\partial r} + \mathcal{B}u = 0,
\end{equation}

\begin{equation}
\mathcal{B}u = \left(\dfrac{1}{2 R} - i k\right)u \label{b_1}
\end{equation}

Note, that due to higher continuity of splines, it is also possible to accommodate ABCs of higher order in the framework of IGA and GIFT.

The weak form of the boundary value problem is written as follows: Find $u^h \in \mathscr{U}^h \subseteq \mathscr{U} = \left\lbrace u \in \mathbb{H}^1 ({\Omega}), u=g \hspace{0.15cm} \textrm{on} \ {\Gamma}_D \right\rbrace$ such that $\forall v^h \in \mathscr{V}^h_0 \subseteq \mathscr{V}_0 = \left\lbrace v \in \mathbb{H}^1 ({\Omega}), v=0 \hspace{0.15cm} \textrm{on} \ {\Gamma}_D \right\rbrace$,

\begin{equation}
a(u^h,v^h) = l(v^h)
\label{eq:weak_form}
\end{equation}

\begin{equation}
a(u^h,v^h) = \int _{\Omega} \triangledown  u^h \triangledown  v^h d\Omega -k^{2} \int _{\Omega}   u^h  v^h d\Omega + \alpha \int _{\Gamma_{R}} u^h v^h d\Gamma + \int _{\Sigma} \mathcal{B}u^h v^h d\Gamma
\label{eq:weak_form_a}
\end{equation}

\begin{equation}
l(v^h) =  i k \int _{\Gamma_{N}}h v^h d\Gamma + \int _{\Gamma_{R}} f v^h d\Gamma
\label{eq:weak_form_l}
\end{equation}

§.§ Non-Uniform Rational B-Splines (NURBS)

A set of $n$ B-spline basis functions of degree $p$ is defined in the parametric space $\xi \in [0,1]$ on a knot vector $\Xi = \{\xi_1,\cdots, \xi_i, \cdots ,\xi_{n+p+1}\}$ (a non-decreasing sequence of knots $\xi_i$ with $\xi_1 = 0$ and $\xi_{n+p+1} = 1$). Then, the $i$-th B-spline basis function of degree $p$, $N_{i, p}(\xi)$ is defined recursively as <cit.>:

\begin{equation}  \label{eq:7}
N_{i,0}(\xi) = \Bigg\{
\begin{array}{lcl}
1, \;\;\;\;\;\;\;\;\;\; \xi_i\leq\xi<\xi_{i+1} \\
0, \;\;\;\;\;\;\;\;\;\; \mbox{ otherwise }
\end{array}
\end{equation}
and for $p\geq1$

\begin{equation}  \label{eq:8}
N_{i,p}(\xi)=\frac{\xi-\xi_i}{\xi_{i+p}-\xi_i} N_{i,p-1}(\xi)+\frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}} N_{i+1,p-1}(\xi)
\end{equation}

The NURBS basis function is written as:

\begin{equation}  \label{eq:9}
R_{i,p}(\xi) = \frac{N_{i,p}(\xi)w_i}{\sum\limits_{\hat{i}=1}^{n} N_{\hat{i},p}(\xi)w_{\hat{i}}}
\end{equation}
where $w_i$ is the $i$-th weight.

A NURBS curve is parameterized by a set of NURBS basis functions and control points $B_i$ as follows:

\begin{equation}  \label{eq:1D-NURBS-parameterization}
\begin{bmatrix}
x(\xi) \\
\end{bmatrix} =\sum\limits_{i=1}^{n} R{_{i,p}} (\xi)B{_i
} \end{equation}

Two-dimensional NURBS basis functions are obtained as a tensor-product of one-dimensional NURBS. For each direction, a knot vector is defined: $\Xi=\{\xi_1,\cdots, \xi_i, \cdots ,\xi_{n+p+1}\}$ and $\Phi =\{\eta_1,\cdots, \eta_i, \cdots ,\eta_{m+q+1}\}$ (where $m$ is the number of basis functions in direction $\eta$ ($\eta\in[0, 1]$) and $q$ is their degree). Then, the corresponding NURBS basis functions are written as

\begin{equation}  \label{eq:NURBS_2D}
R_{ij}(\xi,\eta)=\frac{N_{i,p}(\xi) N_{j,q}(\eta) w_{ij}} {\sum\limits_{\hat{i}=1}^{n} \sum\limits_{\hat{j}=1}^{m} N_{\hat{i},p}(\xi) N_{\hat{j},q}(\eta) w_{\hat{i} \hat{j}}}
\end{equation}
where $w_{ij}$ is the weight associated with the control point $P_{ij}$. The NURBS surface is parameterized as

\begin{equation}  \label{eq:2D-NURBS-parameterization}
\begin{bmatrix}
x(\xi, \eta) \\
y(\xi, \eta)
\end{bmatrix} = \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} R_{ij} (\xi,\eta)P_{ij} \end{equation}

§.§ PHT-splines

Polynomial Splines over hierarchical T-meshes (PHT-Splines) were proposed by Deng et al. <cit.> as a generalization of B-splines over hierarchical T-meshes in order to provide both local refinement and adaptability using a polynomial basis, which is able to parameterize the geometry. The following section briefly covers the formulation of PHT-splines.

Consider $\mathscr{T}$ as a hierarchical T-mesh, $b_{i}(\xi,\eta)$, $i=1,2,...,n$ as a set of Bézier splines (B-splines), $\xi$ and $\eta$ as two parametric coordinates defined on the space $\left[0,1 \right] \times \left[ 0,1\right]$, and $P_{i}$ as the control points. Then, the polynomial spline surface over $\mathscr{T}$ at level 0 is defined by:

\begin{equation}
S(\xi,\eta) = \sum_{i=1}^n b_{i}(\xi,\eta)P_{i}
\end{equation}

A spline basis function is represented locally as a linear combination of Bernstein polynomials. In particular, each B-Spline $b_{i}(\xi,\eta)$ is defined as

\begin{equation}
b_{i} (\xi,\eta) = \sum_{j=1}^{p+1}\sum_{k=1}^{p+1}C^{i}_{jk}\hat{B}_{j,k}(\xi,\eta)
\end{equation}

where $\hat{B}_{j,k}(\xi,\eta) = B_{j}(\xi)B_{k}(\eta)$ is a tensor product of Bernstein polynomials defined on the reference interval $\left[-1,1\right]$ :

\begin{equation}
B_{j}(\xi) = \frac{1}{2^{p}}\binom{p}{j-1}\left(1-\xi\right)^{p-j+1}\left(1+\xi\right)^{j-1}
\end{equation}

The procedure to compute Bézier coefficients $C^{i}_{jk}$ based on a particular NURBS curve is detailed in <cit.>. Then, based on a B-spline surface written over an initial regular T-mesh (level-$k$), we briefly explain how to compute the new coefficients given a refinement on a level-$k+1$:

* Construct a Bézier representation of the basis functions. This is done using Bernstein polynomials and Bézier coefficients (other authors refer them as Bézier ordinates).

* Once a level refinement is applied, the shape functions whose support was lying on refined element must be modified. Given the Bézier representation of the basis function, the process consists in truncating the basis function by zeroing out Bézier coefficients. This process is also explained in the introductory paper to PHT-splines <cit.> and in the work of Garau and Vásquez <cit.>.

* After the truncation of coefficients, the insertion of new basis functions is performed. For each new basis function, the process includes the computation of the corresponding Bézier ordinates and the assignment of a new global basis index. The new basis functions are standard (non-truncated) B-splines computed based on the local knot vector information. In the case of a removed T-junction, the local knot vectors are to be computed based on the neighbour connectivity of the elements stored in the quad-tree structure. Finally, the new basis indices corresponding to the new basis functions are determined.

It is worth noting that this work is based on the open-source package IGAPACK (<https://github.com/canitesc/IGAPack>) which contains all the aforementioned PHT-splines routines. (See also <cit.>)

§.§ Geometry Independent Field approximaTion (GIFT)

Geometry Independent Field approximaTion (GIFT) <cit.> can be seen as an extension of the concept of isoparametric elements in standard IGA to sub- and super-parametric elements. In other words, computational domain can be kept in its original, coarse CAD parameterization, given by NURBS, e.g. eq. (<ref>), while other type of splines can be chosen to approximate the unknown solution. In this work we use PHT-splines, as it was implemented before in <cit.>. The representations for the solution and the geometric quantities are substituted into the weak form (<ref>)-(<ref>)-(<ref>) and after the integration is carried out, the weak form is transformed into a system of linear algebraic equations for unknown control points.

§.§ Shape optimization problem.

The aim of shape optimization problems in acoustic is to improve the acoustic performance of the system by modifying the boundary of the structure. Improving acoustic performance is achieved by minimizing certain quantity of interest $J(u)$ (e.g. level of sound in the protected zone) subjected to some design constrains. Notation $J(u)$ implies that $J$ is function of $u$, i.e. the solution of the Helmholtz boundary value problem (<ref>a)-(<ref>d). One of the main features of spline-based boundary parameterization used in this work is that the boundary is described by a set of control points $\boldsymbol{x}$, which are used as design variables. Hence, the problem is written as follows:

\begin{equation}
\begin{split}
\text{Objective function:   }& \min\limits_{\boldsymbol{x}} J(u), \text{   subject to eq.(\ref{eqn:Boundary_conditions_DNR}a)-(\ref{eqn:Boundary_conditions_DNR}d)}\\
\text{Such that:   }& a_i \leq x_i \leq b_i, \\
\text{             }& m(\boldsymbol{x}) = 0, \\
\text{             }& n(\boldsymbol{x}) \leq 0, \\
\end{split}
\end{equation}

where each design variable $x_i$ is set to take values within interval $[a_i, b_i]$,  and functions $m(\boldsymbol{x})$ and $n(\boldsymbol{x})$ represent the non-linear constraints, such as design limit on the area or volume of the structure.

§.§ Recovery-based adaptive refinement and optimization

This section illustrates the adaptive GIFT shape optimization framework proposed in this study.

In a classical adaptive refinement scheme (Figure <ref>), there are five main steps: Solve, Estimate, Mark, Refine and Check. As the name suggests, in the first step, the BVP is solved by the selected numerical method. Then, the error metrics are computed (e.g. recovery-based error estimator used in this work). The “Mark" step consists in using the marking algorithm to select elements which need to be refined. Subsequently, in the “Refine" step, those selected elements are refined. After that, depending on the convergence criteria, one can either repeat the process or stop it.

Classical scheme for adaptive refinement.

In this work, as stated in the previous section, the BVP is solved using a GIFT formulation with PHT-splines and NURBS. The shape optimization is performed on the physical space, which is parameterized by NURBS basis functions, and therefore, does not affect the PHT-splines basis functions. This means that the optimization process employs a fixed PHT-splines mesh. Consequently, one can sequentially apply optimization and adaptive refinement in a step-by-step way, which means that the optimization is performed over the NURBS control points, without affecting the PHT-splines, and then, after the optimization process has finished, the adaptive refinement is done over the PHT-splines mesh, keeping the NURBS structure unchanged. This scheme is further explained in Algorithm <ref>. The algorithm is controlled by four parameters: $\varepsilon_0$ - initial error tolerance, $\varepsilon_{\text{loop}}$ - error tolerance in each iteration, $\varepsilon_{\text{sol}}$ - error tolerance for the optimization solution, $N_{\text{Max steps}}$ - maximum number of steps. The process is also explained in Figure <ref>.

Adaptive shape optimization pseudo-algorithm

Input parameters: $\boldsymbol{x}_{0}$, $N_{\text{Max steps}}$, $\varepsilon_{0}$, $\varepsilon_{\text{loop}}$, $\varepsilon_{\text{sol}}$
Step 1: create initial NURBS structure based on $\boldsymbol{x}_{0}$ and PHT-splines meshe
Step 2: perform adaptive refinement over the PHT-splines basis, driven by the recovery-based error indicator, controlled by $\varepsilon_{0}$
$i = 0$
$i = i + 1$
perform shape optimization over the NURBS structure, obtaining $\boldsymbol{x}_{i}$ - best solution at $i-$th step
generate new NURBS structure using $\boldsymbol{x}_{i}$
perform recovery-based adaptive refinement based on the optimal solution $\boldsymbol{x}_{i}$ until error estimator reaches the error tolerance $\varepsilon_{\text{loop}}$
$i > N_{\text{Max steps}}$
end while-loop
$\left \| \boldsymbol{x}_{i} - \boldsymbol{x}_{i-1} \right \| < \varepsilon_{\text{sol}}$
end while-loop
update solution and update PHT-splines mesh

Scheme of algorithm 1.

The recovery-based error estimator proposed in <cit.> for the Helmholtz equation is employed.
The “Dörfler strategy’’ is used to mark the elements that need to be refined <cit.>. The strategy consists in sorting the elements according to their error contribution to the total error and then refining the elements, whose error contribution is bigger than a certain percentage of the error estimator. In the subsequent numerical experiments, 50% is utilized for the marking algorithm.

§ NUMERICAL RESULTS

We present three numerical examples to demonstrate the accuracy of the adaptive GIFT shape optimization. We devise one benchmark problem  with a known solution and known local minimum. Next, we analyze the shape optimization problem of a horn. And lastly, we analyze the shape optimization problem of a sound barrier. Both problems have been extensively studied in the literature, mostly in the context of isogeometric boundary element methods (IGABEM) <cit.>. As pointed out previously, NURBS are employed to model the geometry, and PHT-splines of degree $p=3$ are employed to model the sound pressure.

§.§ Plane wave scattering problem by a cylinder in a truncated domain

In the first numerical example, the problem of a plane wave scattered by a cylinder of radius $a$ is considered. For numerical simulations, infinite domain is truncated by an artificial boundary: circle of radius $R$. This truncation boundary is denoted by $\Sigma$ and the domain enclosed by $\Gamma$ and $\Sigma$ is denoted by $\Omega$, as shown in Figure <ref>. The boundary value problem with BGT1 boundary condition consists in finding $u$ such that:

\begin{equation}
\begin{split}
\Delta u + k^2 u &= 0 \text{   in   }\Omega\\
\dfrac{\partial u}{\partial \boldsymbol{n}} &= - \dfrac{\partial u^{\text{inc}}}{\partial \boldsymbol{n}}  \text{   on   }\Gamma\\
\dfrac{\partial u}{\partial r} + \left(\dfrac{1}{2 R} - i k\right)u &= 0\text{   on   }\Sigma
\end{split}\label{BVP-cylinder-BGT1}
\end{equation}

Truncated domain for the cylinder problem.

The incoming wave $u^{\text{inc}}$ travels in the horizontal direction, i.e:

\begin{equation}
u^{\text{inc}} = e^{ikx}\label{u_inc_full}
\end{equation}

In order to derive the analytical solution, the incoming wave is expanded in polar coordinates into a series:

\begin{equation}
e^{ikr\cos\theta} = \sum\limits_{n = 0}^{\infty}\kappa_n i^n J_{n}(kr)\cos n\theta,\label{exp_series}
\end{equation}

\begin{equation}
\kappa_0 = 1, \quad \kappa_n = 2\text{    for } n \geq 1
\end{equation}

and $J_n(k r)$ is Bessel function of the first kind of order $n$.

The detailed derivation of the exact solution $u^{\text{exact}}$ of boundary value problem (<ref>) is given in <cit.>. Here only the main results are recalled.

The solution is sought in the form:

\begin{equation}
u^{\text{exact}} = \sum\limits_{n = 0}^{\infty}(a_n H_n^{(1)}(k r) + b_n H_n^{(2)}(k r))\cos n\theta,\label{u_general_form}
\end{equation}

where $H_n^{(1)}(k r)$ and $H_n^{(2)}(k r)$ are the Hankel functions of the first and second types, respectively. The coefficients $a_n$ and $b_n$ are obtained by substituting eq.(<ref>) to boundary conditions of (<ref>) on $r = a$ and $r = R$, and matching the coefficients of Hankel functions. This procedure leads to the following system of equations for unknown $a_n, b_n$:

\begin{equation}
\left(
\begin{matrix}
A_{11}       & A_{12} \\
A_{21}       & A_{22} \\
\end{matrix}\right)
\left(
\begin{matrix}
a_{n} \\
b_{n} \\
\end{matrix}\right)
\left(
\begin{matrix}
b \\
0 \\
\end{matrix}\right)\label{sys-BGT1}
\end{equation}

\begin{equation}
A_{11} = \dfrac{d}{d r}H^{(1)}_{n}(k r)|_{r = a}, \quad A_{12} = \dfrac{d}{d r}H^{(2)}_{n}(k r)|_{r = a}, \,\,\, b = -\kappa_n i^n \dfrac{d}{d r} J_{n}(k r)|_{r = a}\label{BGT_BC_Neumann}
\end{equation}

\begin{equation}
\begin{split}
A_{21} &= k H^{(1)}_{n - 1}(k R) - \left(i k + \dfrac{(2 n - 1)}{2 R}\right) H^{(1)}_{n}(k R),\\
A_{22} &= k H^{(2)}_{n - 1}(k R) - \left(i k + \dfrac{(2 n - 1)}{2 R}\right) H^{(2)}_{n}(k R)
\end{split}
\end{equation}

The objective function is considered as follows:

\begin{equation}
J(u) = -\int\limits_{\Omega}\left | u \right |^2 dx \label{obj_circ}
\end{equation}

In the numerical results, $R = 2$ and $k = 0.25\pi$ are used. Shape optimization results are presented below for three different cases:

[case 1]: In the first case, we restrict the design space to circular geometries. Hence, the shape of the inner circle can be described with one design variable: its radius. The objective function as a function of radius $a$ is plotted in Figure <ref>. The minimum of the objective function $\min J(u) = -1.9193$ is reached at $a = 1.2689$. In Figure <ref> both the analytical and numerical solutions are presented to validate the implementation.

Cylinder example: objective function $J(u)$ as a function of radius $a$.

[case 2]: In the second case, to describe the shape of the inner boundary, we use six control points, distributed uniformly in the circumferential direction in the first quadrant, parameterized with three design variables ($x_{1}, x_{2}, x_{3}$), as shown in Figure <ref>. Such set up implies that control points are allowed to move only in the radial direction. The shape is created in the first quadrant and then rotated about the center to create shapes in the second, third and forth quadrants. The weights of all control points on the design boundary are set to 1. The design variables vary in the interval: $0.5 \leq x_1, x_2, x_3 \leq 1.8$.  Note that such parameterization cannot represent a circle exactly.

Cylinder example: control points on the inner boundary given by three design variables $(x_1, x_2, x_3)$. In case 2, all weights are equal to 1. In case 3, weights of points $P_2$, $P_3$, $P_4$, $P_5$ are additional design variables $x_4, x_5, x_6$, $x_7$, respectively.

[case 3]: In this study case, we use the same set up as in case 2, but adding weights of the internal points ($P_2$, $P_3$, $P_4$, $P_5$) as design variables ($x_4, x_5, x_6, x_7$, respectively). We assume that $0.5 \leq x_1, x_2, x_3 \leq 1.8$ and $ 0.75 \leq x_4, x_5, x_6, x_7 \leq 2$.

Numerical results for all three cases are shown next. For case 1, in Figure <ref> an initial geometry (a circle with radius $a = 1.5$) and a final geometry (a circle with radius $a = 1.2689$) can be seen. The optimization converges to the exact solution. For cases 2 and 3, an initial geometry is set up as a highly distorted shape, shown in Figure <ref>a. In both design spaces, optimization procedure, starting from this initial state, converges to two local minima, different from the one in the design space of case 1, shown in Figure <ref> and Figure <ref>.

The convergence process is shown in Figure <ref>. Optimization algorithm converges to the only minimum in the design space of case 1 with $\min J(u) = -1.9193$. In a richer space of case 2, the algorithm finds a local minimum with $\min J(u) = -1.985$, and in the design space of case 3, the local minimum with even smaller value of the objective function, $\min J(u) = -2.4461$, is found. Note that the results in terms of the objective function values for uniform and adaptive refinements are quasi-identical. It is an expected result, since at each step, the solution is refined (either adaptively or uniformly) until the error tolerance is lower than $\varepsilon_{\text{loop}} = 10^{-3}$.

Initial geometry, a circle of radius $a = 1.5$

Final geometry, a circle of radius $a = 1.2689$ - uniform mesh

Final geometry, a circle of radius $a = 1.2689$ - adaptive mesh

Cylinder example, [case 1]: initial and final shapes.

Initial geometry, given by $x_1 = 1$, $x_2 = 1.7$, $x_3 = 0.5$

Final geometry, given by  $x_1 = 1.1554$, $x_2 = 1.3788$, $x_3 = 1.3788$. - uniform mesh

Final geometry, given by  $x_1 = 1.1554$, $x_2 = 1.3788$, $x_3 = 1.3788$. - adaptive mesh

Cylinder example, [case 2]: initial and final shapes.

Initial geometry, given by $x_1 = 1$, $x_2 = 1.7$, $x_3 = 0.5$, $x_4 = 1$, $x_5 = 1$, $x_6 = 1$, $x_7 = 1$

Final geometry, given by  $x_1 = 1.0042$, $x_2 = 1.2570$, $x_3 = 1.2570$, $x_4 = 1.8534$, $x_5 = 1.8248$, $x_6 = 1.8248$, $x_7 = 1.8534$ - uniform mesh

Final geometry, given by $x_1 = 1.0042$, $x_2 = 1.2570$, $x_3 = 1.2570$, $x_4 = 1.8534$, $x_5 = 1.8248$, $x_6 = 1.8248$, $x_7 = 1.8534$ - adaptive mesh

Cylinder example, [case 3]: initial and final shapes.

Cylinder example: convergence of the objective function $J(u)$ during the optimization process.

Next, the adaptive shape optimization algorithm <ref> is studied. First, we compare the convergence plots of the recovery-based error estimator for the uniform and adaptive refinement for the initial shape in three cases in Figure <ref>. As it can be seen from the Figure, the solution for the initial shape in case 2(3) benefits slightly from the adaptive refinement due to the fact that, the shape is geometrically more complex than a circle in case 1. However, since this shape does not have any corners, the convergence rate for both, uniform and adaptive refinement, is optimal $(\sim h^{p/2})$.

Cylinder example: convergence of the $H_{1}-$error estimator for the initial shape in terms of DOFs.

In Figure <ref>, comparison in terms of the convergence of the objective function is shown for the initial shapes in case 1 and case 2(3). It can be seen that the results for uniform and adaptive refinements are quasi-identical for any number of DOFs. This can be explained by the fact that, in this particular example with simple regular geometries, objective function (in the form of an integral) converges very fast.

Cylinder example: convergence of the objective function $J(u)$ in terms of DOFs.

Then, uniform and adaptive shape optimization results are compared. Tables <ref>, <ref> and <ref> show the results of those comparisons in terms of DOFs, computational time and the optimal objective function obtained. In the three cases, the initial tolerance employed for the error estimator is $\varepsilon_0 = 10^{-2}$, while the loop tolerances for the error estimator are taken $\varepsilon_{\text{loop}} = 10^{-3}$ and  $\varepsilon_{\text{loop}} = 10^{-4}$. The last row shows the relative difference between the corresponding quantities for the uniform and adaptive refinements, calculated as shown for the time difference below:

\begin{equation}
\text{Difference} = \dfrac{\text{time}_{\text{uniform}} - \text{time}_{\text{adaptive}}}{\text{time}_{\text{uniform}}}\times{100}\%
\end{equation}

Results for cases 2 and 3 show that for both values of $\varepsilon_{\text{loop}}$, adaptive optimization overcomes uniform optimization in terms of DOFs and computational time. The results for case 1, on the other hand, show that adaptive refinement can deliver the optimal result with less DOFs, but with no conclusive difference in terms of computational time. This result can be explained by the fact that the geometry in case 1 is always smooth and simple, therefore there is no significant advantage in using adaptive refinement in this case.

1c     4c$\varepsilon_{\text{loop}} = 10^{-3}$
1lDOFs     1l Time (s)     1l Objective function     1lBest solution
Uniform                  1224     9.568     -1.91936     1.26889

Adaptive                 992     7.673     -1.91936     1.26889
Difference (%)          18.95     19.81     0      0

1l         1l                           1l

1c     4c$\varepsilon_{\text{loop}} = 10^{-4}$
1lDOFs     1l Time (s)     1l Objective function     1lBest solution
Uniform                  4488     23.334     -1.91936     1.26889

Adaptive                 3344     36.356     -1.91936     1.26889
Difference (%)          25.49     -55.81     0     0

Cylinder example: uniform and recovery-based adaptive refinement results for [case 1] with 1 design variable. $\varepsilon_{\text{0}} = 10^{-2}$ and $N_{\text{Max steps}}=10$.

1c     4c$\varepsilon_{\text{loop}} = 10^{-3}$
1lDOFs     1l Time (s)     1l Objective function     1lBest solution
Uniform                  4488     122.787     -1.985      [1.1554, 1.3788, 1.3788]

Adaptive                 1536     44.566     -1.985     [1.1554, 1.3788, 1.3788]
Difference (%)          65.78     63.70     0     0

1l         1l                           1l

1c     4c$\varepsilon_{\text{loop}} = 10^{-4}$
1lDOFs     1l Time (s)     1lObjective function     1lBest solution
Uniform                  17160     188.196     -1.985      [1.1554, 1.3788, 1.3788]

Adaptive                 5304     85.902     -1.985     [1.1554, 1.3788, 1.3788]
Difference (%)          69.09     54.36     0     0

Cylinder example: uniform and recovery-based adaptive refinement results for [case 2] with 3 design variables. $\varepsilon_{\text{0}} = 10^{-2}$ and $N_{\text{Max steps}}=10$.

1c     4c$\varepsilon_{\text{loop}} = 10^{-3}$
1lDOFs     1l Time (s)     1l Objective function     1lBest solution
Uniform                  17160                        861.643                                       -2.4461     [1.0042, 1.2570, 1.2570, 1.8534, 1.8248, 1.8248, 1.8534]

Adaptive                 2976                         234.360                                       -2.4461      [1.0042, 1.2570, 1.2570, 1.8534, 1.8248, 1.8248, 1.8534]
Difference (%)          82.66                            72.80                                     0     0

1l         1l                           1l

1c     4c$\varepsilon_{\text{loop}} = 10^{-4}$
1lDOFs     1l Time (s)     1lObjective function    1lBest solution
Uniform                  67080                        2759.996                                      -2.4461     [1.0042, 1.2570, 1.2570, 1.8534, 1.8248, 1.8248, 1.8534]

Adaptive                 17160                        666.817                                       -2.4461     [1.0042, 1.2570, 1.2570, 1.8534, 1.8248, 1.8248, 1.8534]
Difference (%)          74.42                           75.84                                   0     0

Cylinder example: uniform and recovery-based adaptive refinement results for [case 3] with 7 design variables. $\varepsilon_{\text{0}} = 10^{-2}$ and $N_{\text{Max steps}}=10$.

§.§ Acoustic horn

The acoustic horn problem studied in this section aims to reduce the back reflection of the incoming wave transmission in order to provide an efficient wave distribution towards the far-field. The horn is modeled as a 2D planar symmetry problem. The horn problem domain $\Omega$ is shown in Figure <ref> where $\Gamma_{in}$, $\Gamma_{rigid}$ and $\Gamma_{0}$ refer to the in-going wave inlet, the horn rigid boundary and the symmetry line, respectively. In order to model infinite domain, a truncation artificial boundary $\Sigma$ is defined as a semi-circle with radius $R$.

The horn problem domain.

The boundary value problem for this problem can be written as follows:

\begin{equation}
\begin{split}
\Delta u + k^2u &= 0 \qquad \qquad \, \, \text{in} \; \Omega, \\
\frac{ \partial u}{ \partial \boldsymbol{n}} + \left(ik+\frac{1}{ 2R}\right)u &= 0     \qquad  \qquad \, \text{on} \; \Sigma, \\
\frac{ \partial u}{ \partial \boldsymbol{n}} + iku &= 2ikA_m     \qquad \text{on} \; \Gamma_{in},\\
\frac{ \partial u}{ \partial \boldsymbol{n}} &= 0     \qquad \qquad \, \text{on} \; \Gamma_{rigid},  \\
\frac{ \partial u}{ \partial \boldsymbol{n}} &= 0     \qquad \qquad \, \text{on} \;   \Gamma_{0}
\end{split}
\end{equation}

in which $A_m$ is the incident wave amplitude transferring from $\Gamma_{in}$ into the acoustic domain.

The wave reflection coefficient $\mathcal{R}$ to be minimized in this optimization problem is defined as $\mathcal{R} = B_m/A_m$ on $\Gamma_{in}$, where $B_m$ is the reflected wave amplitude computed as $B_m=|u_{in}-A_m|$ and $u_{in}$ is the predicted acoustic pressure averaged on $\Gamma_{in}$ as $u_{in}= \frac{1}{a} \int_{\Gamma_{in}} u d\Gamma$ <cit.>. No constraint is defined in this optimization problem. Therefore, the objective function can be written directly as follows:

\begin{equation}        \label{eq:obj_horn}
J(u) = \mathcal{R}
\end{equation}

The red line in Figure <ref> represents the part of the horn, which is optimized in order to achieve the minimum of $\mathcal{R}$, where the $x$ and $y$ coordinates of the control points are set as design variables.

The initial geometry of the horn is considered as follows: the horn thickness is $h=2.5 \mathrm{cm}$ and the dimensions $a$ and $b$ are set as $5$ and $30 \mathrm{cm}$, respectively. The first dimension $a$ is constant for a length of $L_1=50 \mathrm{cm}$, after that, a consequent flaring of $L_2=50 \mathrm{cm}$ length is defined. The incident wave amplitude $A_m$ is chosen as $1$ and the radius $R$ is taken as $1 \mathrm{m}$.

For the GIFT implementation, the computational domain is split into five patches, as shown in Figure <ref>. The optimization will be performed only on the location of the control points that define the horn wall on the patch number 2. Note that, if the geometry parameterization is not refined during the solution refinement process, then a change in the horn boundary parameterization will lead to the change in the parameterization of the entire patch 2. Hence, the stiffness matrix for patch 2 needs to be re-calculated every time when the design parameters change during the optimization process. This is an unnecessary computational burden which can be avoided by reducing the continuity of the geometry parameterization in the radial direction. Continuity of the geometry basis is reduced by inserting knot value at $0.5$ twice. Then all changes in the system matrix are localized in the yellow area of patch 2, as shown in Figure <ref>.

Next, we present shape optimization results and then analyze the performance of the adaptive algorithm.

Horn optimization problem: NURBS parameterization.

§.§.§ Shape optimization results

The reflection spectrum of the initial shape is computed using the proposed GIFT model and compared in Figure <ref> against previously published spectra obtained by IGABEM <cit.> and FEM <cit.> where excellent agreement can be seen between the three spectra.

Reflection spectra of the initial horn.

The results obtained from the optimization problem with GIFT are compared against those obtained from both IGABEM with PSO and FEM with Genetic Algorithm reported in <cit.> and <cit.>, respectively. The optimization models of IGABEM in <cit.> use two design variables (the $x$ and $y$ coordinates of one control point) getting benefit from the automatic smoothness of IGA, while FEM models in <cit.> use only one design variable with applying additional curve smoothness using Hermite polynomials. Figure <ref> shows the optimized shapes for all comparative methods for different single frequencies $f = 280, 550, 780, 1000 \mathrm{Hz}$ with sound speed $c = 345\mathrm{m/s}$ and wave-number $k = 5.09940, 10.01667, 14.20546, 18.21213$, respectively. Figure <ref> illustrates the variation of the reflection coefficient $\mathcal{R}$ through the optimization process, while the reflection spectra for initial and optimized horn shapes are included in Figure <ref>. Figure <ref> shows the domain sound pressure values of the initial and optimized horns for frequencies $f = 280$, $550$, $780$ and $1000$ $\mathrm{Hz}$.

Optimized horn shapes using 1 control point - 2 design variables.

Convergence of the objective function $J(u) = \mathcal{R}$ in terms of number of iterations for the horn optimization problem with 1 control point - 2 design variables.

Reflection spectra for the optimized horn shapes using 1 control point - 2 design variables.

Initial sound pressure for $f=280$ Hz

Optimized sound pressure for $f=280$ Hz

Initial sound pressure for $f=550$ Hz

Optimized sound pressure for $f=550$ Hz

Initial sound pressure for $f=780$ Hz

Optimized sound pressure for $f=780$ Hz

Initial sound pressure for $f=1000$ Hz

Optimized sound pressure for $f=1000$ Hz

Sound pressure distribution (absolute value) for the initial and optimized horn shapes using 1 control point - 2 design variables.

In the next study, 2 control points (4 design variables) are studied in the optimization process. The sound speed is again $c = 345$ $\mathrm{m/s}$ and the horn is optimized for four different single frequencies $f = 280$, $400$, $550$ and $780$ $\mathrm{Hz}$. Figure <ref> shows the optimized shapes for all the comparative methods. Figure <ref> shows the convergence plot for each optimized shape, while Figure <ref> presents the reflection spectra plot for each shape. Overall, the results are in agreement with the references for the first three frequencies, on the other hand, for
$f = 780 \mathrm{Hz}$, the obtained results are quite different than the references. Still, Figure <ref> shows improvement in the obtained result versus the initial shape.

Optimized horn shapes using 2 control points - 4 design variables.

Convergence of the objective function $J(u) = \mathcal{R}$ in terms of number of iterations for the horn optimization problem with 2 control points - 4 design variables.

Reflection spectra for the optimized horn shapes using 2 control points - 4 design variables.

§.§.§ Performance of the adaptive shape optimization algorithm

To exploit the capabilities of PHT-splines and the GIFT formulation, an adaptive optimization algorithm, Algorithm <ref>, is studied in this section. First, Figures <ref> and <ref> show the comparison between uniform and recovery-based adaptive refinement for the non-optimized geometry for $f=1000$ $\mathrm{Hz}$. As it can be seen in the figures, adaptive refinement overcomes uniform refinement in both error estimator and reaching the asymptotic value of the objective function.

Horn problem: convergence of the $H_{1}-$error estimator in terms of DOFs.

Horn problem: convergence of the objective function $J(u) = \mathcal{R}$ in terms of DOFs.

Table <ref> shows the comparison between DOFs, computational time (s) and the optimal objective function for two different target tolerances, using $f=1000$ $\mathrm{Hz}$ and one control point (2 design variables). As it can be seen in the table, the adaptive optimization overcomes the uniform optimization in both, number of DOFs and computational time.

Figure <ref> shows the adaptive optimized mesh and the optimized sound pressure obtained for $f=1000$ Hz and 1 control point.

1c     4c$\varepsilon_{\text{loop}} = 10^{-2}$
1lDOFs     1l Time (s)     1lObjective function     1lBest solution
Uniform                  21516       171.225                     $4.449\times 10^{-3}$     [-0.0500, -0.0864]

Adaptive                 2556        36.074                     $4.953\times 10^{-3}$       [-0.0500, -0.0902]
Difference (%)          88.12                           63.49                                   -11.33     3.81
IGABEM <cit.>          -                          -                                   $3.330\times 10^{-3}$

1c     4c$\varepsilon_{\text{loop}} = 10^{-3}$
1lDOFs     1l Time (s)     1l Objective function     1lBest solution
Uniform                  331788     5768.166                 $4.485\times10^{-3}$     [-0.0500, -0.0882]

Adaptive                 11348       186.901                    $4.485 \times 10^{-3}$     [-0.0500, -0.0875]
Difference (%)          96.58                           94.76                                    0      0.69

IGABEM <cit.>          -                          -                                   $3.330\times 10^{-3}$

Horn example: uniform and recovery-based adaptive refinement results. $f$=1000 Hz and 2 design variables are employed.

Adaptive mesh

Optimized sound pressure (absolute value) for $f=1000$ Hz

Horn example: adaptive mesh for the optimized geometry and the domain values. $f$=1000 Hz and 2 design variables are employed.

Table <ref> shows the uniform and adaptive results for 2 control points (4 design variables) and $f=780$ Hz. Similar conclusions can be drawn from this study.

1c     3c$\varepsilon_{\text{loop}} = 10^{-3}$
DOFs        Time (s)      Objective function     Best solution
Uniform                  331788     524.238                  $4.1562\times 10^{-4}$                  [-0.050	-0.364	-0.108	-0.014]
Adaptive                 12484       424.371                    $4.1999\times 10^{-4}$     [-0.056	-0.341	-0.281	-0.001]
Difference (%)          96.24                           91.28                                    -0.968

IGABEM <cit.>          -                          -                                   $1.61\times 10^{-4}$

Horn example: uniform and recovery-based adaptive refinement results. $f$=780 Hz and 4 design variables are employed. $\varepsilon_{0} = 5\times 10^{-2}$

§.§ Sound barrier

In the next example, we consider the shape optimization problem for a sound barrier, as shown in Figure <ref>, where the design objective is to minimize the sound pressure in the protected area behind the barrier by changing the shape of the barrier front side (shown by the red line). The barrier occupies domain $[5, 5.2]\times[0, 3]$. A line source is located at $\textbf{x}_0 = (0, 1)$. The truncation domain is a circle of radius $R = 20$, centered at $(5.1, 0)$ with BGT-1 boundary condition prescribed on it. Figure <ref> shows the truncated domain with the barrier and the sound source. The front side of the barrier is parameterized by NURBS of degree $p=2$ and $N$ equally spaced control points. The end points, $(5,0)$ and $(5, 3)$ are fixed, and the shape is optimized by moving the inner control points along $x-$direction, i.e. by changing their coordinates $x_i$.

The vertical noise barrier problem.

Computational domain for the vertical noise barrier problem with the truncation boundary.

Sound barrier optimization problem: NURBS parameterization.

The full boundary value problem is formulated as follows:

\begin{equation}
\begin{split}
\Delta u + k^2 u &= \delta(\textbf{x} - \textbf{x}_0)
\, \text{   in   }\Omega\\
\dfrac{\partial u}{\partial \boldsymbol{n}} + \mathcal{B}u &= 0 \qquad \qquad \text{   on   }\Sigma\\
\dfrac{\partial u}{\partial \boldsymbol{n}} &= 0 \qquad \qquad \text{   on   }\Gamma_0 \\
\dfrac{\partial u}{\partial \boldsymbol{n}} &= 0 \qquad \qquad \text{   on   }\Gamma \\
\end{split}
\end{equation}

where $\Gamma_0 = \{(x, y): y = 0\}$ is the symmetry line. Operator $\mathcal{B}$ represents BGT-1 condition on $\Sigma$. In this example, instead of solving the Helmholtz equation with the source function given by Dirac delta function,
we decompose the solution into the Green's function (which satisfies Helmholtz equation with Dirac delta function and homogeneous Neumann boundary condition on $\Gamma_0$) and the unknown solution (which satisfies homogeneous Helmholtz equation in the domain and the modified boundary conditions).

Let us assume that

\begin{equation}
u = u_G + \hat{u}  \label{eq:full_u}
\end{equation}

such that

\begin{equation}
\begin{split}
\Delta u_G + k^2 u_G &= \delta(\textbf{x} - \textbf{x}_0)
\text{   in   }\Omega\\
\dfrac{\partial u_G}{\partial \boldsymbol{n}} &= 0 \text{   on   }\Gamma_0 \\
\end{split}
\end{equation}

Then $u_G = G(\textbf{x}, \textbf{x}_0)$ is Green's function in a half-plane, given by <cit.>

\begin{equation}
G(\textbf{x}, \textbf{x}_0) = -\dfrac{i}{4}H_0^{(1)}(k|\textbf{x} - \textbf{x}_0|)-\dfrac{i}{4}H_0^{(1)}(k|\textbf{x} - \bar{\textbf{x}}_0|), \label{Green_function}
\end{equation}

where $\bar{\textbf{x}}_0 = (x_0, -y_0)$. Note that, $G(\textbf{x}, \textbf{x}_0)$ satisfies the Sommerfeld radiation condition at infinity. The remaining part of the solution satisfies:

\begin{equation}\label{eq:barrier}
\begin{split}
\Delta \hat{u} + k^2 \hat{u} &= 0
\quad \quad \quad \, \, \text{   in   }\Omega\\
\dfrac{\partial \hat{u}}{\partial \boldsymbol{n}} &= -\dfrac{\partial u_G}{\partial \boldsymbol{n}} \quad \text{   on   }\Gamma \\
\dfrac{\partial \hat{u}}{\partial \boldsymbol{n}} &= 0 \quad \quad \quad \, \, \text{   on   }\Gamma_0 \\
\dfrac{\partial \hat{u}}{\partial \boldsymbol{n}} + \mathcal{B} \hat{u} &= 0 \quad \quad \quad \, \, \text{   on   }\Sigma\\
\end{split}
\end{equation}

The objective function is defined as

\begin{equation}
J(u) =  \bar{u} u,
\end{equation}

where sound pressure $u$ (and its conjugate $\bar{u}$) is calculated at a grid of $10\times10$ points in domain $[6, 16]\times[0.1, 3]$.
The full optimization problem is formulated as follows:

\begin{equation}
\begin{split}
\min\limits_{x_i} & J(u), \text{    subject to eq.}(\ref{eq:barrier}), (\ref{eq:full_u}), (\ref{Green_function})\\
\text{Such that: } A(u) &\leq 0.6, \\
x_i &\in [4.9, 5.1],
\end{split}
\end{equation}

where $A(u)$ is the area of the barrier.

For GIFT implementation, we split the computational domain into three patches, as shown in Figure <ref>. Similarly to the horn example, we reduce the continuity of the geometry parameterization in the radial direction to localize changes in the stiffness matrix. Continuity of the geometry basis is reduced by inserting knot value at $0.1$ twice. Then all changes in the system matrix are localized in the elements adjacent to the left side of the barrier. This area is shown in yellow color in Figure <ref>.

We start by conducting the convergence study for two values of frequency: $f = 400$ Hz and $f = 1000$ Hz. The results are shown in Figure <ref> together with the reference solutions from <cit.> and <cit.>, where we can see that the objective function converges to its asymptotic value for both, uniform and adaptive refinement of the solution mesh. Figure <ref> shows the convergence of the $H_{1}-$error estimator with respect to the DOFs, for the same two frequency values. Overall, convergence studies show that adaptive refinement reaches the reference objective function slightly faster than the uniform refinement, while the convergence of the error estimator is the same in both cases.

Sound barrier example: convergence of the objective function $J(u)$ in terms of DOFs.

Sound barrier example: convergence of the $H_{1}-$error estimator in terms of DOFs.

§.§.§ Shape optimization results

The optimization results are compared with the reference solutions from the literature. This problem was studied in <cit.> and <cit.> with the isogeometric boundary element method (IGABEM). The results are compared in Figure <ref> for four different numbers of control points (design variables): $N = 5, 10, 15, 20$.

Optimized barrier shapes for frequency 400 Hz and 5, 10, 15, and 20 design variables.

Convergence plots for the objective function and volume constraint are presented in Figures <ref> and <ref>, respectively.

Objective function $J(u)$ vs. number of iterations

Area $A(u)$ vs. number of iterations

Convergence process for the barrier problem.

The sound distribution for the optimized shapes is shown in Figure <ref>.

Initial sound pressure

optimized sound pressure for $N=5$

optimized sound pressure for $N=10$

optimized sound pressure for $N=15$

optimized sound pressure for $N=20$

Sound barrier example: sound pressure distribution (absolute value) for the initial and optimized barriers using 5, 10, 15 and 20 design variables.

§.§.§ Performance of the adaptive shape optimization algorithm

Uniform and adaptive optimization are tested by using Algorithm <ref> for design variables $N=5$, 10, 15 and 20. The studies are performed in terms of DOFs, computational time (s), optimal objective function $J(u)$ and optimal solution. All the results are summarized in Table <ref>. Results show agreement between the uniform and adaptive cases in terms of the objective function and the optimal solution. Regarding computational times, adaptive optimization overcomes uniform optimization for all the studied cases. It is interesting to notice that for cases $N = 15$ and 20, the differences in computational times are 18% and 90%, respectively. On the other hand, in terms of DOFs, uniform optimization seems to use less DOFs compared to adaptive optimization. This is explained by the results from the convergence study in Figure <ref>, since there is no noticeable advantage to use adaptive refinement versus uniform refinement when analyzing the error-estimator versus DOFs. The same is observed in the adaptive meshes which the optimization algorithm generates.

Figure <ref> shows the uniform and adaptive optimal shape for $N=10$, as well as the reference results from <cit.> and <cit.>. Again, uniform and adaptive optimization results as well as the reference results are in a good agreement. Uniform and adaptive final meshes are shown in Figure <ref>.

Optimized barrier shapes for frequency 400 Hz and 10 design variables. Comparison between uniform and adaptive refinement results.

Uniform mesh: optimized shape for $N=5$

Adaptive mesh: optimized shape for $N=5$

Uniform mesh: optimized shape for $N=10$

Adaptive mesh: optimized shape for $N=10$

Uniform mesh: optimized shape for $N=15$

Adaptive mesh: optimized shape for $N=15$

Uniform mesh: optimized shape for $N=20$

Adaptive mesh: optimized shape for $N=20$

Sound barrier example: optimized geometries, uniform and adaptive mesh, zoomed around the barrier area.

4c5 Design Variables
DOFs                            1cTime (s)     1cObjective function     1cBest Sol
1runiform      1r199176     1801.761                     2.818$\times{10}^{-3}$                        [0.000, -0.120, -0.200, 0.000, -0.063]

1radaptive     1r177976     1755.914                     2.814$\times{10}^{-3}$                         [0.000, -0.120, -0.200, 0.000, -0.064]

1rDifference (%)     1r10.64     2.54                     0.14                        -

1rIGABEM <cit.>      1r-     -                     2.801$\times{10}^{-3}$                         -
1l           1l         1l              1l
4c10 Design Variables
DOFs                            1cTime (s)     1cObjective function     1cBest Sol
1runiform      1r199176     4192.800                     1.886$\times{10}^{-3}$                         [0.000, 0.000, -0.023, -0.200, -0.200, 0.000, -0.168, 0.000, -0.057, -0.149]

1radaptive     1r211468     4189.357                     1.877$\times{10}^{-3}$                          [0.000, 0.000, -0.023, -0.200, -0.200, 0.000, -0.170, 0.000, -0.059, -0.150]

1rDifference (%)     1r-6.17     0.08                     0.48                        -

1rIGABEM <cit.>      1r-     -                    1.904$\times{10}^{-3}$                          -
1l           1l         1l              1l
4c15 Design Variables
DOFs                            1cTime (s)     1cObjective function     1cBest Sol
1runiform      1r199176     6098.439                     1.120$\times{10}^{-3}$                          [-0.200, 0.000, 0.000, -0.060, -0.141, -0.200, -0.200, -0.200, -0.020, -0.200, 0.000, 0.000, -0.062, -0.200	, 0.000]

1radaptive     1r232164     4987.424                     1.117$\times{10}^{-3}$                          [-0.200, 0.000, 0.000, -0.059, -0.141, -0.200, -0.200, -0.200, -0.021, -0.200, 0.000, 0.000, -0.062, -0.200, 0.000]

1rDifference (%)     1r-16.56     18.22                    0.27                        -

1rIGABEM <cit.>     1r-     -                    1.145$\times{10}^{-3}$                          -
1l           1l         1l              1l
4c20 Design Variables
DOFs                            1cTime (s)     1cObjective function     1cBest Sol
1runiform      1r199176           89967.211                               3.804$\times{10}^{-4}$                                          [-0.200, 0.000, 0.000, -0.152, 0.000, -0.200, -0.140, -0.185, -0.200, 0.000, -0.013, -0.200, 0.000, -0.200, 0.000, 0.000, -0.200, -0.113, -0.200, 0.000]

1radaptive     1r201564     9203.854                     3.664$\times{10}^{-4}$                        [-0.200, 0.000, 0.000, -0.151, 0.000, -0.200, -0.141, -0.185, -0.200, 0.000, -0.014, -0.200, 0.000, -0.200, 0.000, 0.000, -0.200, -0.112, -0.200, 0.000]

1rDifference (%)     1r-1.20     89.77                    3.75                        -

1rIGABEM <cit.>     1r-     -                    9.630$\times{10}^{-4}$                         -

Sound barrier example: uniform and recovery-based adaptive refinement results. Recovery-based error estimator tolerance is fixed to $\varepsilon_0 = \varepsilon_{\text{loop}} = 5\times{10}^{-2}$.

§ CONCLUSIONS

In this work we show application of Geometry Independent Field approximaTion (GIFT) for problems of acoustic shape optimization. Computational domain, including the optimized boundary, is modeled by NURBS, while the solution is approximated by adaptive PHT-splines. Accuracy of the solution at each step of the optimization process is controlled by the prescribed error tolerance. In three benchmark examples, the obtained results are shown to be in a good agreement with the reference solutions and it is shown that the adaptive algorithm can bring up to 96% reduction in terms of the number of degrees of freedom and 94% in terms of the computational time in comparison with uniform mesh refinement.

§ ACKNOWLEDGMENT

This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. Also The authors acknowledge support from the UNSW Resource Allocation Scheme managed by Research Technology Services at UNSW Sydney.