Multiharmonic analysis for nonlinear acoustics with different scales

TL;DR

This paper develops multiharmonic asymptotic models for nonlinear acoustic wave propagation, capturing boundary layer effects and harmonic interactions, enabling efficient frequency domain simulations compared to traditional time-domain methods.

Contribution

It introduces novel multiharmonic asymptotic expansions for nonlinear acoustics that incorporate viscous boundary layers and nonlinear harmonic interactions, simplifying numerical computations.

Findings

Models accurately approximate wave behavior at different harmonic orders.

Frequency domain models are computationally more efficient than time-domain simulations.

Numerical experiments confirm the models' validity and accuracy.

Abstract

The acoustic wave-propagation without mean flow and heat flux can be described in terms of velocity and pressure by the compressible nonlinear Navier-Stokes equations, where boundary layers appear at walls due to the viscosity and a frequency interaction appears, $\textit{i.e.}$ sound at higher harmonics of the excited frequency $\omega$ is generated due to nonlinear advection. We use the multiharmonic analysis to derive asymptotic expansions for small sound amplitudes and small viscosities both of order $\varepsilon^2$ in which velocity and pressure fields are separated into far field and correcting near field close to walls and into contributions to the multiples of $\omega$. Based on the asymptotic expansion we present approximate models for either the pressure or the velocity for order $0$, $1$ and $2$, in which impedance boundary conditions include the effect of viscous boundary layers and contributions at frequencies $0$ and $2\cdot\omega$ depend nonlinearly on the approximation at frequency $\omega$. In difference to the Navier-Stokes equations in time domain, which has to be resolved numerically with meshes adaptively refined towards the wall boundaries and explicit schemes require the use of very small time steps, the approximative models can be solved in frequency domain on macroscopic meshes. We studied the accuracy of the approximated models of different orders in numerical experiments comparing with reference solutions in time-domain.

Tables

Table 1. Table 1: Pressure approximations p k ε , N subscript superscript 𝑝 𝜀 𝑁 𝑘 p^{\varepsilon,N}_{k} and velocity approximations 𝐯 k ε , N subscript superscript 𝐯 𝜀 𝑁 𝑘 {\mathbf{v}}^{\varepsilon,N}_{k} can be computed separately from each other (except for static mode for k = 0 𝑘 0 k=0 ) and sequentially in the mode index  k 𝑘 k . Then, velocity approximations 𝐰 ε , N superscript 𝐰 𝜀 𝑁 {\mathbf{w}}^{\varepsilon,N} are directly deduced from p k ε , N subscript superscript 𝑝 𝜀 𝑁 𝑘 p^{\varepsilon,N}_{k} and pressure approximations q ε , N superscript 𝑞 𝜀 𝑁 q^{\varepsilon,N} from 𝐯 ε , N superscript 𝐯 𝜀 𝑁 {\mathbf{v}}^{\varepsilon,N} . The velocity and pressure approximations for the excitation frequency ( k = 1 𝑘 1 k=1 ) coincides with the respective solutions for the linear case up to order 1 (see  [ 20 ] ).

order	pressure						velocity
$N=0$			$p_1^{\epsilon ,0}$				${𝐯}_0^{\epsilon ,0}$	$⟸$	${𝐯}_1^{\epsilon ,0}$
$N=1$			$p_1^{\epsilon ,1}$				${𝐯}_0^{\epsilon ,1}$	$⟸$	${𝐯}_1^{\epsilon ,1}$
$N=2$	$p_0^{\epsilon ,2}$	$⟸$	$p_1^{\epsilon ,2}$	$⟹$	$p_0^{\epsilon ,2}$		${𝐯}_0^{\epsilon ,2}$	$⟸$	${𝐯}_1^{\epsilon ,2}$	$⟹$	${𝐯}_2^{\epsilon ,2}$

Equations

$${\partial }_t\mathbf{v}+(\mathbf{v}\cdot \nabla )\mathbf{v}+\nabla p-\nu \Delta \mathbf{v}$$

$${\partial }_{t}p+c^2\mathrm{div}\mathbf{v}+\mathrm{div}(p\mathbf{v})$$

$$\mathbf{v}$$

$$\mathbf{f}(t,\mathbf{x})=\frac{1}{2}\left(\mathbf{f}(\mathbf{x})\exp (-\mathrm{i}\omega t)+\overline{\mathbf{f}(\mathbf{x})}\exp (\mathrm{i}\omega t)\right),$$

$$\mathbf{v}(t+\frac{2\pi }{\omega },\mathbf{x})=\mathbf{v}(t,\mathbf{x}),p(t+\frac{2\pi }{\omega },\mathbf{x})=p(t,\mathbf{x}).$$

$$\mathbf{v}(t,\mathbf{x})$$

$$p(t,\mathbf{x})$$

$${\mathcal{L}}_k(\mathbf{V},P)(\mathbf{x})+{\mathcal{N}}_k(\mathbf{V},P)(\mathbf{x})$$

$$\mathbf{V}$$

$${\mathcal{L}}_k(\mathbf{V},P)$$

$${\mathcal{N}}_0(\mathbf{V},P)$$

$${\mathcal{N}}_k(\mathbf{V},P)$$

$${\mathbf{v}}_0(\mathbf{x})=\frac{\omega }{2\pi }{\int }_0^{2\pi /\omega }\mathbf{v}(t,\mathbf{x})dt,{\mathbf{v}}_k(\mathbf{x})=\frac{\omega }{\pi }{\int }_0^{2\pi /\omega }\mathbf{v}(t,\mathbf{x})\exp (\mathrm{i}k\omega t)dt,k⩾1,$$

$$p_0(\mathbf{x})=\frac{\omega }{2\pi }{\int }_0^{2\pi /\omega }p(t,\mathbf{x})dt,p_k(\mathbf{x})=\frac{\omega }{\pi }{\int }_0^{2\pi /\omega }p(t,\mathbf{x})\exp (\mathrm{i}k\omega t)dt,k⩾1.$$

$${\mathcal{L}}_k^{\epsilon }({\mathbf{V}}^{\epsilon },P^{\epsilon })(\mathbf{x})+{\mathcal{N}}_k({\mathbf{V}}^{\epsilon },P^{\epsilon })(\mathbf{x})={\epsilon }^2\left(\begin{array}{c}{\mathbf{f}}_0(\mathbf{x})+\epsilon {\mathbf{f}}_1(\mathbf{x})+{\epsilon }^2{\mathbf{f}}_2(\mathbf{x})\\ 0\end{array}\right){\delta }_{k=1},$$

$${\mathcal{L}}_k^{\epsilon }({\mathbf{V}}^{\epsilon },P^{\epsilon })=\left(\begin{array}{rl}& -\mathrm{i}k\omega {\mathbf{v}}_k^{\epsilon }-{\epsilon }^2{\nu }_0\Delta {\mathbf{v}}_k^{\epsilon }+\nabla p_k^{\epsilon }\\ & -\mathrm{i}k\omega p_k^{\epsilon }+c^2\mathrm{div}{\mathbf{v}}_k^{\epsilon }\end{array}\right).$$

$$\mathbf{x}(\tau ,s)={\mathbf{x}}_{\partial \Omega }(\tau )-s\mathbf{n}({\mathbf{x}}_{\partial \Omega }(\tau ))$$

$$({\partial }_{\tau }u)(\mathbf{x}):={\partial }_{\tau }u({\mathbf{x}}_{\partial \Omega }(\tau )).$$

$$\kappa ({\mathbf{x}}_{\partial \Omega }(\tau )):=\frac{x_{\partial \Omega ,1}^{\prime }(\tau )x_{\partial \Omega ,2}^{\prime \prime }(\tau )-x_{\partial \Omega ,2}^{\prime }(\tau )x_{\partial \Omega ,2}^{\prime \prime }(\tau )}{(x_{\partial \Omega ,1}^{\prime }(\tau {)}^2+x_{\partial \Omega ,2}^{\prime }(\tau {)}^2{)}^{3/2}}.$$

$${\mathbf{v}}_k^{\epsilon }\sim \sum _{j=0}^{\infty }{\epsilon }^{j+2}\left({\mathbf{v}}_k^j(\mathbf{x})+{\mathbf{v}}_{\mathrm{BL},\epsilon ,k}^j(\mathbf{x})\right),p_k^{\epsilon }\sim \sum _{j=0}^{\infty }{\epsilon }^{j+2}\left(p_k^j(\mathbf{x})+p_{\mathrm{BL},\epsilon ,k}^j(\mathbf{x})\right),$$

$${\mathbf{v}}_{\mathrm{BL},\epsilon ,k}^j(\mathbf{x})$$

$$p_{\mathrm{BL},\epsilon ,k}^j(\mathbf{x})$$

$${\Phi }_{k,\tau }^j(\tau ,0)=-v_{k,\tau }^j(\tau )$$

$${\Phi }_{k,s}^j(\tau ,0)=-v_{k,s}^j(\tau )$$

$$\left(\begin{array}{c}{\mathbf{v}}_k^{\epsilon ,N}\\ q_k^{\epsilon ,N}\end{array}\right),\left(\begin{array}{c}{\mathbf{w}}_k^{\epsilon ,N}\\ p_k^{\epsilon ,N}\end{array}\right)=\sum _{j=0}^N{\epsilon }^{j+2}\left(\begin{array}{c}{\mathbf{v}}_k^j\\ p_k^j\end{array}\right)+O({\epsilon }^{N+3}).$$

$$\Delta p_1^{\epsilon ,0}+\frac{{\omega }^2}{c^2}{p}_1^{\epsilon ,0}$$

$$\nabla p_1^{\epsilon ,0}\cdot \mathbf{n}$$

$$\Delta p_1^{\epsilon ,1}+\frac{{\omega }^2}{c^2}{p}_1^{\epsilon ,1}$$

$$\nabla p_1^{\epsilon ,1}\cdot \mathbf{n}+(1+\mathrm{i})\sqrt{\frac{\nu }{2\omega }}{\partial }_{\tau }^2{p}_1^{\epsilon ,1}$$

$$\displaystyle\Big{(}1-\frac{{\rm i}\omega\nu}{c^{2}}\Big{)}\Delta{p}^{\varepsilon,2}_{1}+\frac{\omega^{2}}{c^{2}}{p}^{\varepsilon,2}_{1}$$

$$\nabla p_1^{\epsilon ,2}\cdot \mathbf{n}+(1+\mathrm{i})\sqrt{\frac{\nu }{2\omega }}{\partial }_{\tau }^2{p}_1^{\epsilon ,2}+\frac{\mathrm{i}\nu }{2\omega }{\partial }_{\tau }(\kappa {\partial }_{\tau }{p}_1^{\epsilon ,2})$$

$$\displaystyle{p}^{\varepsilon,2}_{0}=-\frac{1}{4\omega^{2}}\big{|}{\mathbf{f}}-\nabla{p}^{\varepsilon,2}_{1}\big{|}^{2}$$

$$\Delta p_2^{\epsilon ,2}+\frac{4{\omega }^2}{c^2}{p}_2^{\epsilon ,2}$$

$$\nabla p_2^{\epsilon ,2}\cdot \mathbf{n}$$

$${\mathbf{w}}_1^{\epsilon ,0}$$

$${\mathbf{w}}_2^{\epsilon ,2}$$

Full text

Multiharmonic analysis for nonlinear acoustics with different scales

Anastasia Thöns-Zueva${}^{a,\,b}$, Kersten Schmidt${}^{a,\,b,\,c}$, Adrien Semin${}^{a,\,b}$

$a$: Research center Matheon, 10623 Berlin, Germany

$b$: Institut für Mathematik, Technische Universität Berlin, Straße des 17. Juni 136, 10623 Berlin, Germany

$c$: Brandenburgische Technische Universität Cottbus-Senftenberg, Institut für Mathematik, Platz der deutschen Einheit 1, 03046 Cottbus, Germany

Corresponding author: Anastasia Thöns-Zueva, Institut für Mathematik, Technische Universität Berlin, Berlin, Germany

Address: Technische Universität Berlin, Sekretariat MA 6-4, Straße des 17. Juni 136, D-10623 Berlin

E-mail: zueva@math.tu-berlin.de

Tel: +49 (0)30 314 - 25192

**Abstract

**The acoustic wave-propagation without mean flow and heat flux can be described in terms of velocity and pressure by the compressible nonlinear Navier-Stokes equations, where boundary layers appear at walls due to the viscosity and a frequency interaction appears, i. e.  sound at higher harmonics of the excited frequency $\omega$ is generated due to nonlinear advection. We use the multiharmonic analysis to derive asymptotic expansions for small sound amplitudes and small viscosities both of order $\varepsilon^{2}$ in which velocity and pressure fields are separated into far field and correcting near field close to walls and into contributions to the multiples of $\omega$. Based on the asymptotic expansion we present approximate models for either the pressure or the velocity for order [math], $1$ and $2$, in which impedance boundary conditions include the effect of viscous boundary layers and contributions at frequencies [math] and $2\cdot\omega$ depend nonlinearly on the approximation at frequency $\omega$. In difference to the Navier-Stokes equations in time domain, which has to be resolved numerically with meshes adaptively refined towards the wall boundaries and explicit schemes require the use of very small time steps, the approximative models can be solved in frequency domain on macroscopic meshes. We studied the accuracy of the approximated models of different orders in numerical experiments comparing with reference solutions in time-domain.

**Keywords

**Acoustic wave propagation, Singularly perturbed PDE, Impedance Boundary Conditions, Asymptotic Expansions.

**AMS subject classification

**35C20, 41A60, 42A16, 35Q30, 76D05

1 Introduction

In this article we continue investigating the acoustic equations in the framework of Landau and Lifschitz [13] as a perturbation of the Navier-Stokes equations around a stagnant uniform fluid where heat flux is not taken into account. The aim of this study is to take into account nonlinear advection behaviour as well as viscous effects in the boundary layer near rigid walls. The governing equations in time domain similar to the works of Tam et al. [24, 25], but for the case of isothermal process, i. e.  pressure over density is constant over space, may be written as

[TABLE]

where ${\mathbf{v}}$ is the acoustic velocity, $p=p^{\prime}/\rho_{0}$ with $p^{\prime}$ being the acoustic pressure and $\rho_{0}>0$ being mean density, $c$ is the speed of sound, and $\nu>0$ is the kinematic viscosity. The introduction of $p$ instead of $p^{\prime}$ is only to simplify the equations by removing the constant $\rho_{0}$ and we will regard $p$ as pressure and approximations to $p$ as pressure approximations. In the momentum equation (1.1a) with some known source term ${\mathbf{f}}$ the viscous dissipation in the momentum as well as the advection nonlinear term are not neglected as we consider near wall regions where the derivatives of the acoustic velocity are rapidly increasing and might be crucial. The continuity equation (1.1b) relates the acoustic pressure to the divergence of the acoustic velocity. The system is completed by no-slip boundary conditions.

For gases the viscosity $\nu$ is very small and leads to viscosity boundary layers close to walls. These boundary layers are difficult to resolve in direct numerical simulations. Nevertheless, they have an essential influence on the absorption properties. Mainly based on experiments the physical community has introduced slip boundary conditions for the tangential component of the velocity, also known as wall laws, see for example [11, 16, 17]. For gases with small viscosity the Helmholtz equation can be completed by viscosity dependent boundary conditions [1] to obtain an approximation of high accuracy for which the boundary layers do not have to be resolved by finite element meshes [12].

In the earlier works we studied the linear acoustic equations taking into account viscous effects in the boundary layer near rigid walls. In [21] we derived a complete asymptotic expansion for the problem based on the technique of multiscale expansion in powers of $\sqrt{\eta}$, where $\eta$ is the dynamic viscosity and $\nu=\eta/\rho$. This asymptotic expansion was rigorously justified with optimal error estimates. In [20] we proposed and justified (effective) impedance boundary conditions for the velocity as well as the pressure for possibly curved boundaries.

In case of stable periodic oscillations in nonlinear dynamical systems the harmonic balance principal is used [15, 23, 28]. It is described as a linear combination of a wave of the excitation frequency and its harmonics. The referred method was presented as multiharmonic analysis for the modelling of nonlinear magnetic materials  [3, 4]. Its stability has been demonstrated within the eddy current model. For nonlinear Hamiltonian systems with a simple oscillator a special case of multiharmonic analysis, the so-called modulated Fourier expansion [10], has been well developed. The multiharmonic analysis as a method in the frequency domain is especially attractive for nonlinear acoustics. This application has not previously been investigated, either numerically or with asymptotic expansions. In the current work we restrict ourselves to the case of small sound amplitudes that are of order $O({\nu})$.

The article is ordered as follows. In Sec. 2 we introduce a frequency domain system for the quasi-stationary solution of the nonlinear acoustic wave propagation problem using the multiharmonic analysis. Moreover, the main ideas of the multiscale expansions separating far field and boundary layer contributions are introduced, that lead to the effective systems with impedance boundary conditions, both for the velocity and the pressure, that are finally introduced as main results of the paper. The far field and boundary layer terms of the multiscale expansion and the effective systems are derived in Sec. 3. Finally, in Sec. 4 we verify the effective systems by numerical computations using high-oder finite elements.

2 Multiharmonic analysis, multiscale expansion and approximative models

2.1 Multiharmonic analysis for the nonlinear system

In many acoustic applications the source is of one single frequency $\omega>0$ and so of the form

[TABLE]

with ${\mathbf{f}}:\Omega\to{\mathbb{C}}^{2}$ being complex valued. Then, we assume that the solution $({\mathbf{v}},p)$ of (1.1) tends to a quasi-stationary solution that is periodic in time with a period $T=\tfrac{2\pi}{\omega}$  which we denote by $({\mathbf{v}},p)$ again, i. e.

[TABLE]

This does not mean in general that one obtains a mono-frequency solution of the same form as the source, but as the problem involves only linear and quadratic terms its solution can be written as combination of all the harmonics $\cos(k\omega t),k=0,1,\ldots$ and $\sin(k\omega t),k=1,2,\dots$.

[TABLE]

which is called multiharmonic ansatz [3, 4, 2]. Inserting expansion (2.2) into the time-dependent problem (1.1) and identifying the terms corresponding to $\exp(-{\rm i}k\omega t)$, $k=0,1,\ldots$ leads to the infinite system of non-linear equations in space and frequency domain

[TABLE]

with the vectors {\mathbf{V}}=\big{(}{{\mathbf{v}}}_{0},{{\mathbf{v}}}_{1},\dots\big{)}^{\top} and P=\big{(}{p}_{0},{p}_{1},\dots\big{)}^{\top}, collecting the coefficients of the Fourier ansatz (2.2), the linear differential operators

[TABLE]

and the nonlinear differential operators

[TABLE]

Theorem 2.1**.**

The existence and uniqueness of a quasi-stationary solution $({\mathbf{v}},p)$ of (1.1) satisfying (2.1) is equivalent to the existence and uniqueness of a solution $({\mathbf{v}}_{0},{\mathbf{v}}_{1},\dots,p_{0},p_{1},\dots)$ of (2.3).

Proof.

First, let $({\mathbf{v}},p)$ be the unique solution of (1.1). Defining the vector $({\mathbf{v}}_{0},{\mathbf{v}}_{1},\dots,p_{0},p_{1},\dots)$ by

[TABLE]

the decomposition (2.2) holds and by construction $({\mathbf{v}}_{0},{\mathbf{v}}_{1},\dots,p_{0},p_{1},\dots)$ is solution of (2.3). Moreover, [2, Lemma 3.3] implies that the vector is unique.

Now, let $({\mathbf{v}}_{0},{\mathbf{v}}_{1},\dots,p_{0},p_{1},\dots)$ be the unique solution of (2.3). It is easy to see that $({\mathbf{v}},p)$ defined by (2.2) satisfies (1.1). If (1.1) would admit another solution $({\mathbf{w}},q)$, then its coefficients ${\mathbf{w}}_{k}$, $q_{k}$ defined similarly as in (2.4) would fulfil (1.1) as well, which is a contradiction to the assumption of unicity. ∎

2.2 Asymptotic ansatz for small sound amplitude and viscosity

To investigate acoustic velocity ${\mathbf{v}}$ and acoustic pressure $p$ for small acoustic excitation and viscosity we introduce a small parameter $\varepsilon\in{\mathbb{R}}^{+}$ and replace the acoustic source ${\mathbf{f}}$ by $\varepsilon^{2}\,\sum_{j=0}^{\infty}\varepsilon^{j}{\mathbf{f}}_{j}$, where each term ${\mathbf{f}}_{j}$ is independent of $\varepsilon$, and the viscosity $\nu$ by $\varepsilon^{2}\nu_{0}$ with $\nu_{0}\in{\mathbb{R}}^{+}$. Moreover, we consider the leading order source term $\varepsilon^{2}\,{\mathbf{f}}_{0}$ to be $\operatorname{curl}_{2D}$-free, i. e.  $\operatorname{curl}_{2D}{\mathbf{f}}_{0}=0$, having in mind that ${\mathbf{f}}_{0}=\nabla p_{0}$ where the pressure $p_{0}$ corresponds to a solution of a linear and inviscid wave equation. In addition, we assume for simplicity the source to disappear on the boundary. The impedance boundary conditions with additional terms due to more general source functions will be given in the Appendix A.3.

For these small acoustic excitations the leading part of the solution satisfies a linear equation in the whole domain as considered in [21] and the nonlinearity will come into play on a higher order. The small viscosities on the other hand leads to boundary layers whose thickness becomes proportional to $\varepsilon$. Indicating their dependency on $\varepsilon$ we will label the acoustic velocity ${\mathbf{v}}^{\varepsilon}$ and the acoustic pressure $p^{\varepsilon}$ with a superscript $\varepsilon$. They are described by the system

[TABLE]

with the vectors {\mathbf{V}}^{\varepsilon}=\big{(}{{\mathbf{v}}}^{\varepsilon}_{0},{{\mathbf{v}}}^{\varepsilon}_{1},\dots\big{)}^{\top}, P^{\varepsilon}=\big{(}p^{\varepsilon}_{0},p^{\varepsilon}_{1},\dots\big{)}^{\top} of velocity and pressure coefficients and the linear differential operators

[TABLE]

In the following we specify first the domain and its boundary before we introduce the ansatz for an asymptotic expansion with far field terms and near field correctors and their coupling conditions.

The geometrical setting

Let $\Omega\subset{\mathbb{R}}^{2}$ be a bounded domain with smooth boundary $\partial\Omega$. The boundary shall be described by a mapping ${\mathbf{x}}_{\partial\Omega}:\tau\in\Gamma\to{\mathbb{R}}^{2}$ from a one-dimensional reference domain $\Gamma\subset{\mathbb{R}}$. We assume the boundary to be $C^{\infty}$ such that points in some neighbourhood $\Omega_{\Gamma}$ of $\partial\Omega$ can be uniquely written as

[TABLE]

where ${\mathbf{n}}$ is the outer normalised normal vector and $s$ the distance from the boundary.

Without loss of generality we assume $|{\mathbf{x}}_{\partial\Omega}^{\prime}(\tau)|=1$ for all $\tau\in\Gamma$. The orthogonal unit vectors in these tangential and normal coordinate directions are ${\mathbf{e}}_{\tau}=-{\mathbf{n}}^{\bot}$, where we use the notation ${\mathbf{u}}^{\bot}=(u_{2},-u_{1})^{\top}$ for a turned vector clockwise by $90^{\circ}$, and ${\mathbf{e}}_{s}=-{\mathbf{n}}$. This allows us to write the tangential derivative $\nabla u({\mathbf{x}})\cdot{\mathbf{e}}_{\tau}$ of a function $u\in C^{1}(\Omega)$ with abuse of notation as

[TABLE]

Moreover, the curvature $\kappa$ on the boundary $\partial\Omega$ is given by

[TABLE]

Asymptotic ansatz.

Within this article we consider the acoustic source of the same order as the boundary layer, i. e. ${\mathbf{f}}={\mathbf{f}}^{\varepsilon}=O(\varepsilon^{2})$. In the linear model the resulting acoustic velocity and pressure are of the same order and for the considered nonlinear model the same is true. The solution ${\mathbf{V}}^{\varepsilon}$, $P^{\varepsilon}$ of (2.5) should be approximated by a two-scale asymptotic expansion in the framework of Vishik and Lyusternik [27] and for each coefficient we take the ansatz

[TABLE]

where ${\mathbf{v}}^{j}_{k}$ and $p^{j}_{k}$ are the far field velocity and pressure of order $j$ and ${\mathbf{v}}^{j}_{{\mathrm{BL},\varepsilon},k}$ and $p^{j}_{{\mathrm{BL},\varepsilon},k}$ represent the respective near field velocity and pressure. They are seeked in scaled coordinate $S(s)=\tfrac{s}{\varepsilon}$ of the local normalised coordinate system (2.6) in the form

[TABLE]

taking into account the fact that the boundary layer thickness scales linearly with $\varepsilon$. For the desired decay properties we require the near field terms $\Phi^{j}_{k,\tau}(\tau,S)$, $\Phi^{j}_{k,s}(\tau,S)$ and $\Pi^{j}_{k}(\tau,S)$ as well as their higher derivatives to vanish with $S\to\infty$. The subscript $\cdot_{\mathrm{BL},\varepsilon}$ stands for “boundary layer” expressing the nature of the near field terms that they are essentially defined in a small layer close the boundary. Indeed the equality (2.9) can be assumed to be true only in an $O(1)$ neighbourhood of the boundary in which the local coordinate system (2.6) is defined. Outside this neighbourhood the expression on the right hand sides of (2.9), that decaying exponentially in an $O(\varepsilon)$ distance, is multiplied with a smooth cut-off function such that the product is exactly zero where the local coordinate system is not defined.

In the linear case [21] the near field velocity turned out to be divergence free such that there is no boundary layer for the pressure. Due to the coupling of the velocity and pressure by the nonlinear terms this property can not be assumed in general. However, we will see in our analysis that the near field pressure terms $\Pi^{j}_{{\mathrm{BL},\varepsilon},k}$ vanish at least up to order 2 and up to this order the resulting near field terms for the frequency $\omega$ of the excitation remain exactly the same as for the linear system.

Coupling of far and near field by the no-slip boundary conditions

By the homogeneous Dirichlet boundary condition the tangential trace \big{(}{\mathbf{v}}^{\varepsilon}_{k}+{\mathbf{v}}^{{\mathrm{BL},\varepsilon}}_{k}\big{)}\cdot{\mathbf{e}}_{\tau} and normal trace \big{(}{\mathbf{v}}^{\varepsilon}_{k}+{\mathbf{v}}^{{\mathrm{BL},\varepsilon}}_{k}\big{)}\cdot{\mathbf{n}} vanish for any $k$ and separately in the orders in $\varepsilon$, cf (1.1c), therefore the traces of the far field have to fulfil the conditions

[TABLE]

The far and near field terms will be derived order by order up order 2 in Section 3 as well as the effective systems with impedance boundary conditions that we will present already in the following subsection.

2.3 Effective systems with impedance boundary conditions

In this section we present effective models of order $N=0$, $1$ and $2$ for approximative far field solutions ${\mathbf{V}}^{\varepsilon,N}=({\mathbf{v}}^{\varepsilon,N}_{0},{\mathbf{v}}^{\varepsilon,N}_{1},\ldots)^{\top}$, $P^{\varepsilon,N}=(p^{\varepsilon,N}_{0},p^{\varepsilon,N}_{1},\ldots)^{\top}$, in which the nonlinear and viscous behaviour in the layers close to the boundary are incorporated with impedance boundary conditions. The steps for deriving the systems for approximative velocity and pressure will follow in Section 3.4. Contrary to the original system (1.1) in time domain or its multiharmonic approximation (2.5), for which all modes of velocity and pressure couple, the approximative pressure and velocity coefficients decouple for all modes $k>0$. Therefore, we introduce separately systems for pressure coefficients $p^{\varepsilon,N}_{k}$ only, where associated velocity coefficients ${\mathbf{w}}^{\varepsilon,N}_{k}$ are defined afterwords as a function of the pressure, and systems for velocity coefficients ${\mathbf{v}}^{\varepsilon,N}_{k}$ only, where associated pressure coefficients $q^{\varepsilon,N}_{k}$ follow directly. Only for the static mode $k=0$ we have coupled velocity and pressure systems. In general, the directly defined pressure coefficients $p^{\varepsilon,N}_{k}$ and the pressure coefficients $q^{\varepsilon,N}_{k}$ computed from the velocity may differ as well as the two velocity approximations ${\mathbf{v}}^{\varepsilon,N}_{k}$ and ${\mathbf{w}}^{\varepsilon,N}_{k}$. We also distinguish the two for the static mode $k=0$ even so here velocity and pressure coefficients are defined in a coupled system as the right hand side of this system depends on the different approximations.

Derived from the asymptotic expansion, both approximative far field solutions for velocity and pressure order $N$ shall be close to the respective far field expansion of order $N$, i. e. for the frequency mode $k$ we expect that

[TABLE]

Even equally important the asymptotic regime of small sound amplitudes leads to an iterative procedure to obtain the coefficients for different modes $k$ (see Table 1). In general, the coefficients for $k=1$ can be defined independently and the neighbouring modes for $k=0$ and $k=2$ follow. Moreover, up to order $2$ there are no modes for $k>2$ as indicating that the response at the higher harmonics $3\,\omega,4\,\omega,\ldots$ are more than two orders in $\varepsilon$ smaller than the excitation amplitude. In general, for approximation of order $N$ we have only the modes $k=0,1,\ldots,\lceil\frac{N+1}{2}\rceil$.

2.3.1 Systems for the pressure

Here we present the approximative models for the far field pressure. This is different to the original equations in which no boundary conditions for the pressure, but for both velocity components, are imposed for each order. Approximative velocities can be computed a-posteriori (see Sec. 2.3.2).

Order $O(\varepsilon^{2})$

The approximative model for the pressure in frequency of the excitation $1\cdot\omega$ is given by a linear system

[TABLE]

All the terms $p_{k}^{\varepsilon,0}$ for $k\neq 1$ are zero, meaning that the limit acoustic pressure is exactly as in the linear case.

Order $O(\varepsilon^{3})$

The approximative model in frequency $1\cdot\omega$ is given by the linear system

[TABLE]

Again, all the terms $p_{k}^{\varepsilon,1}$ for $k\neq 1$ are zero and the resulting acoustic pressure approximation is exactly as in the linear case. The impedance boundary conditions (2.13b) are of Wentzell type. See [5, 18] for the functional framework and variational formulation.

Order $O(\varepsilon^{4})$

For frequency $1\cdot\omega$ the pressure of order 2 is solution of

[TABLE]

Even for $N=2$ the nonlinear terms do not affect the pressure approximation in frequency of excitation, which coincides with the approximations in the linear case and are also obtained via a system decoupled from the velocity. However, in this order of approximation the first time other frequency modes come into play, namely that for the frequency $0\cdot\omega$, a so called acoustic streaming [14], and for the frequency $2\cdot\omega$. The acoustic pressure at frequency $0\cdot\omega$ is explicitly defined by the algebraic equation

[TABLE]

and the one at frequency $2\cdot\omega$ by the Helmholtz equation

[TABLE]

For a well-posed definition of $p^{\varepsilon,2}_{2}$ the source function ${\mathbf{f}}$ has to be continously differentiable and also the pressure approximation $p^{\varepsilon,2}_{1}$ needs higher regularity. Note, that the right hand side of (2.16a) can be simplified using the fact that $p^{\varepsilon,2}_{1}$ is solution to an Helmholtz problem (see Appendix A.4). For a numerical approximation with $C^{0}$-continuous finite elements, for which this regularity is only attained approximately, the right hand side can be evaluated as the projection of the pressure gradient $\nabla p^{\varepsilon,2}_{1}$ to continuous vector fields.

2.3.2 Post-processing of velocity from systems for the pressure

When the far field pressure is computed we may obtain a-posteriori approximations to the far field velocity to the respective order.

The far field velocities at frequency $1\cdot\omega$ are defined at the different approximation orders by

[TABLE]

and those at frequency $2\cdot\omega$ at order $2$ by

[TABLE]

For the frequency $0\cdot\omega$ a far field velocity approximation ${\mathbf{w}}^{\varepsilon,1}_{0}$ of order 1 can be obtained as solution of linear Stokes system similarly to (2.21) in the following subsection that is directly for a velocity approximation, however, using ${\mathbf{w}}^{\varepsilon,1}_{1}$ on its right hand side. Likewise, a far field velocity approximation ${\mathbf{w}}^{\varepsilon,2}_{0}$ for order 2 can be defined by a nonlinear Navier-Stokes like system as (2.23) that depends on ${\mathbf{w}}^{\varepsilon,2}_{1}$.

The far field velocity can be used as approximation away from the boundary and has to be corrected by a near field velocity approximation (see Sec. 2.3.5).

2.3.3 Systems for the velocity

Here we propose approximative models directly for the far field velocity. For each order an approximative pressure can be computed afterwards (see Sec. 2.3.4) as well as a near field velocity approximation (see Sec. 2.3.5).

Order $O(\varepsilon^{2})$

The limit model is given by a linear system in frequency of excitation

[TABLE]

and all other terms ${\mathbf{v}}_{k}^{\varepsilon,0}$, $k\neq 1$ are zero. So the limit acoustic velocity coincides with the one in the linear case.

Order $O(\varepsilon^{3})$

In frequency of excitation the approximative model is given by

[TABLE]

and there is a non-zero acoustic streaming velocity at frequency $0\cdot\omega$ that satisfies the Stokes system

[TABLE]

The purely real right hand side of (2.21) implies that its solution $({{\mathbf{v}}}^{\varepsilon,1}_{0},q^{\varepsilon,3}_{0})$ is purely real. Note that $q^{\varepsilon,3}_{0}$ is not only a Lagrange multiplier but a higher order approximation of the pressure at zero frequency.

Order $O(\varepsilon^{4})$

The approximative model in frequency $1\cdot\omega$ is defined by

[TABLE]

[TABLE]

Again, the solution $({{\mathbf{v}}}^{\varepsilon,2}_{0},q^{\varepsilon,4}_{0})$ of (2.23) is purely real and $q^{\varepsilon,4}_{0}$ is a pressure approximation of higher order, where q_{0}^{\varepsilon,4}=-\frac{1}{4}\big{|}{{\mathbf{v}}}^{\varepsilon,2}_{1}\big{|}^{2}+O(\varepsilon^{5}). At frequency $2\cdot\omega$ a velocity approximation satisfies the Helmholtz equation

[TABLE]

Remark 2.2**.**

Note, that $({{\mathbf{v}}^{\varepsilon,2}_{0}}\cdot\nabla){\mathbf{v}}^{\varepsilon,2}_{0}=({{\mathbf{v}}^{\varepsilon,1}_{0}}\cdot\nabla){\mathbf{v}}^{\varepsilon,1}_{0}+O(\varepsilon^{5})$. Therefor instead of solving the nonlinear system (2.23a)–(2.23c) one could, first, find ${{\mathbf{v}}^{\varepsilon,1}_{0}}$ by solving the linear system (2.21a)–(2.21c), and then substitute it in (2.23a)–(2.23c) for ${{\mathbf{v}}^{\varepsilon,2}_{0}}$ in the advection term. That will lead again to a linear system.

2.3.4 Post-processing of pressure from systems for the velocity

When the far field velocity approximation is computed we may obtain a-posteriori an associated far field pressure approximation for the frequencies $1\cdot\omega$ and $2\cdot\omega$. The approximations for frequency $1\cdot\omega$ are given by

[TABLE]

and the approximation of order 2 for frequency $2\cdot\omega$ by

[TABLE]

Moreover, a pressure approximation of order 2 at frequency $0\cdot\omega$ is given by

[TABLE]

2.3.5 Post-processing of a near field velocity

Close to the wall the far field velocity approximations ${\mathbf{V}}^{\varepsilon,N}$ have to be corrected by boundary layer functions in tangential as well as normal direction

[TABLE]

where $\chi$ is an admissible cut-off function (see [21]) that takes the constant value $1$ in some subset of $\Omega_{\Gamma}$, $s$ is the distance function to the boundary, i. e. there exists for each point ${\mathbf{x}}\in\Omega_{\Gamma}$ a base point ${\mathbf{x}}_{\partial\Omega}\in\partial\Omega$ such that ${\mathbf{x}}={\mathbf{x}}_{\partial\Omega}+s({\mathbf{x}}){\mathbf{n}}^{\bot}({\mathbf{x}}_{\partial\Omega})$, and the operators ${\mathbf{H}}^{\ell}:(C^{\infty}(\partial\Omega))^{\infty}\to(C^{\infty}(\Omega_{\Gamma}))^{\infty}$ with ${\mathbf{H}}^{\ell}=(0,{\mathbf{h}}^{\ell}_{1},{\mathbf{h}}^{\ell}_{2},\ldots)^{\top}$ that are acting on the tangential velocity traces $V^{\varepsilon,N}_{\tau}=(v^{\varepsilon,N}_{0,\tau},v^{\varepsilon,N}_{1,\tau},\ldots)^{\top}$ with $v^{\varepsilon,N}_{k,\tau}({\mathbf{x}}_{\partial\Omega}):={\mathbf{v}}^{\varepsilon,N}_{k}({\mathbf{x}}_{\partial\Omega})\cdot{\mathbf{n}}^{\bot}({\mathbf{x}}_{\partial\Omega})$ where we note that $v^{\varepsilon,N}_{k,\tau}=0$ for $k>\lceil\frac{N+1}{2}\rceil$. Hence, to define ${\mathbf{V}}^{{\mathrm{BL},\varepsilon},N}$ for $N=0,1,2$ we state the operators

[TABLE]

where we note that $\kappa$, ${\mathbf{n}}^{\bot}$, ${\mathbf{n}}$ and $v^{\varepsilon,N}_{1,\tau}$ are functions of the base point ${\mathbf{x}}_{\partial\Omega}$ of ${\mathbf{x}}$ and $\partial_{\tau}$ is the tangential derivative defined in (2.7).

3 Derivation of terms of multiscale expansion and effective systems

In Sec. 2.2 we have introduced the ansatz of the two-scale expansion (2.8), which expresses an approximation to the exact solution as a two-scale decomposition into far field terms, modelling the macroscopic picture of the solution, which are corrected in the neighbourhood of the boundary by near field terms. To separate the two scales we use the technique of multiscale expansion as described in Sec. 2.2, which defines the near field terms in a local normalised coordinate system (2.6) such that they decay rapidly away from the wall and are set to zero where the local coordinate system is not defined (using a cut-off function). In the following we define the terms of asymptotic expansion (2.8) order by order.

3.1 Correcting near field

In this section we will give the near field equations and their solutions up to order 2. They are derived such that the near field velocity and pressure expansions (2.9) inserted into (2.5) leave a residual as small as possible in powers of $\varepsilon$ and that the sum of tangential far and near field velocity vanishes at the boundary. The general form of the near field equations of any order and in any frequency can be found in the Appendix A.2.

The near field terms of order $O(\varepsilon^{2})$.

The near field equation for $j=0$ in frequency $1\cdot\omega$ yields

[TABLE]

It is easy to see that its unique solution together with the coupling condition for far and near fields (2.10) and decay condition for the near field is given by

[TABLE]

This is the dominating boundary layer term close to the wall.

The near field terms of order $O(\varepsilon^{3})$.

The near field equations for $j=1$ in frequency $1\cdot\omega$ are given by

[TABLE]

which unique solution, using the terms in (3.1) together with the coupling condition, is

[TABLE]

The near field terms of order $O(\varepsilon^{4})$.

The near field equations for $j=2$ in frequency $1\cdot\omega$ are given by

[TABLE]

which unique solution, using the terms in (3.1) and (3.1) together with the coupling condition, is

[TABLE]

In frequency $2\cdot\omega$ the first non trivial terms appear for $j=2$ with the near field equations given by

[TABLE]

Its unique solution together with the coupling condition is given by

[TABLE]

For frequency $0\cdot\omega$ the unique solution is the trivial solution at least up to order $j=2$, i. e. the boundary layer disappears.

3.2 Far field velocity terms

In the following section we will derive the terms of asymptotic expansion for the far field velocity up to order 2. The resulting expressions in frequency $1\cdot\omega$ are exactly the expressions for the linear case which are derived and analysed in [21]. The expressions for frequencies $0\cdot\omega$ and $2\cdot\omega$ are only due to the nonlinear advection term and do not appear for the linear case. The general form of the far field equations of any order and in any frequency can be found in the Appendix A.1.

Approximation of order $O(\varepsilon^{2})$.

The limit model for the far field velocity in frequency $1\cdot\omega$ is given by

[TABLE]

The far field approximation for frequency $0\cdot\omega$ is given by the stationary incompressible Navier-Stokes equations

[TABLE]

Approximation of order $O(\varepsilon^{3})$.

The first correcting terms, i. e. for $j=1$, for frequency $1\cdot\omega$ are given by

[TABLE]

and for frequency $0\cdot\omega$ the far field approximation solves the Stokes system

[TABLE]

where we have used ${\mathbf{v}}^{0}_{0}=0$. The stationary velocity term ${\mathbf{v}}^{1}_{0}$ is coupled with the stationary pressure term $p^{3}_{0}$, however, the system is linear.

Approximation of order $O(\varepsilon^{4})$.

The next correcting terms, i. e. for $j=2$, for frequency $1\cdot\omega$ are given by

[TABLE]

where we used (3.7) and the fact that $\operatorname{curl}_{2D}{{\mathbf{v}}}^{0}_{1}=0$. By the assumption on the source function $\operatorname{curl}_{2D}{\mathbf{f}}_{0}=0$ the term $\tfrac{\nu_{0}}{\omega^{2}}\operatorname{{\bf curl}}_{2D}\operatorname{curl}_{2D}{\mathbf{v}}^{0}_{1}$ in (3.11a) disappears. For the system for the frequency $0\cdot\omega$ the far field pressure term $p^{0}_{1}$ is needed, which is obtained a-posteriori from the far field velocity ${\mathbf{v}}^{0}_{1}$ as

[TABLE]

For frequency $0\cdot\omega$ far field second correcting terms are

[TABLE]

where we used $\operatorname{curl}_{2D}{\mathbf{v}}^{2}_{1}=0$ which is due to the fact that $\operatorname{curl}_{2D}{\mathbf{v}}^{0}_{1}=0$.

Moreover, the far field velocity at frequency $2\cdot\omega$ is obtained from

[TABLE]

3.3 Far field pressure up to $2^{\mathrm{nd}}$ order

As it was mentioned before, the far field approximations in frequency $\omega$ for $j=0,1,2$ are exactly the results for the linear problem in [21]. Accordingly, we can rewrite the equations in terms of the far field pressure with the suitable boundary conditions.

Approximation of order $O(\varepsilon^{2})$.

The limit model is given by

[TABLE]

Approximation of order $O(\varepsilon^{3})$.

The first correcting terms are given by

[TABLE]

Approximation of order $O(\varepsilon^{4})$.

The next correcting terms, i. e. for $j=2$, for frequency $\omega$ are given by

[TABLE]

When the far field pressure terms for the frequency $\omega$ are computed we may obtain a posteriori the far field velocity terms by

[TABLE]

The far field approximation for frequency $0\cdot\omega$ is given by

[TABLE]

and in frequency $2\cdot\omega$ by

[TABLE]

3.4 Deriving the effective systems with impedance boundary conditions

In the previous sections we have derived the terms of the asymptotic expansions (2.8) up to order 2, which we can assemble to obtain pressure and velocity approximations of these orders. To obtain pressure approximations of order 1 two Helmholtz systems have to be solved, for order 2 these are four Helmholtz systems. To obtain velocity approximations of order 1 we need to solve three PDEs, and for order 2 these are already six. In general, the number of terms in the asymptotic expansion increase like $\frac{1}{4}N^{2}$ with the order $N$, and, hence, the number of systems to solve. In this section, we derive the effective systems given in Sec. 2.3 that are written directly for approximative solutions of order 0, 1 and 2. The approximative solutions show the same accuracy as the asymptotic expansions (2.8) but for a less computational effort as all terms for $j\leqslant N$ would have been computed at once. Here, the number of systems to solve increases only linearly with $N$ and to obtain pressure and velocity approximations of order 2 only two or three systems, respectively, have to be solved. The main idea is to combine the equations satisfied by each far field term of (2.8) and to neglect the next order terms. In this way we obtain equations satisfied by the pressure coefficients $p_{k}^{\varepsilon,N}$, where associated velocity coefficients ${\mathbf{w}}_{k}^{\varepsilon,N}$ are defined afterwords as a function of the pressure (see Sec. 3.4.1 for the first order and in Sec. 3.4.2 for the second order model), or equations satisfied by the velocity coefficients ${\mathbf{v}}_{k}^{\varepsilon,N}$, where associated pressure coefficients $q_{k}^{\varepsilon,N}$ are defined afterwords as a function of the velocity (see Sec. 3.4.3 and Sec. 3.4.4 for the first and second order model, respectively).

3.4.1 Derivation of $1^{\textrm{st}}$ order effective system

for the pressure

The derivation for the $1^{\textrm{st}}$ order effective system (2.13) for the pressure is exactly as for the linear case [20]. Adding the system (3.15) for $p^{1}_{0}$ and $\varepsilon$ times the system (3.16) for $p^{1}_{1}$ we obtain a system for the first order asymptotic expansion $\widetilde{p}_{1}^{\varepsilon,1}:=\varepsilon^{2}(p_{1}^{0}+\varepsilon p_{1}^{1})$

[TABLE]

Neglecting in (3.21b) the $O(\varepsilon^{4})$ term and replacing $\varepsilon\sqrt{\nu_{0}}$ by $\sqrt{\nu}$ gives (2.13) with Wentzel boundary conditions on the domain boundaries. Now, adding (3.18) for $j=0$ and $\eqref{eq:pj1:vj1}$ for $j=1$ multiplied by $\varepsilon$ we obtain for the first order asymptotic expansions $\widetilde{{\mathbf{v}}}^{\varepsilon,1}_{1}:={\mathbf{v}}^{0}_{1}+\varepsilon{\mathbf{v}}^{1}_{1}$

[TABLE]

where ${\mathbf{f}}^{\varepsilon,1}=\varepsilon^{2}({\mathbf{f}}_{0}+\varepsilon{\mathbf{f}}_{1})$. Noting that ${\mathbf{f}}={\mathbf{f}}^{\varepsilon,1}+O(\varepsilon^{4})$ and if $\widetilde{p}^{\varepsilon,1}_{1}=p^{\varepsilon,1}_{1}+O(\varepsilon^{4})$ holds (this is the case if $\frac{\omega^{2}}{c^{2}}$ is not a Neumann eigenvalue of $-\Delta$, see  [20]) and neglecting the $O(\varepsilon^{4})$ terms, we find that (2.17) defines a first order velocity approximation ${\mathbf{w}}_{1}^{\varepsilon,1}$. Adding the system (3.8) for ${\mathbf{v}}_{0}^{0}$ and $\varepsilon$ times the system (3.10) for ${\mathbf{v}}_{0}^{1}$ and neglecting the $O(\varepsilon^{4})$ terms we find that (2.21) defines a first order approximation stationary velocity ${\mathbf{w}}_{0}^{\varepsilon,1}$ that depends on ${\mathbf{w}}^{\varepsilon,1}_{1}$ and incorporates with a Lagrange multiplier $q_{0}^{\varepsilon,3}$ that is $O(\varepsilon^{4})$. Finally, we can reconstruct the pressure and the velocity in time by

[TABLE]

with the neglected term in the reconstructions (3.22) being in $O(\varepsilon^{4})$.

3.4.2 Derivation of $2^{\textrm{nd}}$ order effective system

for the pressure

Similarly, taking $\eqref{eq:p:0}+\varepsilon\eqref{eq:p:1}+\varepsilon^{2}\eqref{eq:p:21}$, neglecting the $O(\varepsilon^{3})$ term and using that $\nu=\varepsilon^{2}\nu_{0}$ leads to the $2^{\textrm{nd}}$ order effective system (2.14) for the pressure at frequency $\omega$. As well, we obtain a posteriori the far field velocity approximation ${\mathbf{w}}_{1}^{\varepsilon,2}$ defined by (2.17), combining $\varepsilon^{j}\eqref{eq:pj1:vj1}$ for $j=0$, $j=1$ and $j=2$ and neglecting the $O(\varepsilon^{5})$ term. Using that in the expansion $p_{0}^{0}=p_{0}^{1}=0$, taking (3.19) and neglecting the $O(\varepsilon^{3})$ term leads to

[TABLE]

and so to (2.15). Similarly, using that ${\mathbf{w}}^{\varepsilon,2}_{1}=\varepsilon^{2}{\mathbf{v}}^{0}_{1}+O(\varepsilon^{3})$, $p^{\varepsilon,2}_{1}=\varepsilon^{2}p^{0}_{1}+O(\varepsilon^{3})$ and $\frac{\omega^{2}}{c^{2}}p^{\varepsilon,2}=\varepsilon^{2}\operatorname{div}{\mathbf{v}}^{0}_{1}+O(\varepsilon^{3})$, the latter being a consequence of (3.18) for $j=0$ and (3.15), we find the $2^{\textrm{nd}}$ order effective system (2.16) for the pressure contribution at frequency $2\cdot\omega$. Then, using the equality (3.20c) and that $p^{\varepsilon,2}_{2}=\varepsilon^{{4}}p^{2}_{2}+O(\varepsilon^{{5}})$ (assuming that $\frac{2\omega}{c}$ is not a Neumann eigenvalue of $-\Delta$) and $\frac{{\rm i}}{\omega}({\mathbf{f}}-\nabla p^{\varepsilon,1})=\varepsilon^{2}{\mathbf{v}}^{0}_{1}+O(\varepsilon^{5})$ we obtain the equation (2.18) for the velocity approximation ${\mathbf{w}}^{\varepsilon,2}_{2}$ in terms of $p^{\varepsilon,2}_{1}$ and $p^{\varepsilon,2}_{2}$. Adding the system (3.8) for ${\mathbf{v}}_{0}^{0}$, $\varepsilon$ times the system (3.10) for ${\mathbf{v}}_{0}^{1}$ and $\varepsilon^{2}$ times the system (3.13) we find that (2.23) defines a first order approximation stationary velocity ${\mathbf{w}}_{0}^{\varepsilon,2}$ with a Lagrange multiplier $q_{0}^{\varepsilon,4}=p_{0}^{\varepsilon,2}+O(\varepsilon^{5})$. Finally, we can reconstruct the pressure and the velocity in time by

[TABLE]

with the neglected term in the reconstructions (3.24) being in $O(\varepsilon^{5})$.

3.4.3 Derivation of $1^{\textrm{st}}$ order effective system for

the velocity

Taking $\eqref{eq:v:0:w1}+\varepsilon\eqref{eq:v:1:w0}$ and neglecting the $O(\varepsilon^{4})$ term leads to the $1^{\textrm{st}}$ order effective system (2.20) for the velocity component to the frequency $\omega$.

Now, adding $\varepsilon^{4}\eqref{eq:v:0:w0}$ and $\varepsilon^{5}\eqref{eq:v:1:w1}$ and using that ${\mathbf{v}}^{0}_{0}=0$ we find that $\widetilde{{\mathbf{v}}}^{\varepsilon,1}_{0}:=\varepsilon^{3}{\mathbf{v}}^{1}_{0}$, $\widetilde{p}^{\varepsilon,1}_{0}:=\varepsilon^{4}(p^{2}_{0}+\varepsilon p^{3}_{0})$ solve

[TABLE]

Neglecting the $O(\varepsilon^{6})$ term on the right hand side leads to (2.21). Finally, we get similarly the equality (2.25) for the a-posteriori computed pressure approximation $q^{\varepsilon,1}_{1}$, where (3.12) and

[TABLE]

are used.

3.4.4 Derivation of $2^{\textrm{nd}}$ order effective system for

the velocity

The derivation of the systems (2.22) for ${\mathbf{v}}^{\varepsilon,2}_{1}$ and (2.23) for ${\mathbf{v}}^{\varepsilon,2}_{0}$ is similar to the respective first order systems as well as equality (2.25) for $q^{\varepsilon,2}_{1}$. The system (2.24) for the velocity contribution ${\mathbf{v}}^{\varepsilon,2}_{2}$ at frequency $2\cdot\omega$ is a direct consequence of (3.14) and the equation(2.16) for the pressure $q^{\varepsilon,2}_{2}$ is a direct consequence of (3.14c). Then, the system (2.23) for the velocity contribution ${\mathbf{v}}^{\varepsilon,2}_{0}$ at frequency $0\cdot\omega$ is derived similarly to (2.21) of order 1 using the systems (3.10) and (3.13).

4 Numerical results

We verify the derived approximative models with impedance boundary conditions on a ring domain $\Omega$ centered at $(0,0)$ whose inner radius is $R_{1}=1.5$ and outer radius is $R_{2}=2.0$. For this we choose several values for $\varepsilon$, where the viscosity $\nu=\varepsilon^{2}$ (i. e., $\nu_{0}=1$) and the source ${\mathbf{f}}$ takes the decomposition ${\mathbf{f}}(t,{\mathbf{x}})=\left(\varepsilon^{2}{\mathbf{f}}_{0}({\mathbf{x}})+\varepsilon^{3}{\mathbf{f}}_{1}({\mathbf{x}})\right)\cos(\omega t)$ with $\omega=15$, where the dominating part $\varepsilon^{2}{\mathbf{f}}_{0}$ is $\operatorname{curl}_{2D}$-free. More precisely, we take ${\mathbf{f}}_{0}={\operatorname{Re}}\nabla p_{0}$ with

[TABLE]

with the polar coordinates $(r,\phi)$ in the ring, $\lambda=4$ and $k=2.28945$ computed numerically such that the Neumann trace $\nabla p_{0}\cdot{\mathbf{n}}=0$ on the boundary $\partial\Omega$ of the ring. In this way, $p_{0}\in H^{1}(\Omega)$ is solution of the Helmholtz equation

[TABLE]

Hence, the normal component of the source ${\mathbf{f}}_{0}$ vanishes on $\partial\Omega$. As the tangential component of ${\mathbf{f}}_{0}$ does not vanish we use the formulations with additional terms ${\mathbf{f}}\cdot{\mathbf{n}}^{\bot}$ that are given in Appendix A.3. Moreover, the second term of the source is a bubble function with ${\mathbf{f}}_{1}({\mathbf{x}})=(R_{1}^{2}-r({\mathbf{x}})^{2})(R_{2}^{2}-r({\mathbf{x}})^{2})\binom{1}{1}$ that vanishes in both components on $\partial\Omega$.

We have computed numerically approximative solutions of different order by high order finite elements with curved cells using the numerical C++ library Concepts [9, 19, 7], where we use the formulation for the pressure. To estimate the modelling error of these approximative solution we compute numerically a reference solution in time-domain using a modified Crank-Nicolson scheme in which the nonlinear advection terms are discretized explicitly (see [26, 29] for similar schemes for incompressible fluids). The time-domain formulation with time step $\Delta t>0$ is in both variables, the pressure and the velocity, and given by

[TABLE]

where $({\mathbf{v}}^{\ell},p^{\ell})$ is a numerical approximation to $({\mathbf{v}}(\ell\Delta t,\cdot),p(\ell\Delta t,\cdot))$, ${\mathbf{f}}^{\ell}:={\mathbf{f}}(\ell\Delta t,\cdot)$ and ${\mathbf{v}}^{\ell+\frac{1}{2}}$, $p^{\ell+\frac{1}{2}}$, ${\mathbf{f}}^{\ell+\frac{1}{2}}$ denote the averages

[TABLE]

As initial velocity and pressure we use the solution of the linear system (without the nonlinear advection terms) and simulate for $15$ periods to $50$ depending on $\varepsilon$ to obtain an accurate approximation to the quasi-stationary solution that is periodic in $t$. To resolve the boundary layers of order $\varepsilon$ numerically we use the hp-adaptive strategy of Schwab and Suri [22], where we use a mesh with curved cells with a high aspect ratio (see the meshes in Fig. 2) where the size normal to the boundary behaves linear in $\varepsilon$ (or $\sqrt{\nu}$), $\sqrt{\omega}$ and the polynomial order $p$. In the experiments we have used a uniform polynomial order $p=8$ and a time steps $\Delta t$ between $2\cdot 10^{-4}$ for $\varepsilon=10^{-2}$ and $2\cdot 10^{-3}$ for $\varepsilon=10^{-1}$. Even so not necessary, we use the same mesh and polynomial order for the approximative models of order $N=0,1,2$. Note, that the computation of the reference solution is by far more expensive than the computation of the approximative models.

In Fig. 2 the approximative pressures distributions

[TABLE]

that are composed of the modes $p^{\varepsilon,N}_{k}$ are shown for $N=0,1,2$ in comparison with the reference solution that is obtained in time domain. Depending on the magnitude of the viscosity and the source a good agreement is achieved with a high enough order of the approximative solution that is $N=0$ or $N=2$ in the two examples.

Using the inverse Fourier transform of the reference pressure in the respective last period we have computed the $L^{2}(\Omega)$-norm of the contributions to the frequencies [math], $\omega$, $2\,\omega$ and $3\,\omega$ and obtained a very well agreement with the contributions of the approximative solution of second order (see Fig. 3).

Finally, we have studied the modelling error of the approximative models of order $N=0,1,2$ in dependence of the parameter $\varepsilon$ in the $L^{2}(\Omega)$-norm in the last period of the in time-domain computed reference solution (see Fig. 4). We clearly observe convergence orders of $1$ and $2$ of the relative $L^{2}$-modelling error for the approximative model of order [math] and $1$, which are both linear and have only contributions at the frequency $\omega$. For the approximative model of order $2$ that has non-linear contributions for the frequencies [math] and $2\cdot\omega$ we find numerically a convergence rate larger than two and much lower error levels for the considered values of $\varepsilon$.

5 Conclusion

For the acoustic wave-propagation in the presence of viscous boundary layers and frequency interaction due to nonlinear advection terms approximative models up to order 2 with impedance boundary conditions have been introduced. They are based on a multiscale expansion and multiharmonic analysis of the compressible nonlinear Navier-Stokes equations for small sound amplitudes of $O(\varepsilon^{2})$ and viscosities of $O(\varepsilon^{2})$, where $\varepsilon$ is a small parameter. In the approximative models the contributions to the excitation frequency $\omega$ and its harmonics can be computed sequentially after each other. As the approximative models are for macroscopic pressure or velocity fields no adaptive mesh refinement is necessary, which implies in the original model of the Navier-Stokes equations a reduction of time steps in explicit and semi-implicit schemes. In numerical experiments using finite element discretisation of the frequency domain approximative model and the original instationary compressible Navier-Stokes equations a good agreement has been shown as well as a convergence of the approximative solutions to the reference solution.

The derivation is for small sound amplitudes and it is of interest to extent the results to higher sound amplitudes than $O(\varepsilon^{2})$ which exhibit higher frequency interaction in the viscous layers. Moreover, the nonlinear Navier-Stokes equations and frequency interaction is of high interest for the modelling of liners where for periodically perforated plates the method of surface homogenization has been developed to obtain approximative impedance transmission conditions [6, 8].

Acklowledgements

The authors gratefully acknowledge the financial support by the research center Matheon through the Einstein Center for Mathematics Berlin (project MI–2) and are thankful to the fruitful exchange with the DLR Berlin.

Appendix A

A.1 Deriving the far field equations

The far field terms will be defined in physical coordinates in the whole domain $\Omega$ where we assume $\partial\Omega$ to be $C^{\infty}$. Inserting the expansion (2.8) into the system (2.5) for ${\mathbf{v}}^{\varepsilon,M}$, $p^{\varepsilon,M}$ for a particular coordinate ${\mathbf{x}}\in\Omega$ and letting $\varepsilon$ tend to zero, the near field terms concentrate closer and closer to the wall and vanish on ${\mathbf{x}}$. Collecting terms of the same order in $\varepsilon$ results in the far field equations:

[TABLE]

for $k\in{\mathbb{N}}^{+}$, where ${\mathbf{v}}^{-1}_{k}={\mathbf{v}}^{-2}_{k}={\mathbf{0}}$, and ${{\mathbf{f}}_{j}}=0$ for $j>1$. The far field equations will be completed by boundary conditions, which are specified in Sec. 2 for $j=0,1,2$.

A.2 Deriving the near field equations

The following near field equations in local coordinates derived under a condition, that the near field expansion (2.9) inserted into (2.5) leaves a residual as small as possible in powers of $\varepsilon$, which is at least of order $\varepsilon^{N+1}$

[TABLE]

where the coefficients are following

[TABLE]

Coefficients related to the nonlinear terms in the momentum equation given by

[TABLE]

A.3 Impedance boundary conditions for the far field with the source on the

boundary

In case if the source function ${\mathbf{f}}$ does not disappear on the boundary impedance boundary conditions contain additional terms in the frequency of the excitation mode $1\cdot\omega$. For the far field pressure they are

[TABLE]

and for the far field velocity

[TABLE]

The systems for frequencies $0\cdot\omega$ or $2\cdot\omega$ do not change.

A.4 Reducing the order of derivation in the system for frequency $2\cdot\omega$

Using vector calculus identities and the relation between far field pressure and velocity at frequency $1\cdot\omega$ we can rewrite equation (2.16a) as

[TABLE]

where ${\mathbf{J}}$ denotes the Jacobi matrix. Here, only second derivatives of $p^{\varepsilon,2}_{1}$ appear, but no third derivatives.