Equation level matching: An extension of the method of matched asymptotic expansion for problems of wave propagation

TL;DR

This paper proposes an equation-level matching method as an alternative to traditional matched asymptotic expansions, simplifying the process and extending applicability to complex wave propagation problems.

Contribution

It introduces a novel equation-level matching approach that avoids solving intermediate asymptotic equations, enabling analysis of wave problems where traditional methods fail.

Findings

Allows matching without explicit solutions of asymptotic equations

Effective in wave problems with differing time behaviors

Simplifies the process of deriving approximate solutions

Abstract

We introduce an alternative to the method of matched asymptotic expansions. In the "traditional" implementation, approximate solutions, valid in different (but overlapping) regions are matched by using "intermediate" variables. Here we propose to match at the level of the equations involved, via a "uniform expansion" whose equations enfold those of the approximations to be matched. This has the advantage that one does not need to explicitly solve the asymptotic equations to do the matching, which can be quite impossible for some problems. In addition, it allows matching to proceed in certain wave situations where the traditional approach fails because the time behaviors differ (e.g., one of the expansions does not include dissipation). On the other hand, this approach does not provide the fairly explicit approximations resulting from standard matching. In fact, this is not even its aim, which to produce the "simplest" set of equations that capture the behavior.

Equations

$$\left[ϵ\frac{d^2}{d\text{/}{x}^2}-\frac{d}{d\text{/}x}+1+ϵ\right]u=0\text{/},0<x<1\text{/},0<ϵ\ll 1\text{/},$$

$$u\sim\mbox{$\sum_{n=0}^{\infty}$}\,\epsilon^{n}\,\tilde{u}_{n}(x)\/,\quad\mbox{with}\quad\tilde{\mathcal{L}}\tilde{u}_{n}=[d_{x}^{2}+1]\tilde{u}_{n-1}\/,$$

$$\left[\frac{1}{ϵ}\frac{d^2}{d\text{/}{y}^2}-\frac{1}{ϵ}\frac{d}{d\text{/}y}+1+ϵ\right]u=0\text{/}.$$

$$u\sim\mbox{$\sum_{n=0}^{\infty}$}\,\epsilon^{n}\,\hat{u}_{n}(y)\/,\quad\mbox{with}\quad\hat{\mathcal{L}}\hat{u}_{n}=-\hat{u}_{n-1}-\hat{u}_{n-2}\/,$$

$${\tilde{u}}_0=e^x\text{/},{\tilde{u}}_1=2x{e}^x\text{/},{\tilde{u}}_2=(2x^2+4x)e^x\text{/},\dots$$

$${\hat{u}}_0=a_0{v}_0\text{/},{\hat{u}}_1=a_0{v}_1+a_1{v}_0\text{/},{\hat{u}}_2=a_0{v}_2+(a_0+a_1)v_1+a_2{v}_0\text{/},\dots$$

$$v_0=1-e^y\text{/},v_1=y+y{e}^y\text{/},v_2=y+\frac{1}{2}{y}^2+(y-\frac{1}{2}{y}^2)e^y\text{/},\dots$$

$$u\sim\mbox{$\sum_{n=0}^{\infty}$}\,\epsilon^{n}\,u_{n}(x;\epsilon)\/,\quad\mbox{with}\quad{\mathcal{L}}\,u_{n}=-u_{n-1},$$

$$d_x^2{u}_{n,j-1}-d_x{u}_{n,j}+u_{n,j}=-u_{n-1,j}\text{/},$$

$$d_y^2{u}_{n,j}-d_y{u}_{n,j}+u_{n,j-1}=-u_{n-1,j-1}\text{/}.$$

$$u=(e^{{\lambda }_{1}x}-e^{{\lambda }_2(x-1)+{\lambda }_1})/(1-e^{{\lambda }_1-{\lambda }_2})\text{/},$$

$$u=e^x-e^{(x-1)/ϵ}+\mathcal{O}(ϵ).$$

$$u_0=(e^{{\mu }_{1}x}-e^{{\mu }_2(x-1)+{\mu }_1})/(1-e^{{\mu }_1-{\mu }_2})\text{/},$$

$${\vec{u}}_t+(\vec{F}(\vec{u}){)}_x$$

$${\lambda }_t+\varphi \left(\vec{u}\right){\lambda }_x$$

$$\vec{F}={\vec{F}}_0+A_1\vec{u}+{\vec{A}}_2(\vec{u}\text{/},\vec{u})+{\vec{A}}_3(\vec{u}\text{/},\vec{u}\text{/},\vec{u})+\dots \text{/},A_1=\nabla (\vec{F})(\vec{0})\text{/},\text{vspace}\ast -0.3em$$

$$\vec{u}\sim ϵ{\vec{u}}_1(\chi \text{/},\tau )+{ϵ}^2{\vec{u}}_2(\chi \text{/},\tau )+\dots \text{mbox}and\lambda \sim {\lambda }_0(\chi \text{/},\tau )+ϵ{\lambda }_1(\chi \text{/},\tau )+\dots \text{/},$$

$$\mathcal{O}({ϵ}^1)(A_1-c)({\vec{u}}_1{)}_{\chi }$$

$$\mathcal{O}({ϵ}^2)(A_1-c)({\vec{u}}_2{)}_{\chi }$$

$$\dots$$

$$\mathcal{O}({ϵ}^n)(A_1-c)({\vec{u}}_n{)}_{\chi }$$

$${\vec{u}}_1={\sigma }_1\vec{r}\text{/},$$

$$({\sigma }_1{)}_{\tau }+(a{\sigma }_1^2{)}_{\chi }=\vec{l}\cdot \vec{W}\left({\sigma }_1\vec{r}\text{/},{\lambda }_0\right)\text{/},$$

$$\mathcal{O}({ϵ}^0)({\varphi }_0-c)({\lambda }_0{)}_{\chi }$$

$$\dots$$

$$\mathcal{O}({ϵ}^n)({\varphi }_0-c)({\lambda }_n{)}_{\chi }$$

$$\mathcal{L}\vec{\Delta }=\vec{0}\text{/},$$

$${\vec{u}}_T+(\vec{F}(\vec{u}){)}_X$$

$${\lambda }_T+\varphi (\vec{u}){\lambda }_X$$

$$\vec{u}$$

$$\text{mbox}and\lambda$$

$$\mathcal{O}({ϵ}^1)(A_1-c)({\vec{u}}_1{)}_{\xi }$$

$$\mathcal{O}({ϵ}^2)(A_1-c)({\vec{u}}_2{)}_{\xi }$$

$$\dots$$

$$\mathcal{O}({ϵ}^n)(A_1-c)({\vec{u}}_n{)}_{\xi }$$

$${\vec{u}}_1={\sigma }_1\vec{r}\text{mbox}and{\vec{u}}_n={\sigma }_n\vec{r}+{\vec{v}}_n\text{mbox}forn\ge 2\text{/},$$

Full text

1 Introduction

The roots of the method of matched asymptotic expansion can be traced back to the seminal work by Prandtl on viscous boundary layers [1]. With a revolutionizing idea, Prandtl realized that the effects of viscosity on an object moving through a fluid are felt only in a thin layer near the object, where the fluid velocity must rapidly match that of the object’s surface. In Prandtl’s physical picture the flow encompasses two distinct domains: an outer region, away from the object, where viscosity is negligible, and a thin transition layer, where the fluid quickly changes its velocity to match that of the moving object. A few years later G. I. Taylor elucidated some puzzles surrounding shock waves by employing a modification of Prandtl’s boundary layer concept — an internal, moving, thin layer [2] where the fluid rapidly change values.

The method of of matched asymptotic expansions was given a solid mathematical foundation in the 1950’s, through the work of Kaplun and Lagerstrom [3, 4, 5, 6]. Kaplun introduced a precise definition for the “matching” procedure for asymptotic approximations with overlapping domain of validity. For a review of his ideas see [7], as well as any of the many textbooks in asymptotic expansions that describe it, such as [8, 9, 10, 11, 12, 13, 14]. A rough, quick description, of the method follows. As an example, assume that the objective is to obtain the solution to a boundary value problem for an ODE, or a PDE, which has a small parameter in it. Then, first one obtains several expansions for the solution of the problem, under various scalings for the independent variables. Each expansion is valid in some region of space, and may have parameters that need to be determined (e.g., to satisfy the boundary conditions). If the solutions have overlapping regions of validity as the small parameter vanishes, then they can be “matched” by re-expanding each of them in terms of some intermediate scaling, and equating the result term by term. This provides relationships between the parameters that allow their determination, thus arriving at a complete description of the solution. A final step, not always taken, is to combine the various expressions for the solution (each valid in some region) to produce a “composite” or “uniform” solution, which describes the solution everywhere. In some sense, the composite solution is the “simplest” approximation that captures the whole behavior of the solution. The aim of the method introduced in this paper can be said to be: provide the “simplest” equation, or set, that capture the whole behavior of the solutions. The two concepts are related, but they are not equal, as we will see though examples.

The method of matched asymptotic expansions is very powerful, and has been successfully used to solve many problems in applications. For example: the theory of high activation energy asymptotics for flames [15], viscous flows past solid objects [16, 17, 13, 18, 19], critical layers in parallel shear flows [20], transonic flows [21], freezing/melting interfaces and (more generally) heat transfer problems [22, 23], physics of plasmas [24], electro-chemistry and electro-osmosis [25, 26], etc. However, it has limitations, particularly for problems that involve wave motion. The matching process requires detailed information about the solution, which may not be accessible. It is also, mainly, focused on “one” solution, not whole classes of possible time evolutions resulting from varying initial conditions — though it is very well suited for the calculation of eigenmodes. For problems where the time evolution is affected by, say, dissipation that occurs mainly on a boundary layer [but not everywhere], the time evolution that the various expansions yield is different for different regions — hence a “matching” as described above does not seem possible (for an example, see § 4, in particular § 4.1).

To ameliorate the difficulties outlined above, we propose an extension of the method of matched asymptotic expansions, which relies on matching by using equations, not solutions: i.e., equation level matching. To explain the idea, imagine that two [or more] different expansions (corresponding to different scalings) have been produced, each valid in some region. Each expansion is characterized by the sequence of equations that determine the terms in the expansion. In the approach proposed here, instead of solving these equations, another expansion is constructed (the uniform equations expansion), which must have the following property:

[TABLE]

The process is best illustrated with examples, starting with a very simple one § 2. The following points will become clear throughout this paper:

The “output” of this method are equations whose solutions provide uniform approximations to the problem solution. These uniform approximations are similar to the “uniform (or composite) solutions” [8, 11, 12] that the standard method produces [at least for problems where both approaches work], but not exactly the same (see the example in § 2). 2. 2.

The uniform equation expansion is not a “standard” expansion, in the sense that the small parameter appears in the equations that characterize each level of the approximation. This must be so for the expansion to actually be uniform. This makes the approach more complicated, but it gives it the added flexibility needed to deal with time dependence in wave problems. 3. 3.

The method is not a replacement for the standard method. For one, in problems where both approaches work, the standard method usually is simpler. In addition, its objective is to show that a particular set of equations provides a uniform expansion. Hence, for wave problems at least, it is a way to obtain canonical equations that model the simplest setting where a particular set of phenomena matter. But it does not actually provide any solutions. Finally, while designed to deal with problems where the standard approach has difficulties, there may very well be situations where the opposite is true — though we have not investigated this possibility, we suspect that it happens.

In the following sections we first present the method through a simple example § 2, and then show its application to more complex problems of wave propagation (a weakly nonlinear detonation model in § 3, and boundary layer dissipation for acoustic waves in § 4). We make no attempt at presenting a “general formulation” of the method for nonlinear problems (as we do for linear problems), but a specific example is discussed.

2 Asymptotic matching of equations: ODE simple BVP

We illustrate the method with an ordinary differential equation (ODE) boundary value problem (BVP) whose solution has a boundary layer. That is

[TABLE]

with boundary conditions (BC) $u(0)=1\/$ and $u(1)=0\/$. As $\epsilon\downarrow 0\/$, $u\/$ develops a steep layer 111 We skip many details here, as we assume that the readers are familiar with the standard theory of matching. near $x=1\/$. A very simple problem, meant as an illustrative example only.

2.1 Standard matching

The regular or outer expansion (can only satisfy the left BC) is

[TABLE]

where $d_{x}=\frac{d}{d\/x}$, $\tilde{\mathcal{L}}=d_{x}-1\/$, and $\tilde{u}_{n}=0\/$ for $n<0\/$. In terms of the “inner” variable $y=(x-1)/\epsilon$, the equation becomes

[TABLE]

The inner or singular expansion is then

[TABLE]

where $d_{y}=\frac{d}{d\/y}$, $\hat{\mathcal{L}}=d_{y}^{2}-d_{y}\/$, and $\hat{u}_{n}=0\/$ for $n<0\/$.

The expansion in (3) is valid everywhere, except for too close to $x=1\/$ (i.e.: as long as $d_{x}\,u=\mathcal{O}(1)\/$). The expansion in (5) is valid as long as $y\/$ is not too large. Both expansions have an overlap region, which can be captured by the variable $z=(x-1)/\delta\/$, where $\epsilon\ll\delta\ll 1\/$ — e.g., $\delta=\sqrt{\epsilon}\/$. The standard matched asymptotic expansions procedure is to write both expansions in terms of $z\/$, re-expand in terms of the new small parameters, and require that they agree order by order — thus figuring out the values of any free constants that may appear in the expansions. Here (3) yields

[TABLE]

where $u(0)=1\/$ has been enforced. Similarly, upon enforcing $u(1)=0\/$, (5) yields

[TABLE]

where the $a_{n}\/$ are constants and

[TABLE]

The matching procedure then yields $a_{0}=e\/,a_{1}=2\,e\/,a_{2}=6\,e\/$, etc.

2.2 Equation matching

The matching procedure in § 2.1 relies on being able to solve the expansion’s equations — or, at least, have detailed information about their solutions. What if this is not possible, as can be the case for complex, time-dependent problems (e.g., wave propagation)? Let us pretend this is the case here. Then, instead of solve and match, we propose to construct another expansion, including both the inner and outer problems as limits. Hence we will operate at the equation level only. The matching expansion corresponds to the composite or uniform solution in standard matching, hence we call it the uniform expansion.

In the simple case in this section the uniform expansion is

[TABLE]

where ${\mathcal{L}}=\epsilon\,d_{x}^{2}-d_{x}+1\/$, $u_{n}=0\/$ for $n<0\/$, $u_{0}(0)=1\/$, and $u_{0}(1)=u_{n}(0)=u_{n}(1)=0\/$ for $n\geq 1\/$. The operator ${\mathcal{L}}\/$ follows by looking at the terms involving $\tilde{u}_{n}\/$ in (3), those involving $\hat{u}_{n}$ in (5), and writing the operator that involves all of them. This process is the equation level analog of the standard matching composite solution, which is obtained by adding the inner and outer solutions, and subtracting the common terms. However, analog does not mean equivalent — see § 2.3

Next we check that property (1) is satisfied. First, expand each $u_{n}$ in (9) using the outer scaling $u_{n}\sim\sum_{j=0}^{\infty}\,\epsilon^{j}\,u_{n,j}(x)\/$. This yields

[TABLE]

where $u_{n,j}=0\/$ if any subscript is negative. Define $U_{m}=\sum_{n+j=m}\,u_{n,j}\/$. Then, from (10) it follows that $\{U_{m}\}\/$ satisfies (3). Next expand each term in (9) using the scaling for (5): $u_{n}\sim\sum_{j=0}^{\infty}\,\epsilon^{j}\,u_{n,j}(y)\/$. This yields

[TABLE]

Then, if $U_{m}=\sum_{n+j=m}\,u_{n,j}\/$, $\{U_{m}\}\/$ satisfies (5). Thus (9) contains both the inner and outer expansions.

For this example the leading order, ${\mathcal{L}}\,u_{0}=(\epsilon\,d_{x}^{2}-d_{x}+1)\,u_{0}=0\/$ is not substantially simpler than the original problem (2). This is not surprising when starting from a simple equation; for then there is nothing fundamentally simpler that can approximate the full behavior. One may then ask the question: is there any advantage for the approach in § 2.2 versus the one in § 2.1? If the aim is to obtain explicit approximations, certainly not, at least for examples as simple as this. However, in terms of uniform approximations, there are advantages:

The process by which composite solutions are obtained with the approach in § 2.1 is not entirely simple — particularly at higher orders. On the other hand, the approach in § 2.2 yields uniform approximations directly. But this ignores the fact that solving the equations level matching equations is, generally, harder than solving the ones from standard matching. 2. 2.

The uniform approximations produced by the equation level matching approach tend to be “better”, in the sense that they are not only more accurate, but remain valid for a larger range in the small parameter. This is discussed in § 2.3 for the example here, and later in § 4 for the acoustics example.

Of course, in the context of an example as simple as this it is hard to make meaningful comparisons. For this example writing the exact solution handily beats both techniques. The meaningful differences arise for wave problems, when the target is not a specific solution, but a simplified model equation for the physics.

2.3 Comparison between methods

Since an analytical solution to (2) exists, it is useful to compare the approximations and their errors. First: the exact solution to (2) is

[TABLE]

where $\lambda_{j}\/$ solves $\epsilon\,\lambda^{2}-\lambda+(1+\epsilon)=0\/$. Note that $\lambda_{2}=1/\epsilon-\lambda_{1}\/$, where $\lambda_{1}=1+2\,\epsilon+4\,\epsilon^{2}\/+\mathcal{O}(\epsilon^{3})\/$ for $\epsilon\ll 1\/$. Second, the standard matched asymptotic composite solution is, to leading order

[TABLE]

Finally, the leading order for the equation level matching gives

[TABLE]

where $\mu_{j}\/$ solves $\epsilon\,\mu^{2}-\mu+1=0\/$. It is easy to show that the characteristic values $\mu\/$ and $\lambda\/$ are related by $\lambda_{1}=\mu_{1}+\nu\/$ and $\lambda_{2}=\mu_{2}-\nu\/$, where $\nu=\epsilon+2\epsilon^{2}+\dots$. Hence the characteristic values for ${\mathcal{L}}\/$ (the $\mu_{j}\/$) approximate those of the full problem up to errors that are small — this is why $u_{0}\/$ can provide a uniform approximation.

Note that

(14) is substantially more complicated than (13), but not alarmingly so. The extra complication is compensated by increased accuracy, and (more important) qualitative validity even for not small $\epsilon$-values. See figures 1–2. 2. 2.

The higher order terms grow rapidly in complication. However, if the objective is to capture the essential behavior at leading order (often the case in wave problems), higher order terms are not important.

A comparison between the exact solution (12), and (14) is show in figure 1 for two different values of $\epsilon\/$. An important feature is that, besides being a uniformly valid approximation for $\epsilon\/$ small, it behaves qualitatively correct even for large values of $\epsilon\/$. Validity for a wide $\epsilon$-range is an attractive property for complex physical situations, where the aim is to obtain as simple a model as possible, for the purpose of understanding the observed behaviors — often the case in wave research.

Figure 2 shows that equation level matching leads to errors which are both smaller, and remain small for large $\epsilon\/$. Note: that the error remains small even for huge $\epsilon$-values (it keeps decreasing beyond $\epsilon=1\/$, and asymptotes $\cot(1)-\arccos(\sin(1))\approx 0.071\/$); this is a consequence of the simplicity of the example, and should not be taken seriously

3 Example from waves theory: weakly nonlinear

detonations

In this example we use the method to study weakly nonlinear waves in a reactive gas. The goal is to rationally treat dissipation as a singular perturbation problem, and derive the appropriate model incorporating dissipative effects. For reactive gas dynamics background information see, e.g., [27, 28, 29]. The main technical issue to address is: in the limit of high activation energy and weak heat release, kinetic theory predicts exponentially small diffusion coefficients. 222 For inert gases it is possible to consider length scales where dissipative effects can be incorporated via a weakly nonlinear wave expansion — e.g., derivation of Burgers’ equation. But chemical reactions introduce a spacial scale, the reaction length, which is much larger than the scale where dissipation plays a role. Hence in weakly nonlinear wave perturbation theory these effects are left out [30, 31, 32]. However, there is some current interest on the influence of transport effects (viscosity, thermal conductivity, and species diffusion) on the stability of detonations [33, 34, 35]. Thus an investigation of these effects within a weakly nonlinear model seems appropriate.

3.1 The mathematical model

Here we will consider an abstract version of the 1-D reacting equations for compressible gas dynamics. Specifically, in adimensional variables:

[TABLE]

where $\vec{u}\/$ is the vector of perturbations (from a constant rest state) to the “fluid” conserved densities, $\vec{F}\/$ is the vector of fluxes, $0\leq\lambda\leq 1\/$ is a reaction progress variable, $\phi\/$ is the “particle speed”, $0<\epsilon\ll 1\/$ is the inverse of the “activation energy”, and $\delta>0\/$ is a parameter characterizing the transport effects — assumed exponentially small in $\epsilon\/$. $\vec{W}\/$ is the vector of sources caused by the “reaction” (the prefactor $\epsilon^{2}\/$ indicates “weak heat release”), while $w\geq 0\/$ is the reaction rate. Finally, $\mathcal{D}(\partial_{x})=D\,\partial_{x}^{2}\/$ ($D\/$ a constant square matrix) is a diffusion operator: all the solutions to the linearized ($\vec{u}\/$ infinitesimal), non-reacting (no $\lambda\/$), problem decay in time.

We assume that $\vec{F}\/$ is smooth, and that the square matrix $\nabla(\vec{F})\/$ has distinct real eigenvalues (the system is strictly hyperbolic). Expanding for $\vec{u}\/$ small

[TABLE]

where $\vec{A}_{2}$ is a symmetric vector valued bilinear function, etc.

The motivation for the assumed scalings is: we seek for a situation where the waves are weakly nonlinear. Because of the reaction terms, such waves are possible only for a suitably small heat release. Further, for the weak waves to couple with the reaction, sensitive dependence of the reaction terms on the fluid variables is needed, as provided by the high activation energy assumption.

3.2 Weakly nonlinear detonation waves expansion

Let $c\/$ be an eigenvalue of $A_{1}\/$ such that $c\neq\phi_{0}=\phi(\vec{0})\/$. Let $\vec{l}\/$ and $\vec{r}\/$ be left and right eigenvectors of $A_{1}\/$ corresponding to $c\/$, normalized by $\vec{l}\cdot\vec{r}=1\/$. Then we propose a weakly nonlinear traveling wave expansion:

[TABLE]

where $\chi=x-c\,t\/$ and $\tau=\epsilon\,t\/$. Assume that the wave moves into the rest state ($\vec{u}\/$ vanishes for $\chi\to\infty$). Substituting into (15), and collecting equal powers of $\epsilon$,

[TABLE]

In (21):   (i) $n\geq 2\/$;   (ii) $\vec{N}_{n}=\sum_{j\geq 2}\;\;\sum_{n_{1}+\dots+n_{j}=n}\vec{A}_{j}\left(\vec{u}_{n_{1}}\/,\,\dots\vec{u}_{n_{j}}\right)\/$,  with $n_{q}\geq 1\/$;   (iii) $\vec{W}_{q}\/$ is the $q$-th term in the expansion of $\vec{W}\/$ — note that $\vec{W}_{q}$ depends on $\vec{u}_{1}\/,\,\dots\vec{u}_{q+1}\/$ and $\lambda_{0}\/,\,\dots\lambda_{q}\/$ only. Because of the assumption on the size of $\delta\/$, no dissipative terms appears at any order.

The $\mathcal{O}(\epsilon)\/$ equations yield

[TABLE]

where $\sigma_{1}=\sigma_{1}(\chi\/,\,\tau)\/$ is scalar valued. For the $\mathcal{O}(\epsilon^{2})$ equations we write $\vec{u}_{2}=\sigma_{2}\,\vec{r}+\vec{v}_{2}$, where $\vec{v}_{2}$ is a linear combination of the right eigenvectors of $\vec{A}_{1}$ corresponding to eigenvalues different from $c\/$ (thus $\vec{l}\cdot\vec{v}_{2}=0\/$). Then (20) has a solution if and only if the right hand side is orthogonal to $\vec{l}$. Thus

[TABLE]

where $a=\vec{l}\cdot\vec{A}_{2}(\vec{r}\/,\,\vec{r})\/$. In general, at any order $\vec{u}_{n}=\sigma_{n}\,\vec{r}+\vec{v}_{n}\/$, with $\vec{l}\cdot\vec{v}_{n}=0\/$. Then $\vec{v}_{n}\/$ follows from (21), and $\sigma_{n}\/$ from the solvability condition at the next order — the right had side of (21) must be orthogonal to $\vec{l}\/$ at all orders (see § 3.2.1).

Finally the expansion for $\lambda\/$ is governed by

[TABLE]

where the BC $\lambda=0\/$ is imposed for $\chi\to\infty\/$ (ahead of the wave). In (25) $\phi_{q}\/$ and $w_{q}\/$ are the $q$-th terms in the expansions for $\phi\/$ and $w\/$. Note that $w_{n}\/$ involves dependence on $\vec{u}_{n+1}\/$ and $\lambda_{n}\/$, via the terms $(\vec{u}_{n+1}\cdot\nabla_{u})\,w(\vec{u}_{1}\/,\,\lambda_{0})$ and $\lambda_{n}\,w_{\lambda}(\vec{u}_{1}\/,\,\lambda_{0})\/$.

3.2.1 The higher orders

The leading order equations controlling the expansion are (23) and (24). Next we describe the higher order equations. Substitute $\sigma_{1}\to\sigma_{1}+\Delta\sigma\/$ and $\lambda_{0}\to\lambda_{0}+\Delta\lambda\/$ into (23) and (24), where $\Delta\sigma\/$ and $\Delta\lambda\/$ are infinitesimals. Then $\Delta\sigma\/$ and $\Delta\lambda\/$ satisfy linear homogeneous equations (with coefficients that depend on $\sigma_{1}\/$ and $\lambda_{0}\/$)

[TABLE]

where $\vec{\Delta}\/$ is the vector with components $\Delta\sigma\/$ and $\Delta\lambda\/$. The higher order equations are forced versions of this equation: ${\mathcal{L}}\,\vec{\Delta}_{n}=\vec{f}_{n}\/$, $n\geq 2\/$, where $\vec{f}_{n}\/$ depends on the lower order terms only, and $\vec{\Delta}_{n}$ is the vector with components $\sigma_{n}\/$ and $\lambda_{n-1}\/$. In specific examples one can see that these equations can develop secularities for $\chi\/$ or $\tau\/$ large. The expansion in this section requires $\chi=o(1/\epsilon)\/$ and $\tau=o(1/\epsilon)\/$.

3.3 Burgers’ shock waves expansion

Change variables in (15–16), to $T=t/\delta\/$ and $X=x/\delta\/$. Then

[TABLE]

Expand

[TABLE]

where $\xi=X-c\,T\/$ and $\mu=\epsilon\,T\/$. Then

[TABLE]

The solution to these equations follows the same pattern as for (19–21). Write

[TABLE]

with $\vec{l}\cdot\vec{v}_{n}=0\/$. Then the $\sigma_{n}\/$ are determined by the solvability conditions: the right hand sides in the equations above must be orthogonal to $\vec{l}\/$. In particular:

[TABLE]

where $a\/$ is as in (23) and $\nu=\vec{l}\cdot D\,\vec{r}>0\/$ is a constant — $\nu>0\/$ follows because $D\/$ is a diffusion matrix. Note that (33) is Burgers’ equation, the canonical equation describing viscous weak shocks [36, 37, 38, 39, 40].

Similarly

[TABLE]

These equations mean that no reaction occurs in this limit.

3.4 Equation level matching

Next we extend the process introduced in § 2 to this section’s problem. Consider:

• (a)

The expansion in § 3.2, given by equations (19–21) and (24–25).

• (b)

The expansion in § 3.3, given by equations (29–31) and (34–35).

Then take the “union” of these expansions to produce a unified expansion reducing to those in § 3.2 and § 3.3 in the appropriate set of variables. Thus matching occurs at the equation level. Hence, because the equations match, we know that the solutions also match (in the sense of standard matched asymptotic expansions), without the need to actually compute (or even know) the solutions.

We propose the unified expansion

[TABLE]

where $\vec{l}\cdot\vec{v}_{n}=0\/$. Here, unlike (18) or (28), the terms include an explicit dependence on the small parameter $\delta\/$. A term by term formal expansion is not possible. The equations must be obtained via matching, as follows:

[TABLE]

where we use the same notation as in § 3.2 and § 3.3. Similarly

[TABLE]

This expansion matches the one in § 3.2 because 333 Recall that $\delta$ is smaller than any power of $\epsilon\/$. the effect of the terms multiplied by $\delta\/$ in the equations occurs beyond all the $\vec{u}_{n}\/$ and $\lambda_{n}\/$, so that the equations reduce to those in (19–21) and (24–25). Rewriting the equations in terms of the variables $\xi=\chi/\delta\/$ and $\mu=\tau/\delta\/$, it is easy to see (same argument) that the expansion here matches the one in § 3.3 — the equations reduce to those in (29–31) and (34–35).

Finally, note that the leading order uniform expansion is given by

[TABLE]

The equations in [41] reduce to the ones above in the plane-wave case. However, in [41] the non-physical assumption $\delta\sim\epsilon$ was introduced. A similar approach to incorporate transport effects in a more elaborate setting (including 2-D effects, species diffusion and decoupling of the thermal effects) can in be found in [42].

The dissipative term in (43–44) adds structure to the shocks, which in the inviscid case are just point discontinuities moving in space-time. Although this is only of secondary importance when studying the steady traveling wave solutions to (43–44), the finite shock width appears to play a more important role in the highly unstable regime where detonations are typically found [33, 34]. More importantly, including dissipation in a modified version of the weakly nonlinear theory presented in this paper can trigger new types of subsonic traveling wave solution not present in the inviscid theory [43]. Further physically relevant questions occur in relationship with the role of dissipation for multi-dimensional wave interactions. This is the subject of current study by the authors, and will be reported elsewhere.

4 Example: dissipation from the acoustic boundary layer

One topic we hope this method is useful for is: incorporate into wave models dissipation and/or dispersion, when these effects occur at higher order in standard asymptotic wave theory. In these situations conventional matching can fail as well because, for example: (i) The behavior is too complex, and there are no solutions available for the matching process. (ii) The time behavior for the “inner” and “outer” expansions differs. In this section we study a simple example of the situation in (ii), and consider the effects of boundary layer dissipation in acoustics. In § 4.1 we include a brief description of the standard methods in this subject.

For simplicity, we work in 2-D, where the linearized isentropic 444 In particular this excludes the thermal boundary layer. This is meant as an illustrative example. Navier-Stokes equations for flow a channel, $0<\tilde{y}<L\/$ and $-\infty<\tilde{x}<\infty\/$ are

[TABLE]

Here $\rho_{0}=$ ambient density, $c_{0}=$ sound speed, and $\nu=$ kinematic viscosity. The tildes denote dimensional variables. Nondimensionalize: $\tilde{\rho}=\rho_{0}\,\rho\/$, $\tilde{u}=c_{0}\,u\/$, $\tilde{v}=\epsilon\,c_{0}\,v\/$, $\tilde{x}=\lambda\,x\/$, $\tilde{y}=L\,y\/$, and $\tilde{t}=(\lambda/c_{0})\,t\/$, where $\lambda\/$ is a typical wavelength and $\epsilon=L/\lambda\/$. Then

[TABLE]

where $\delta=\nu/(c_{0}\,L)\/$ and $\Delta=\partial_{y}^{2}+\epsilon^{2}\,\partial_{x}^{2}\/$. The equations apply for $-\infty<x<\infty\/$ and $0<y<1\/$, with boundary conditions: $u=v=0\/$ at both $y=0\/,\,1\/$.

4.1 Historical perspective

Here we consider expansions for (48–50) in the regime $0<\delta\ll\epsilon\ll 1\/$, where $\epsilon=L/\lambda\/$ is the long wave parameter. In fact, we will assume that $\delta=O(\epsilon^{m})=\mathcal{C}\,\epsilon^{m}\/$, where $\mathcal{C}>0\/$ is a constant and $m=2\,p+1>1\/$ is an odd integer. 555 This is to simplify the algebra. Note that for some wind instruments (e.g., flute) $p=1\/$ is reasonable.

It should be obvious that a “regular” (powers of $\epsilon\/$) or “outer” expansion for (48–50) is inviscid at leading order, with the effects of viscosity appearing only as forcing terms at higher order — in fact, causing secular growth. 666 An easy way to see why secularities arise, is to do a regular expansion for the equation $\eta_{t}+\eta_{x}=\epsilon\,\eta_{x\/x}\/$. One could, in principle, eliminate these secularities via a multiple times expansion. But this would be pointless, as the total dissipation is usually dominated by the boundary layer dissipation. On this last point: because the regular expansion can satisfy the BC for $v\/$, but not for $u\/$, a boundary layer expansion is needed to complete the picture. Unfortunately, because the outer/regular expansion does not incorporate any boundary layer effects, it has the “wrong” time dependence — lacks the decay that dissipation in the boundary layer induces. Thus the regular expansion cannot actually be made to match with the boundary layer expansion in the traditional sense — there are no intermediate scaled variables in which “inner” and “outer” solutions can be re-expanded, so that they match term by term. The reason is that the standard approach to matched asymptotics is biased towards steady state solutions, and “arbitrary” time dependences do not fit within it too well. Various solutions to these issues have been developed (e.g., see [44]). A few examples are:

Calculate the eigenmodes for the full 2-D or 3-D problem, including the boundary conditions at the walls. Then approximate the exact “dispersion relation” thus obtained using the small parameters [45, 46]. This approach goes back to the beginnings of the subject (acoustics). It precedes the development of the mathematical theory of matched asymptotic expansions [3, 4, 5, 6], or even the introduction of the concept of boundary layer by L. Prandtl. Two disadvantages are: (i) Very labor intensive. (ii) No extension to nonlinear problems. 2. 2.

Separate time, and use matched asymptotic expansions to calculate the eigenfunctions and eigenvalues. A variation of the approach in item 1. It has the advantage of going directly for the desired simplified dispersion relation, avoiding some messy calculations. But it is still restricted to the linear problem. 3. 3.

Calculate the boundary layer dissipation per unit area for a time harmonic field, space independent and with flow parallel to the boundary. Then use the result to correct the Helmholtz equation — see art. 328 in [47]. This is a clever shortcut for the process in item 1. But it is limited to single frequency waves, with no clear extension to nonlinear problems. It is also slightly inconsistent: because of the dissipation, the waves outside the layer are not exactly harmonic (unless an external forcing is applied). Thus one should actually compute the dissipation by the layer of a non-harmonic forcing. 4. 4.

Chester [48] generalized the approach in item 3, and computed the dissipation produced by an arbitrary time dependent velocity field, using a Laplace Transform approach. This is then converted into a drag per unit length along the tube, and inputted into a derivation of the governing equations by using conservation principles — assuming longitudinal dependence, only, in the solution. This approach is applicable to “not too nonlinear” situations — since then the boundary layer dissipation should still be a linear process.

The (slight generalization) of the method of matched asymptotic expansions introduced here allows us to bypass work-arounds such as ones above. It also produces a final result which applies for a larger set of parameter regimes. In particular, we expect it to capture the transition from thin-boundary layer to fully viscous, as the wave-length grows (or the tube diameter is reduced) — this is work in progress. Of course, the extended model includes (in the appropriate regime) the prior ones.

4.1.1 Acoustic outer/regular expansion

Start with (48–50) and substitute expansions of the form $\rho=\sum_{n=0}^{\infty}\epsilon^{n}\,\tilde{\rho}^{n}$, etc. Then collect equal powers of $\epsilon\/$

[TABLE]

where $\tilde{\rho}^{n}=\tilde{u}^{n}=\tilde{v}^{n}=0\/$ if $n<0\/$, and $\tilde{v}_{n}=0\/$ for $y=0\/$ or $y=1\/$ — the BC for $u\/$ cannot be satisfied. Note that:

The leading order equations are $\tilde{\rho}^{0}_{t}+\tilde{u}^{0}_{x}+\tilde{v}^{0}_{y}=\tilde{u}^{0}_{t}+\tilde{\rho}^{0}_{x}=\tilde{\rho}^{0}_{y}=0\/$. Average these equations over $y\/$, and use the BC to obtain the 1-D acoustic equations  $\tilde{\rho}^{0}_{t}+\tilde{u}^{M}_{x}=\tilde{u}^{M}_{t}+\tilde{\rho}^{0}_{x}=0\/$,  where $\tilde{u}^{M}=\int_{0}^{1}\tilde{u}^{0}d\/y\/$. 2. 2.

The higher orders in this expansion develop secular behavior in time, through resonances with the lower orders. These are related to bulk dissipation, since this expansion ignores the boundary layers. It may be possible to eliminate these secularities using multiple time scales. Here we will ignore them, as they occur on time scales much longer than those associated with the boundary layer dissipation — which we will incorporate, see § 4.1.3.

4.1.2 The acoustic boundary layer expansion

The boundary layers at $y=0\/,\,1\/$ can be treated in exactly the same way, so we only show the calculations for $y=0\/$. Start with (48–50) and change variables to $V=\epsilon\,v\/$ and $Y=(\sqrt{\epsilon/\delta})\,y\/$. Then

[TABLE]

where $\Delta_{b}=\partial_{Y}^{2}+\xi^{2}\,\partial_{x}^{2}\/$ and $\xi=\sqrt{\epsilon\,\delta}=\sqrt{\mathcal{C}}\,\epsilon^{p+1}\/$. Note that when $p\/$ is not an integer, $\xi\/$ involves fractional powers. This is not a problem, but it requires that extra terms be added to the regular expansion, to make matching possible (a well known phenomenon in matched asymptotic expansions). To gain simplicity, here we selected the relationship between $\delta\/$ and $\epsilon\/$ to avoid the effect.

Now substitute into (54) expansions of the form $\rho=\sum_{n=0}^{\infty}\epsilon^{n}\,\bar{\rho}^{\,n}\/$, etc., which leads to

[TABLE]

where $\bar{\rho}^{n}=\bar{u}^{n}=\bar{v}^{n}=0\/$ if $n<0\/$, and $\bar{u}^{n}=\bar{v}^{n}=0\/$ for $Y=0\/$.

These equations need BC for $Y\/$ large. In standard matching, these would follow from matching with (51–53). But this is impossible because the expansions’ time behaviors differ. Specifically: (51–53) do not incorporate any decay, while (55–57) do, via $\bar{u}^{n}_{Y\/Y}\/$ in (56). Though equation level matching is possible, as shown next.

4.1.3 The uniform expansion, expansion level matching

We use the process explained in § 2 to implement (1). Since the problem here is linear, modulo technical details, the procedure is the same. Re-write (48)–(50) in terms of the linear operator $\vec{\mathcal{L}}\/$ which is the “union” of the linear operators applied to the $n^{\rm th}\/$ order terms (on the left) in (51–53), (55–57), and the analog $y=1\/$. That is:

[TABLE]

where $\vec{u}=(\rho\/,\,u\/,\,v)\/$ and the $\mathcal{L}_{j}\/$ are the components of $\vec{\mathcal{L}}\/$, defined by the formulas. Then the equation level uniform expansion is $\vec{u}\sim\sum_{n=0}^{\infty}\,\epsilon^{n}\,\vec{u}^{\,n}\/$, with

[TABLE]

where $\vec{u}^{n}=0\/$ if $n<0\/$, and $u^{n}=v^{n}=0\/$ at both $y=0\/,\,1\/$. Showing that (51–53), (55–57), and the analog for the $y=1\/$ boundary, match [in the sense of (1)] via (61–63), can now be done exactly as in § 2.2. The calculations are a bit more cumbersome, but this is the only difference.

Let us now examine the leading order uniform equations

[TABLE]

where $u=v=0\/$ at both $y=0\/,\,1\/$, and we have dropped the superscript $0\/$ for notational simplicity. We can eliminate $v\/$, by averaging the first equation and using the BC

[TABLE]

Here both $\rho\/$ and $u^{M}\/$ are functions of $x\/$ and $t\/$ only. We can go a bit further, by taking the average of the second equation as well,

[TABLE]

where $[u_{y}]\/$ denotes the jump in $u_{y}\/$ from $y=0\/$ to $y=1\/$. The equations above are the equations for acoustic waves in a channel, with boundary layer dissipation incorporated. Just as with the example in § 2, it appears as if the validity of (64) extends beyond its original intended range. For example, it encompasses compressible Poiseuille Flow, given by

[TABLE]

or its multiples. Hence (64) may cover the whole range, from nearly inviscid waves, to heavily dissipated waves, to the transition where there are no more waves. We leave the investigation of these question for another publication.

Of course for (66) to be truly one-dimensional, $[u_{y}]\/$ (which incorporates the dissipation from the boundary layer) must be expressed in terms of $\rho\/$ and $u^{M}\/$. As shown next in § 4.1.4 and 4.1.5, this can be done under some restrictions.

4.1.4 Normal modes

Let $0<\mu=\sqrt{\delta/\epsilon}\ll 1\/$. Then (65) takes the form

[TABLE]

where $u=0\/$ for $y=0\/,\,1\/$. To show that these equations include the results from the “traditional” approach, we look for solutions with exponential time dependence — proportional to $e^{-\ell^{2}\,t}\/$ ($\ell\/$ is a complex number). Then the second equation in (68) reduces to an ode BVP for $u\/$, with solution:

[TABLE]

where

[TABLE]

The first equation in (68) then reduces to 777 Note that $M\/$ is not singular for $\ell_{*}=0\/$, and that $M\/$ is a function of $\ell_{*}^{2}\/$.

[TABLE]

where

[TABLE]

Hence a longitudinal dependence proportional to $e^{i\,k\,x}\/$, $k\/$ a real constant, yields the “dispersion” relation

[TABLE]

It is easy to see that, for $\ell_{*}\/$ with a large imaginary part, $\ell_{*}^{3}\,M(\ell_{*})\sim-\ell_{*}+2\,i\,\sigma\/$, where the error is exponentially small and $\sigma=\text{sign}(\text{Im}(\ell^{*}))\/$. Using this approximation in (71), and solving for $\ell\/$ in terms of $k\/$, yields

[TABLE]

where $p_{r}=(\pm 1+i\,\sigma)/\sqrt{2}\/$, $\sigma=\pm 1\/$, is any one of the four roots of $p_{r}^{4}=-1\/$. Hence

[TABLE]

Since the time evolution is via $e^{-\ell^{2}\,t}\/$, equation (72) shows that the boundary layer corrections have two effects:

(i) a small frequency correction to the waves, given by $\mu\,\sqrt{|k|/2}\/$.

(ii) Dissipation, with dissipation coefficient $\mu\,\sqrt{|k|/2}\/$.

This approximation requires $|k|\gg\mu^{2}\/$. It is not valid for very long waves, while (68) remains valid beyond the wave regime (very small $k\/$) — see (67). Note also that (71) has other modes/solutions in the regime $|k|\gg\mu^{2}\/$, which are excluded by (72). These modes are not wave-like, and decay much faster than the ones retained.

4.1.5 Reduction to a 1-D in space problem

Assume solutions such that (72) applies — i.e., enforce a very long wave frequency cut-off, and note that (72) is equivalent (up to the order displayed) to

[TABLE]

Introduce now the pseudo-differential operators (on functions of $x\/$) defined by

[TABLE]

With these definitions

[TABLE]

yields the dispersion relation (73). Note that (75) is (66) with $\mu\,[u_{y}]={\mathcal{L}}_{2}\,u^{M}-{\mathcal{L}}_{1}\,\rho\/$. At this moment it is (to us) unclear if an equation like (75) can be written that is fully equivalent to (66) — no frequency cut-off. Equations such as (75) can only be formulated for problems where the boundary conditions in $x\/$ are compatible with a Fourier Expansion, such as full line or periodic. This limits their usefulness.

A final note: Equations similar to (75) can be found in many papers and books dealing acoustic boundary layer dissipation (e.g., see [44, 48]). However, quite often, the pseudo-differential operators employed operate in time, not space. For example, convolution operators of the form ${\mathcal{L}}_{*}\,u=\int_{0}^{\infty}K(s)\,u(t-s)\,d\/s\/$, or defined via how they operate on exponential functions of the form $e^{i\,\omega\,t}\/$. There are a several fundamental problems with such definitions: (i) The situation under consideration involves decay of the solutions in time, hence any possible extension of them backwards in time will either blow up exponentially (at best) or not be even possible (dissipation problems are often ill-posed for negative time). Thus the operators so defined may have no meaning. (ii) Even if a meaning exists, the problem is no longer a problem for which an IVP makes sense, one needs to know the whole past to go forwards in time.

5 Conclusion

We presented a variation of matched asymptotic expansions where the matching occurs at the equation level — not their solutions. This has the advantage that it does not require much information about the inner and outer solutions in order to perform the matching. Of course, the method also has the obvious disadvantage that the equations obtained typically cannot be solved analytically, and thus numerical tools must be employed.

Acknowledgments

The research of R. R. Rosales was partially supported by NSF grants DMS-1318942 and DMS-1614043.