Authors: Sofie Plöckl Julian Held Tolga Gökalp Laura Winter

1. Introduction

The Wave Based Method (WBM) is an approach for solving acoustic and vibration problems that emerged as an alternative to the Finite Element Method (FEM) in the late 1990s. While FEM divides the problem domain into nodes and uses shape functions to interpolate between them, with solution accuracy depending on mesh refinement, WBM takes a fundamentally different approach. It employs wave functions that inherently satisfy the underlying differential equations throughout the domain, only introducing errors at the boundaries. These boundary errors are then minimized through weighted residual or least squares approaches, with the weighting factors serving as the degrees of freedom to be solved for. This means WBM's accuracy isn't mesh-dependent leading to significantly smaller system matrices and better convergence properties. This makes WBM particularly advantageous for mid-frequency problems where FEM becomes computationally prohibitive due to the fine mesh requirements needed to control pollution errors. The method's development progressed from simple 2D acoustic problems to handling complex geometries through domain decomposition techniques, coupling with FEM, and extension to 3D and vibro-acoustic problems. Key challenges include handling non-convex domains and optimizing computational efficiency.

This report outlines the enhancement of an existing MATLAB designed for single convex domain problems. The work focuses on implementing two key extensions: domain partitioning into convex subdomains and a multi-level approach as depicted in the Figure above. These developments enable the analysis of more complex geometries, including domains with internal holes. Through these implementations, the code's capabilities have been expanded to handle a broader range of acoustic problems.

2. The Wave Based Method

2.1. Basic Concept

The WBM is an indirect Trefftz approach with uses weighted wave functions to solve a boundary value problem. Those wave functions directly satisfy the underlying differential equations, while violating the boundary conditions. By using a weighted residual scheme the errors are forced to become zero, which is captured in the weighting factors.

The steady-state pressure response inside a domain $\begin{array}{l}\Omega\end{array}$ can be captured by the Helmholtz equation:

(1)	$\begin{array}{l}\displaystyle (\nabla + k^2)p(\mathbf{r}) = -j\rho\omega q\delta(R)\end{array}$

$\begin{array}{l}\displaystyle \forall \mathbf{r} \in \Omega\end{array}$

The chosen set of wave functions $\begin{array}{l}\Phi_a\end{array}$ are homogenous solutions of the governing Helmholtz equation (1).

(2)	$\begin{array}{l}\displaystyle \Phi_a(x,y)=e^{-j(k_{a,x}x+k_{a,y}y)\end{array}$

with $\begin{array}{l}k^2_{a,x}+k^2_{a,y}=k^2\end{array}$

The systems response is already well captured by the shape functions, as shown in the following picture:

The particular solution can be described as a free field pressure response:

(3)	$\begin{array}{l}\displaystyle \hat{p}(\mathbf{r})=\frac{\rho\omega q}{4}H^{(2)}_0(kR)\end{array}$

with the zero-order Hankel function of the second kind $\begin{array}{l}H^{(2)}_0\end{array}$ .

By adding the particular solution (3) to the $\begin{array}{l}n_a\end{array}$ wave functions (2) the total pressure field can be approximated by:

(4)	$\begin{array}{l}\displaystyle \hat{p}(\mathbf{r})= \sum_{a=1}^{n_a}\Phi_a(\mathbf{r})p_a+\hat{p}_q(\mathbf{r})\end{array}$

This pressure approximation directly satisfies the underlying differential equations and therefor introduces errors at the boundaries. Different boundary conditions can be modeled with the WBM:

(5)	$\begin{array}{l}\begin{matrix} L_v(p) &= &\frac{p}{\bar{Z}}\\ L_v(p) &= &\bar{v}_n\\ p &= &\bar{p} \end{matrix}\end{array}$

at $\begin{array}{l}\Gamma_Z\end{array}$

at $\begin{array}{l}\Gamma_v\end{array}$

at $\begin{array}{l}\Gamma_p\end{array}$

2.2. Weighted Residuals

In order to apply a weighted residual formulation regarding the boundary conditions, some approximation errors are induced, yielding the following residual error functions:

$\begin{array}{l}\begin{matrix} R_Z &= &L_v(\hat{p})-\frac{\hat{p}}{\bar{Z}}\\ R_v &= &L_v(\hat{p}) - \bar{v}_n\\ R_p &= &\hat{p} - \bar{p} \end{matrix}\end{array}$

at $\begin{array}{l}\Gamma_Z\end{array}$

at $\begin{array}{l}\Gamma_v\end{array}$

at $\begin{array}{l}\Gamma_p\end{array}$

Now we are using the Galerkin approach to force the residual functions to become zero by applying the same shape functions as before as weighting functions $\begin{array}{l}\tilde{p}\end{array}$ :

$\begin{array}{l}\displaystyle \int_{\Gamma_Z} \tilde{p}R_Z \,ds + \int_{\Gamma_v} \tilde{p}R_v \,ds - \int_{\Gamma_p} L_v(\tilde{p})R_p \,ds = 0\end{array}$

Since only finite sized models are amenable to numerical implementation and the boundary conditions are continuous, the boundary conditions can only be satisfied in an approximate way. Therefor we are using numerical Gauss integration along the boundaries.

When each wave function $\begin{array}{l}\Phi_a\end{array}$ is used as the test function $\begin{array}{l}\tilde{p}\end{array}$ , we get $\begin{array}{l}n_a\end{array}$ linear equations that make up the wave model:

$\begin{array}{l}\displaystyle \mathbf{A}^{(WR)}\mathbf{p}=\mathbf{b}^{(WR)}\end{array}$

With $\begin{array}{l}\mathbf{A}^{(WR)}\end{array}$ being the system matrix (also called stiffness matrix):

$\begin{array}{l}\displaystyle \mathbf{A}^{(WR)} = \int_{\Gamma_Z} \mathbf{\Phi} \left(L_v(\mathbf{\Phi}^T) - \frac{\mathbf{\Phi}^T}{Z}\right) \,ds + \int_{\Gamma_v} \mathbf{\Phi}L_v(\mathbf{\Phi}^T) \,ds - \int_{\Gamma_p} L_v(\mathbf{\Phi})\mathbf{\Phi}^T \,ds\end{array}$

Properties of $\begin{array}{l}\mathbf{A}^{(WR)}\end{array}$ :

- small compared to finite element models (shape functions already describe the solution quite well)
- fully populated
- complex and implicitly frequency dependent
- symmetric

and $\begin{array}{l}\mathbf{b}^{(WR)}\end{array}$ being the system vector:

$\begin{array}{l}\displaystyle \mathbf{b}^{(WR)} = - \int_{\Gamma_Z} \mathbf{\Phi} \left(L_v(\hat{p}_q) - \frac{\hat{p}_q}{Z}\right) \,ds + \int_{\Gamma_v} \mathbf{\Phi} (\bar{v}_n - L_v(\hat{p}_q)) \,ds - \int_{\Gamma_p} L_v(\mathbf{\Phi}) (\bar{p} - \hat{p}_q) \,ds\end{array}$

3. Subdomain-Partitioning-Approach

The convergence of the WBM is mathematically guaranteed when applied to convex domains. However, engineering applications frequently involve non-convex geometries, necessitating the integration of domain decomposition techniques into the WBM framework to ensure computational feasibility and solution accuracy. The following chapter explores the underlying equations in terms of coupling as well as the equivalent modifications of the existing WBM one domain code.

3.1. Coupling Conditions

Consider the original 2D acoustic problem, which is now split into two domains $\begin{array}{l}\Omega_\alpha\end{array}$ and $\begin{array}{l}\Omega_\beta\end{array}$ . The boundary $\begin{array}{l}\Gamma_b\end{array}$ of each subdomain $\begin{array}{l}\Omega_b\end{array}$ is the union of the boundary conditions and the interface $\begin{array}{l}\Gamma_{bi}\end{array}$ , at which the following conditions hold:

(6)	$\begin{array}{l}\begin{matrix} p_\alpha & = & p_\beta\\ L_{\beta v}(p_\beta) & = & -L_{\alpha v}(p_\alpha) \end{matrix}\end{array}$

with $\begin{array}{l}L_{b v}=\frac{j}{\rho\omega}\frac{d}{dn_b}\end{array}$

at $\begin{array}{l}\Gamma_{\alpha i}\end{array}$

at $\begin{array}{l}\Gamma_{\beta i}\end{array}$

Those Interface conditions (6) guarantee pressure and velocity continuity along the coupled edge through a direct approach, which means that no auxiliary interface variables are introduced. Additionally one continuity condition is imposed on each subdomain to ensure a well-posed problem. Another option to model the coupling is done by using the impedance. As this parameter links the velocity and the pressure, only one prerequisite needs to be considered.

After applying the weighted residual formulation the following coupling matrices $\begin{array}{l}\mathbf{C}\end{array}$

(7)

$\begin{array}{l}\begin{matrix} \mathbf{C_{\alpha\alpha}} & = & -\int_{\Gamma_{\alpha i}} L_{\alpha v}(\Phi_\alpha)\Phi_\alpha^T ds_\alpha\\ \mathbf{C_{\alpha\beta}} & = & \int_{\Gamma_{\alpha i}} L_{\alpha v}(\Phi_\alpha)\Phi_\beta^Tds_\alpha\\ \mathbf{C_{\beta\alpha}} & = & \int_{\Gamma_{\beta i}}\Phi_\beta L_{\alpha v}(\Phi_\alpha^T)ds_\beta\\ \mathbf{C_{\beta\beta}} & = & \int_{\Gamma_{\beta i}}\Phi_\beta L_{\beta v}(\Phi_\beta^T)ds_\beta \end{matrix}\end{array}$

as well as coupling vectors $\begin{array}{l}\mathbf{c}\end{array}$ can be evaluated:

(8)

$\begin{array}{l}\begin{matrix} \mathbf{c_{\alpha\alpha}} & = & \int_{\Gamma_{\alpha i}} L_{\alpha v}(\Phi_\alpha)\hat{p}_{\alpha q}ds_\alpha\\ \mathbf{c_{\alpha\beta}} & = & -\int_{\Gamma_{\alpha i}} L_{\alpha v}(\Phi_\alpha)\hat{p}_{\beta q}ds_\alpha\\ \mathbf{c_{\beta\alpha}} & = & -\int_{\Gamma_{\beta i}}\Phi_\beta L_{\alpha v}(\hat{p}_{\alpha q}) ds_\beta\\ \mathbf{c_{\beta\beta}} & = & -\int_{\Gamma_{\beta i}}\Phi_\beta L_{\beta v}(\hat{p}_{\beta q}) ds_\beta \end{matrix}\end{array}$

with $\begin{array}{l}s_b\end{array}$ being the boundary coordinate in the tangential direction of subdomain $\begin{array}{l}\Omega_b\end{array}$ .

3.2. Extended System Matrices

The wave model is constructed by ensuring that $\begin{array}{l}\hat{p}_b\end{array}$ satisfies both the boundary conditions (5) and the interface conditions (6) in an integral sense. This wave model formulation encompasses the coupling between the two acoustic subdomains:

$\begin{array}{l}\begin{bmatrix} \mathbf{A}_\alpha + \mathbf{C}_{\alpha\alpha} & \mathbf{C}_{\alpha\beta} \\ \mathbf{C}_{\beta\alpha} & \mathbf{A}_\beta + \mathbf{C}_{\beta\beta} \end{bmatrix} \begin{bmatrix} \mathbf{p}_\alpha \\ \mathbf{p}_\beta \end{bmatrix} = \begin{bmatrix} \mathbf{b}_\alpha + \mathbf{c}_{\alpha\alpha} + \mathbf{c}_{\alpha\beta} \\ \mathbf{b}_\beta + \mathbf{c}_{\beta\alpha} + \mathbf{c}_{\beta\beta} \end{bmatrix}\end{array}$

Properties of system matrix:

fully populated
complex and implicitly frequency dependent
non-symmetric, due to Coupling Matrices (7)

3.3. Modification of the Code

The presented flowchart depicts the basic algorithm of the modified code. To begin with the input data phase where necessary parameters (geometry, external loading, physical quantities) and model specifications (coupled edge) are initialized. The subsequent processing phase encompasses several sequential operations. At its core, the Wave Based Model creation process is executed through calculating the local stiffness matrix and local force vector for each subdomain. Following these computations, the procedure advances to assemble the global stiffness matrix and the global load vector. Note that in this step the calculation of the coupling conditions in global coordinates is performed. The algorithm then solves the resultant global wave model and concludes with the visualization of the computed results.

3.4. Validation

As the final part of this section, we are validating algorithmic modifications against the reference codebase. This process involves evaluating interface conditions through comprehensive examination of shared boundary properties to ensure consistent behavior across coupled domains. Additionally a convergence analysis using $\begin{array}{l}L_2-norm\end{array}$ calculations across various truncation parameters is performed, establishing solution stability and accuracy thresholds. Finally, frequency response characteristics to assess system behavior across the spectral domain are analyzed, providing insights into the dynamic performance of the system.

Continuity Validation

$\begin{array}{l}p_\alpha=p_\beta\end{array}$

This graph shows the continuity of absolute pressure across the shared interface.
Pressure should match between Subdomain 1 and Subdomain 2.
Any discontinuity may indicate issues with the coupling formulation.

$\begin{array}{l}L_{\beta v}(p_{\beta}) = -L_{\alpha v}(p_{\alpha})\end{array}$

A small mismatch at the boundary coordinates of the coupled edge can be observed, which arises from the performed numerical Gauss integration at the limits of the coupled edge.

$\begin{array}{l}L_{\beta v}(p_{\beta}) = -L_{\alpha v}(p_{\alpha})\end{array}$

Continuity ensures that energy transfer between subdomains is physically accurate.
Discontinuities could indicate incorrect normal derivative calculations.
Any discontinuity may indicate issues with the coupling formulation.

Convergence

$\begin{array}{l}L_2 = \sqrt{\sum_i |p_i|^2}\end{array}$

$\begin{array}{l}p_i\end{array}$ represents the computed pressure values at different points.
The sum is taken over all points in the domain.
Adding more wave functions does not significantly change the solution.
Our Code is limited to a specific truncation number, because of the applied Gauss integration with a maximum of 100 Gauss points. As a higher truncation number relates to higher wave oscillations, we would need a lot more Gauss points to capture high frequency behavior.

Frequency Response Function Analysis

A frequency response function verification is performed for 2 points, one in each subdomain. The result is compared to the reference code.
FRF point out natural frequencies of a system
In our case the response should match, since both solution are calculated with WBM

3.5. Conclusion

The presented Subdomain-Partitioning-Approach within the Wave Based Method (WBM) framework successfully incorporates coupling conditions to handle non-convex geometries, ensuring computational feasibility and solution accuracy. Our implemented code efficiently calculates local stiffness matrices and load vectors for each subdomain, assembles the global matrices, and solves the resultant wave model with visualization capabilities.

Future Improvements:

Implementing point load applications across different subdomains to assess load distribution and interaction effects.
Investigating geometries with holes to extend the applicability of the method to more complex engineering problems.
Expanding the code to handle multiple subdomains or highly complex shapes.
Enhancing the numerical integration by increasing the number of Gauss points, enabling higher truncation numbers and capturing high-frequency behaviors more accurately.

These advancements will enhance the robustness and flexibility of the WBM code, making it suitable for a broader range of acoustic and engineering applications.

4. Multi-Level-Approach

The Subdomain-Partitioning approach reaches its limits for complicated domains, where a more sophisticated domain partitioning needs to be performed. This usually leads to a strong increase in the number of necessary convex subdomains and thus the number of continuity conditions that need to be fulfilled between them. Moreover there are some problem settings in which no efficient partitioning can be achieved due to geometry, e.g. domains with circular obstructions where there is no sufficiently small number of convex subdomains that can be assigned at the boundaries. To be able to handle such tasks, the Multi-Level Approach was introduced, which allows to describe complex domain geometries and multiple obstacles in a cavity by a small number of subdomains.

In the Multi-Level-Approach a multi-scattering problem is divided into a set of single-scatterer problems, each of which are called "level". The general concept is visualized in the illustration below. In this case only one scatterer is considered, but one is able to extend the procedure analogously to a multiple scattering problems by adding one unbounded scatterer level for each considered obstruction. The problem domain is decomposed into one bounded level considering solely the outer boundary of the domain and neglecting the impact of the scatterer (left image) and one unbounded domain mimicking the influence of the obstacle (right image). For every obstruction inside the global domain a single artificial source is defined which is located in the center of the inclusion. With this a solution field for the artificial source in the unbounded domain is obtained. In the case of multiple scatterers the number of unbounded domains is going to increase by the amount of additional scatterers. Each of these decomposed domains (bounded, unbounded) are constructed independently defined by their own specific set of wave functions (see chapter 4.1). The goal is to assemble the global system of equations consisting of both bounded and unbounded contributions. Through solving this system the vector of unknowns is determined, which contains the weighting factors for the wave functions of all levels. The total solution field is the sum of the homogeneous solution and the particular solution. The homogeneous solution is the superposition of the weighted response fields of every level and is combined with the particular solution depending to the load of the system. Utilizing this methodology it is possible to describe the response of a complicated system with a strongly reduced number of unknowns.

4.1. Sets of wave functions

Following the classic Wave Based Method the shape functions for the bounded domain are defined as follows:

(9)	$\begin{array}{l}\displaystyle \mathbf{\Phi}_w^{}(\mathbf{r}(x,y)) = \mathbf{\Phi}_{w_r}(x,y) = \left\{ sin \left( k_{xw_r}x \right),cos\left(k_{xw_r}x) \right\} e^{-ik_{yw_r}y}\end{array}$

and

(10)	$\begin{array}{l}\displaystyle \mathbf{\Phi}_w^{}(\mathbf{r}(x,y)) = \mathbf{\Phi}_{w_s}(x,y) = e^{-ik_{xw_s}x} \left\{ sin \left( k_{yw_s}x \right),cos\left(k_{yw_s}x) \right\}\end{array}$

(11)	$\begin{array}{l}\displaystyle k_r = \left( \frac{r\pi}{L_x}, \pm \sqrt{k^2_j - k^2_{r,x}} \right), \forall r = 0,1,\ldots,n_r\end{array}$

(12)	$\begin{array}{l}\displaystyle k_s = \left( \pm \sqrt{k^2_j - k^2_{s,y}} ,\frac{s\pi}{L_y}\right), \forall s = 0,1,\ldots,n_s\end{array}$

The set of functions for the unbounded domain represents the solution of a single source in space and are defined as follows, where $\begin{array}{l}H_w^{(2)}\end{array}$ denotes the Hankelfunction of the second kind:

(13)	$\begin{array}{l}\displaystyle \mathbf{\Phi}_w^I(\mathbf{r}(r,\theta)) = \mathbf{\Phi}_{w_c}(r,\theta) = H_w^{(2)}(k_jr)cos(w\theta)\end{array}$

and

(14)	$\begin{array}{l}\displaystyle \mathbf{\Phi}_w^I(\mathbf{r}(r,\theta)) = \mathbf{\Phi}_{w_s}(r,\theta) = H_w^{(2)}(k_jr)sin(w\theta)\end{array}$

$\begin{array}{l}\displaystyle w= 0,1,2,\ldots,w_{max}\end{array}$

$\begin{array}{l}\displaystyle w = 1,2,3,\ldots,w_{max}\end{array}$

The descriptive wave functions of the bounded and unbounded domains not only differ in their construction but also in the corresponding reference frame. The bounded wave functions are defined in the cartesian coordinate system while the unbounded domain is set up in a polar coordinate system with its origin in the center of the scatterer. This has to be considered during calculation of the linking matrices and assembly of the global system. In this example the global system is defined in the cartesian space so a transformation of coordinates has to be performed for the functions related to the unbounded domain.

4.2. Global Wave Model

The global wave model is set up by defining the system matrices of each level $\begin{array}{l}\mathbf{A_i}\end{array}$ , the coupling matrices $\begin{array}{l}\mathbf{C}_{i,j}\end{array}$ between level i and j, the vector of unknowns $\begin{array}{l}\mathbf{p}_w\end{array}$ and the load vector $\begin{array}{l}\mathbf{F}_w\end{array}$ considering the external loading of the system.

(15)

$\begin{align} \begin{bmatrix} \mathbf{A}_1 & \mathbf{C}_{1,i} & \cdots & \mathbf{C}_{1,n_j} \\ \mathbf{C}_{i,1} & \mathbf{A}_i \\ \vdots & & \ddots & \vdots \\ \mathbf{C}_{n_j,1} & & \cdots & \mathbf{A}_{n_j} \end{bmatrix} \begin{bmatrix} \mathbf{p}_w^{(1)} \\ \mathbf{p}_w^{(i)} \\ \vdots \\ \mathbf{p}_w^{(n_j)} \end{bmatrix} = \begin{bmatrix} \mathbf{F}_1 \\ \mathbf{F}_i \\ \vdots \\ \mathbf{F}_{n_j} \end{bmatrix} \end{align}$

$\begin{array}{l}\mathbf{A}_i\end{array}$ :

$\begin{array}{l}\mathbf{p}_w^i\end{array}$ :

$\begin{array}{l}\mathbf{C}_{i,j}\end{array}$ :

$\begin{array}{l}\mathbf{F}_i\end{array}$ :

$\begin{array}{l}n_j\end{array}$ :

stiffness matrix level i

weighting factors level i

coupling matrix between levels i,j

force vector level i

number of scatterers

The construction rules for the stiffness matrix $\begin{array}{l}\mathbf{A}_i\end{array}$ for each level, the coupling matrix $\begin{array}{l}\mathbf{C}_{i,j}\end{array}$ between level i and j and the force vector $\begin{array}{l}\mathbf{F}_i\end{array}$ of each level are defined as:

(16)

$\begin{align} \mathbf{A}_{i} = \int_{\Gamma_v^{i}} \mathbf{\Phi}_i \left( \mathcal{L}_v \left( \mathbf{\Phi}_i^T \right) - \bar{v}_n \right) d \Gamma + \int_{\Gamma_Z^{i}} \mathbf{\Phi}_i \left( \mathcal{L}_v \left( \mathbf{\Phi}_i^T \right) - \frac{\mathbf{\Phi}_i^T}{\bar{Z}_n} \right) d \Gamma + \int_{\Gamma_p^{i}} -\mathcal{L} \left( \mathbf{\Phi}_i\right) \left( \mathbf{\Phi}_i^T - \bar{p} \right) d\Gamma \end{align}$

(17)

$\begin{align} \mathbf{C}_{i,j} = \int_{\Gamma_v^{i,j}} \mathbf{\Phi}_i \left( \mathcal{L}_v \left( \mathbf{\Phi}_j^T \right) - \bar{v}_n \right) d \Gamma + \int_{\Gamma_Z^{i,j}} \mathbf{\Phi}_i \left( \mathcal{L}_v \left( \mathbf{\Phi}_j^T \right) - \frac{\mathbf{\Phi}_j^T}{\bar{Z}_n} \right) d \Gamma + \int_{\Gamma_p^{i,j}} -\mathcal{L} \left( \mathbf{\Phi}_i\right) \left( \mathbf{\Phi}_j^T - \bar{p} \right) d\Gamma \end{align}$

(18)

$\begin{align} \mathbf{F}_i = \int_{\Gamma_v^i} \mathbf{\Phi}_i \left( \mathcal{L}_v(\hat{p}_q) - \bar{v}_n \right) d \Gamma + \int_{\Gamma_Z^i} \mathbf{\Phi}_i \left( \mathcal{L}_v (\hat{p}_q) - \frac{\hat{p}_q}{\bar{Z}} \right) ds + \int_{\Gamma_p^i} -\mathcal{L}_v (\mathbf{\Phi}_i) \left( \bar{p} - \hat{p}_q \right) d\Gamma \end{align}$

The weighted residuals in this formulas as well as the particular solution term are defined in chapter 2 above. Equation (17) shows that the linking of two layers is accomplished by evaluating the residuals of one level and weighting them with the weighting functions of the other.

This system of equations is solved for the vector of unknowns $\begin{array}{l}\mathbf{p}_w\end{array}$ which contains the weighting factors for the respective wave functions. The total solution field is obtained by superposition of the homogeneous solution, consisting of the weighted wave functions, and the particular solution:

(19)	$\begin{align} \mathbf{p} (\mathbf{r}) = \mathbf{p}_w^{(1)} \mathbf{\Phi}_1 (\mathbf{r}) + \sum_{i = 1}^{n_j}\mathbf{p}_w^{(i)}\mathbf{\Phi}_i(\mathbf{r}) + \hat{\mathbf{p}}_q(\mathbf{r}) \end{align}$

4.3. Procedure of Matlab code

The main task is the implementation of the Multi-Level-approach in Matlab by expanding the existing introductory code by Hannes Englert at the TUM Chair of Structural Mechanics with the previously discussed system quantities. In the following the general procedure the code is following can be found:

Define input data (domain and obstacle definition, boundary conditions etc.)
Creation of Gauss points along all boundaries
Creation of wave numbers (necessary for bounded set of wave functions)
Computation of stiffness matrices $\begin{array}{l}\mathbf{A}_i\end{array}$ of each level i
1. Evaluation of residuals of level i along boundaries
2. Weighting with the same set of functions i
Computation of linking matrices $\begin{array}{l}\mathbf{C}_{i,j}\end{array}$ of levels i,j
1. Evaluation of residuals of level j
2. Weighting with the other set of functions i
Assembly to global system matrix
Computation of load vector $\begin{array}{l}\mathbf{F}_w^{(i)\end{array}$ for each level
1. Evaluation of residuals of load functions along boundaries
2. Weighting with each set of wave functions
Assembly to global load vector
Solve system of equations for vector of unknowns (weighting factors) $\begin{array}{l}\mathbf{p}_w\end{array}$
Postprocessing
1. Evaluation of weighted wave functions of all levels in entire domain (homogeneous solution)
2. Evaluation of load function in entire domain (particular solution)
3. Obtain total solution field as superposition of homogeneous solution of all levels i and particular solution

Apart from obtaining the solution field using the Multi-Level-Approach of the Wave-Based-Method the script also tackles the task of depicting results and boundary values.

4.4. Validation

As a validation of the modified code the compliance of the prescribed boundary conditions of the pre-existing introductory code and the new code are compared. As the new code currently is limited to an imposition of a vanishing normal velocity at the boundaries, both codes were executed for this type of boundary condition using the same configuration of the unobstructed level. For the following graphs the bounded domain was discretised by $\begin{array}{l}n_1 = n_2 = 30\end{array}$ and the unbounded level by $\begin{array}{l}n_3 = 3\end{array}$ . A sketch of the implemented domain is shown in section 4.6 (profile A).

Comparing the velocity distribution on the outer boundaries of the introductory code with the Mulit-Level-Approach a strong similarity both in form and scaling is achieved. When looking at the inner boundaries the velocity deviates stronger from zero but shows similar scaling to the outer boundaries. From this observation a sufficient compliance with the boundary constraints is assumed. Nevertheless, as especially at the inner boundaries significant oscillations and deviations occur, the results of such simulations should always be assessed under precautions.

4.5. Convergence studies

The precision of this numerical simulation strongly depends on the number of functions used to describe the responsive behaviour of each level. In literature one can find truncation rules for the different sets of functions which state $\begin{array}{l}n_r \ge T \frac{k_jL_x}{\pi}\end{array}$ and $\begin{array}{l}n_s \ge T \frac{k_jL_y}{\pi}\end{array}$ for the bounded level and $\begin{array}{l}w_{max} \geq 2RTk_j = Tk_j\sqrt{L_{x,inclusion}^2 + L_{y,inclusion}^2}\end{array}$ for an unbounded level containing a rectangular inclusion respectively. With the truncation parameter $\begin{array}{l}T= 1\ldots 6\end{array}$ and a desired similarity in spatial resolution a range of suitable configurations can be defined. To examine the behaviour of the code for an increasing number of wave functions a series of exemplary configurations were executed for the specifications below.

	run 1	run 2	run 3
$\begin{array}{l}T_{bounded}\end{array}$	$\begin{array}{l}0.5\end{array}$	$\begin{array}{l}5\end{array}$	$\begin{array}{l}14\end{array}$
$\begin{array}{l}n_r\end{array}$	$\begin{array}{l}3\end{array}$	$\begin{array}{l}30\end{array}$	$\begin{array}{l}82\end{array}$
$\begin{array}{l}n_s\end{array}$	$\begin{array}{l}2\end{array}$	$\begin{array}{l}15\end{array}$	$\begin{array}{l}41\end{array}$
$\begin{array}{l}T_{unbounded}\end{array}$	$\begin{array}{l}0.3\end{array}$	$\begin{array}{l}1\end{array}$	$\begin{array}{l}3\end{array}$
$\begin{array}{l}w_{max}\end{array}$	$\begin{array}{l}1\end{array}$	$\begin{array}{l}3\end{array}$	$\begin{array}{l}8\end{array}$

Firstly the influence of the increasing resolution on the boundary conditions was investigated:

The above figures depict that the first refinement step yields significant improvement regarding the compliance if the normal velocity to the imposed constraints on the outer boundaries. Especially the more distant edges show pretty nice approximations. The differences between the second and the third run however show no relevant enhancement especially regarding the tremendously increased number of functions and thus degrees of freedom which directly relate to the computation effort.

A look onto the progression of the normal velocities at the boundaries of the inclusion leads to a similar conclusion. The first refinement step especially improves the scaling i.e. the amplitude of deviations. The second step shows an interesting effect, where contrary to the behaviour on the outer boundaries the quality improves significantly amplitude-wise. Furthermore the number of slightly deflecting oscillations rises.

From these observations no disadvantage for arbitrarily high numbers of wave functions is apparent. This changes taking into account the influence on the pressure field response.

Again, comparing run 1 and run 2 the quality of the presentation of the system response increases significantly. In the first run the influence of the inclusion is noticeably overrepresented while in the second run a well balanced pressure field is achieved. The third run then establishes the necessity of an upper limit for the truncation parameters. Very high numbers of wave functions introduce very high numerical oscillations into the system and yield implausible results as well as multiple singularities along the edges. These effects are also discernible in the depiction of the imaginary parts of the total pressure field.

4.6. Results

Next the results of the Multi-Level-Approach for various profiles were examined. Following the findings of the previous chapter the numeric configuration was set to $\begin{array}{l}T_{bounded}= 5\end{array}$ and $\begin{array}{l}T_{unbounded} = 1\end{array}$ and the following spatial structures:

The first quantity of interest is the real part of the total pressure response:

It is apparent that in profile A the inclusion does not have to much influence on the real part of the pressure response. However looking at profile B the vicinity of the inclusion to the source causes distinctly higher oscillations and amplitudes throughout the whole domain demonstrating the increasing impact of the inclusion on the whole field with decreasing distance to the source. In the analysis of profile C the area in front of the big inclusion displays relatively small amplitudes and little oscillations. At the boundaries of the inclusion however the model creates really high peaks and singularities which increase for an increasing length of the edge with the majority of the peaks appearing in the middle of the edge. Behind the inclusion unexpectedly high amplitudes occur which do not coincide with the anticipated physical outcome. This already hints at limitations and problems this particular numerical model might entail.

These observations become even clearer in the depictions of the imaginary part. Comparing profiles A and B the change in influence of the inclusion versus the loading is discernible. The decline in quality of the representation in profile C also becomes apparent comparing A and C as the area in front of the inclusion shows nearly identical behaviour while the inclusion boundaries itself and the respective shadow area display artificial numeric oscillations and singularities:

4.7. Conclusion and outlook

From the tests conducted with the modified Matlab code it can be stated that this methodology works nicely for certain conditions. In general, this model is able to simulate the scattering behaviour of a rectangular inclusion. To expand the applicability of the code some additional implementation possibilities are suggested:

Introduce further boundary conditions (impedance, pressure) into the construction of the system matrices to cover a broader variety of problems
Consideration of multiple scatterers with variable cross sections (e.g implementation of round scatterers since this represents a deficiency in the classical wave based method and is one of the advantages of the multi-level-approach)

Furthermore an inspection of the numerical impact of the inclusion and introduction of artificial deflections, especially for very large/ small obstacles, is recommended to determine limitation parameters in which plausibility of the model is preserved.

5. References

Van Hal, B.; Hepberger, A.; Priebsch, H.-H.; Desmet, W.; Sas, P.: High performance implementation and conceptual development of the wave based method for steady-state dynamic analysis of acoustic problems. In: Proceedings of ISMA 2002 – Volume II

Pluymers, B.: Wave based modelling methods for steady-state vibro-acoustics [PhD]. Katholieke Universiteit Leuven, Leuven 2006

Bergen, B.: Wave Based Modelling Techniques for Unbounded Acoustic Problems [PhD]. Katholieke Universiteit Leuven, Leuven 2011

Van Genechte, B.; Vergote, K.; Vandepitte, D.; Desmet, W.: A multi-level wave based numerical modelling framework for the steady-state dynamic analysis of bounded Helmholtz problems with multiple inclusions. In: Computer methods in applied mechanics and engineering 2010-06, Vol. 199(29)

Englert, H.: Wave Based Method: Introductory Example. (Matlab script). For: Technical University of Munich Department of Civil Engineering, Chair of Structural Mechanics

Seitenhierarchie

Implementation of a multi-level-approach with the Wave Based Method