<figure><img src="image_1.png" /></figure>

HAL
open science

# Simulation of Acoustic Wave Propagation in Anisotropic Media Using Dynamic Programming Technique

Nikolai Botkin, Varvara Turova

► To cite this version:

Nikolai Botkin, Varvara Turova. Simulation of Acoustic Wave Propagation in Anisotropic Media Using Dynamic Programming Technique. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.36-51, 10.1007/978-3-662-45504-3_5 . hal-01286451

HAL Id: hal-01286451

https://inria.hal.science/hal-01286451

Submitted on 10 Mar 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simulation of acoustic wave propagation
in anisotropic media
using dynamic programming technique

Nikolai Botkin* and Varvara Turova

Technical University of Munich, Center for Mathematics
Boltzmannstr. 3, 85748 Garching b. Munich, Germany
{botkin,turova}@ma.tum.de

**Abstract.** It is known that the Hamiltonian of the eikonal equation for an anisotropic medium may be nonconvex, which excludes the application of Fermat's minimum-time principle related to minimum-time control problems. The idea proposed in this paper consists in finding a conflict control problem (differential game) whose Hamiltonian coincides with the Hamiltonian of the eikonal equation. It turns out that this is always possible due to Krasovskii's unification technique. Having such a differential game allows us to apply dynamic programming methods to computing the value function of the game, and therefore to describe the propagation of wave fronts. This method is very appropriate for the simulation of wave patterns in surface acoustic wave biosensors. Numerical computations given in this paper prove the feasibility of the approach proposed.

**Keywords:** WKB-approximation, Hamilton-Jacobi equations, viscosity solutions, differential game, unification

# 1 Introduction

The paper concerns the development of methods for modeling the propagation of acoustic waves in anisotropic media. This investigation is very important for many applications such as acoustic sensors whose operating principle is based on the excitation and detection of acoustic waves of very high frequency in piezoelectric crystals.

For anisotropic media, the WKB (Wentzel-Kramers-Brillouin) approximation yields eikonal equations whose Hamiltonians are neither convex nor concave in the impulse variable. Therefore, the well-known Fermat principle of wave propagation fails in this case. Moreover, the propagation occurs in such a way as if an antagonistic opponent aims to slow down the movement of the wave fronts. Thus, we arrive at the idea to use methods of differential games in the analysis of wave propagation. If the negative of the Hamiltonian of a differential game approximates the Hamiltonian of the eikonal equation, then the value function of

<footnote>* Corresponding author.</footnote>
the game approximates the phase function satisfying the eikonal equation. The authors have developed effective and precise algorithms for solving Hamilton-Jacobi equations arising from differential games, which yields an effective tool for the numerical investigation of eikonal equations.

If a differential game is chosen appropriately, the level sets of its value function represent the wave fronts, and optimal trajectories are associated with the propagation of rays. Thus, it makes possible to describe very complicated behavior of rays using game-theoretic classification of the so-called singular surfaces that can attract, repulse, and break the trajectories. For example, the caustic-like behavior of rays can be interpreted as the attraction of neighboring optimal trajectories to a singular surface. The monograph [1] extends the classical method of characteristics by introducing the so-called generalized characteristics that are related to the above mentioned trajectories and singular surfaces.

The main objective of this paper is the numerical simulation of propagation of bulk and surface acoustic waves in anisotropic monocrystals and multi-layered structures used in surface acoustic wave sensors. It is demonstrated that the propagation fronts can be found very precisely even in the case of very complicated geometry of wave emitters. Numerical results are presented for the case of bulk and surface waves characterized by non-convex slowness surfaces.

This investigation is inspired by the cooperation with professor A. A. Melikyan (deceased) from the Institute for Problems in Mechanics, Moscow, Russia.

## 2 Wave velocity in piezoelectric crystals and Eikonal equation

Assume that the indexes i, j, k, l run from 1 to 3 and use the summation convention over the repeated indexes. Let $u_1$, $u_2$, and $u_3$ be the displacements in $x_1$, $x_2$, and $x_3$ directions, respectively; $\varphi$ is the electric potential such that the electric field $E_i$ is given by the relation $E_i = \partial\varphi/\partial x_i$. Electro-elasticity equations for a piezoelectric anisotropic crystal read:

$$ \rho u_{itt} - C_{ijkl} \frac{\partial^2 u_l}{\partial x_j \partial x_k} + e_{kij} \frac{\partial^2 \phi}{\partial x_k \partial x_j} = 0, \quad (1) $$

$$ \epsilon_{ij} \frac{\partial^2 \phi}{\partial x_i \partial x_j} + e_{ikl} \frac{\partial^2 u_l}{\partial x_i \partial x_k} = 0. \quad (2) $$

where $\rho$, $\epsilon_{ij}$, $e_{ikl}$, and $C_{ijkl}$ denote the density, the material dielectric tensor, the stress piezoelectric tensor, and the elastic stiffness tensor, respectively.

The WKB (high frequency) approximation (see e.g. [2]) uses the ansatz

$$ u_j = u_j^0(t, x) \cdot \varepsilon e^{iS(t, x)/\varepsilon}, \quad \phi = \phi^0(t, x) \cdot \varepsilon e^{iS(t, x)/\varepsilon}, \quad (3) $$

where $\varepsilon = \omega^{-1}$ is a small parameter ($\omega$ is the frequency), $S(t, x)$ is the phase function, and $u_j^0(t, x)$ and $\phi^0(t, x)$ are functions defining the wave polarization. The symbol $i$ in the exponent denotes the imaginary unit.
Substitution of (3) into (1) and (2) and collection of the terms of order $1/\epsilon$ yield the equations

$$ \left( -\rho S_t^2 \delta_{il} + C_{ijkl} \frac{\partial S}{\partial x_j} \frac{\partial S}{\partial x_k} \right) u_l^0 - e_{kij} \frac{\partial S}{\partial x_k} \frac{\partial S}{\partial x_j} \phi^0 = 0, \quad (4) $$

$$ -e_{ikl} \frac{\partial S}{\partial x_i} \frac{\partial S}{\partial x_k} u_l^0 - \epsilon_{ij} \frac{\partial S}{\partial x_i} \frac{\partial S}{\partial x_j} \phi^0 = 0. \quad (5) $$

The condition of nontrivial solvability of the system (4) and (5) leads to the eikonal equation

$$ \det \left[ \frac{1}{\rho} \left( \frac{C_{ijkl} \frac{\partial S}{\partial x_j} \frac{\partial S}{\partial x_k} - e_{kij} \frac{\partial S}{\partial x_k} \frac{\partial S}{\partial x_j}}{-e_{ikl} \frac{\partial S}{\partial x_i} \frac{\partial S}{\partial x_k} - \epsilon_{ij} \frac{\partial S}{\partial x_i} \frac{\partial S}{\partial x_j}} \right) - \left( \frac{\delta_{il}}{0_{1 \times 3}} \frac{|0_{3 \times 1}|}{0} \right) S_t^2 \right] = 0, $$

which can be rewritten as

$$ S_t - |\nabla S| c_{\alpha} \left( \frac{\nabla S}{|\nabla S|} \right) = 0, \quad (6) $$

where $c_\alpha(n)$, $\alpha = 1, 2, 3$, are eigenvalues of the problem

$$ \det \left[ \frac{1}{\rho} \left( \frac{C_{ijkl} n_j n_k}{-e_{ikl} n_i n_k} - \frac{-e_{kij} n_k n_j}{-\epsilon_{ij} n_i n_j} \right) - \left( \frac{\delta_{il}}{0_{1 \times 3}} \frac{|0_{3 \times 1}|}{0} \right) c^2 \right] = 0. $$

Here, $n_1, n_2, n_3, |n| = 1$, are components of the normalized wave vector (the direction of propagation). Therefore, for each vector $n$, there are three types of waves propagating in this direction. Each of them has its own velocity $c_\alpha$ and the corresponding nontrivial solutions, $u_l^0$ and $\phi^0$, of (4) and (5) defining the wave polarization.

Figure 1a shows the phase velocity surface for LiTaO$_3$ piezoelectric crystals. This surface is obtained as the set of points of the form $c_\alpha(n) \cdot n$, where $n$ belongs to a grid on the surface of the unit sphere. The index $\alpha$ corresponds to a quasi shear wave where the displacements are near orthogonal to the propagation direction. Figure 1b presents the so-called slowness surface consisting of points of the form $c_\alpha^{-1}(n) \cdot n$. It is easy to prove that the slowness surface can be described as $\{p \in \mathbb{R}^3 : c_\alpha(p/|p|)|p| = 1\}$. The nonconvexity of the slowness surface shows that the Hamiltonian, $c_\alpha(p/|p|)|p|$, of equation (6) is non-convex in $p$.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Characteristic surfaces for lithium tantalate (LiTaO₃): (a) Phase velocity surface; (b) Slowness surface. The coordinate frames show the crystalline axes X, Y, and Z.</figcaption></figure>

**3 Surface acoustic wave biosensors**

Biosensors serve for the measurement of small amounts of biological substances in liquids. Usually a biosensor can be considered as a multi-layered structure (see Fig. 2) whose bottom layer is the ST-cut of piezoelectric α-quartz. Acoustic shear waves are excited here by means of a high-frequency voltage applied to electrodes placed on the ST-cut surface. The waves are transmitted into an isotropic guiding layer deposited on the top of the quartz substrate. The top gold layer is covered by DNA or RNA molecules, aptamers, that are able to specifically bind protein molecules from the contacting liquid. Binding protein molecules results in additional mass loading, which causes a phase shift in the electric signal measured by the output electrodes.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Structure of acoustic biosensor: (a) Spatial representation; (b) Vertical cross-section with thicknesses of the layers.</figcaption></figure>

A high sensitivity regarding to the added mass is achieved due to the usage of shear horizontally polarized guided waves (Love waves) because of their low
interaction with the contacting fluid. The input and output electrodes are located between the substrate and the guiding layer. To obtain purely shear polarized modes, the direction of the wave propagation is chosen to be orthogonal to the crystalline X-axis (see Fig. 3).

<figure><img src="image_5.png" /><figcaption>Fig. 3. Piezoelectric α-quartz crystal: (a) Direction of crystalline axes; (b) Orientation of ST-cut.</figcaption></figure>

The next three subsections consider a mathematical model of the biosensor and two methods of numerical investigation of acoustic Love waves including their phase velocity, decay with depth, polarization, etc.

## 3.1 Mathematical model of biosensor

The governing equations for the displacements and the electric potential in the quartz substrate are given by formulae (1) and (2).

The gold layer is conductor so that there is no electric field inside it. The electric field inside the guiding layer is also neglected because of its low dielectric permeability. Therefore, the electric potential vanishes, and the gold and guiding layers are described by the equation of the form

$$ \bar{\rho}u_{i tt} - \bar{C}_{ijkl} \frac{\partial^2 u_l}{\partial x_j \partial x_k} = 0. \qquad (7) $$

In the fluid layer, the Stokes and mass conservation equations hold:

$$ \rho_0 v_i t - \nu \Delta v_i - \left(\zeta + \frac{\nu}{3}\right) \frac{\partial}{\partial x_i} \text{div} v + \frac{\partial}{\partial x_i} \mathcal{P} = 0, \qquad (8) $$

$$ \gamma \mathcal{P}_t + \frac{\partial}{\partial x_i} v_i = 0, $$

where $v_i$ are components of the velocity, $\mathcal{P}$ is the pressure, $\rho_0$ is the fluid density at a reference pressure $P_0$, $\nu$ and $\zeta$ are the dynamic and volume viscosities of the fluid, respectively, and $\gamma = \frac{1}{\rho_0} \frac{\partial \rho}{\partial \mathcal{P}}|_{P_0}$ is the compressibility of the fluid.
A special homogenization technique developed in [3] is used to treat the aptamer-fluid structure. This bristle structure is replaced by an averaged material whose properties are derived as the number of bristles goes to infinity, their thickness tends to zero, and the height remains constant. The resulting new layer whose thickness is equal to the height of the aptamer emulates the aptamer-fluid structure. The governing equation for this layer is given by the relation (see [3] and [4])

$$ \hat{\rho}u_{i tt} - \hat{C}_{ijkl} \frac{\partial^2 u_l}{\partial x_j \partial x_k} - \hat{P}_{ijkl} \frac{\partial^2 u_{lt}}{\partial x_j \partial x_k} = 0, \quad (9) $$

where the term containing the tensor $\hat{P}$ describes the viscous damping coming from the liquid part of the aptamer-fluid structure. The term containing $\hat{C}$ represents elastic stresses. The density $\hat{\rho}$ is a weighted combination of the densities of the fluid and the aptamer. The tensors $\hat{P}$ and $\hat{C}$ are computed with FE-method using an analytical representation of solutions of the so-called cell equation arising in homogenization theory.

The conditions on the interfaces between the layers are carefully considered in [5] and [6]. Briefly, the continuity of the displacements and the equilibrium of the normal pressures must hold on the interface between every two neighboring solid layers (the averaged aptamer-fluid layer is considered as solid). Moreover, the electric displacement and the tangent component of the electric field in the substrate must be zero on the interface between the quartz substrate and the guiding layer. The conditions on the interface between the aptamer layer and the fluid include the no-slip assumption and the equilibrium of the pressures.

## 3.2 Finite element modeling

The FE-model extends the above described basic model by accounting for two alternated groups of electrodes (see Figure 2a) and a damping area around the side and bottom faces to suppress the wave reflection thereon.

The electrodes are typically made of gold. Therefore, they can be accounted for by the linear elasticity equation of type (7).

Accounting for the damping is done by adding the term $-\mathrm{div}(\beta(x)\nabla u_{it})$ to equations (2), (7), and (9), where $\beta(x)$ is a piecewise-linear function which is equal to zero outside of the damping region and grows up to some value $\beta_0 > 0$ towards the side and bottom faces.

The FE-approach provides accurate results because of accounting for the exact parameters of the sensor such as the shape of the electrodes, their position, mass, electro-conductivity properties. This allows us to estimate important characteristics of the biosensor and effects caused by scattering of waves (see [5] and [6] for simulation results).

The main difficulty of this approach is very high resource-consuming because of a very small wavelength. A large number of finite elements in $x_1$-direction is required to resolve the wave structure. The number of degrees of freedom lies in
the range of 10⁶-10⁷, which makes impossible, e.g. to compute the phase wave surface with appropriate accuracy.

## 3.3 Harmonic analysis (Dispersion relations)

The approach related to the harmonic analysis is developed in [4] and provides a method for the construction of travelling wave solutions feasible in the biosensor structure under the assumption of its unboundedness in the lateral and downward directions. This assumption is very realistic because real biosensor chips are imbedded up to the surface in very viscous damping media that suppresses the reflection of waves on the the side and bottom faces.

The algorithm is described here quite briefly (see [7] and [4] for more details). It is assumed that all the layers are infinite in $x_1$ and $x_2$ directions, the (top) fluid layer and the (bottom) substrate layer are semi-infinite in $x_3$ direction. The electrodes are not taken into account.

We are looking for solutions describing plain waves propagating in $x_1$ direction. This means that the displacements in the solid layers, the velocities in the fluid, and the electric potential in the substrate are of the form:

$$u_i(x_1, x_3) = a_i(x_3) \cos(\kappa x_1 - \omega t) + b_i(x_3) \sin(\kappa x_1 - \omega t), \quad (10)$$

$$v_i(x_1, x_3) = c_i(x_3) \cos(\kappa x_1 - \omega t) + d_i(x_3) \sin(\kappa x_1 - \omega t), \quad (11)$$

$$\varphi(x_1, x_3) = f(x_3) \cos(\kappa x_1 - \omega t) + g(x_3) \sin(\kappa x_1 - \omega t), \quad (12)$$

where $\kappa$ is the wave number, and $\omega$ is the circular frequency which is equal to the frequency of the voltage applied to the input electrodes in our case. Substitute (10) and (12) into (1) and (2) for the substrate; (10) into (7) and (9) for non-piezoelectric layers and for the aptamer layer; and (11) into (8) for the fluid layer. Collecting all coefficients of cos and sin yields a system of ordinary linear differential equations for the coefficients $a_i, b_i, c_i, d_i, f, g$ in each layer. Solving these systems for every layer, we obtain the representation of the functions $a_i, b_i, c_i, d_i, f, g$ in the following form (only the expression for the function $a = (a_1, a_2, a_3)$ is given here because the form of the other functions is similar):

$$a(x_3) = \sum_j D^j h^j e^{\lambda^j \kappa x_3}, \quad (13)$$

where $D^j$ are arbitrary coefficients, $\lambda^j$ and $h^j$ are eigenvalues and eigenvectors of the matrix of the corresponding system of differential equations. For the semi-infinite fluid and substrate layers, only terms decreasing towards $x_3$ for the fluid and towards $-x_3$ for the substrate, i.e. terms with negative $Re\lambda^j$ for the fluid and positive $Re\lambda^j$ for the substrate, are kept.

Every layer has its own set of coefficients $D^j$, eigenvalues $\lambda^j$, and eigenvectors $h^j$. To find any particular travelling wave solution in the whole structure we need to determine the coefficients $D^j$ for each layer, which is being done by substituting the expressions of the form (13) for the functions $a_i, b_i, c_i, d_i, f, g$
into (10)–(12) and then the resulting functions $u_i, v_i, \phi$ into the interface conditions outlined at the end of subsection 3.1 (see [6] for exact description of the interface conditions). Since all of the interface conditions are linear relations, the computation yields a homogeneous system of linear equations for the unknown coefficients $D^j$. Denote by $G(\omega, \kappa)$ the matrix of this system. Fix the circular frequency $\omega$ and denote the unknown phase velocity by $V = \omega/\kappa$ to consider $G$ as a function of $V$. The phase velocity is feasible if and only if the system has a nontrivial solution, which is equivalent to the condition $\det |\bar{G}^T(V)G(V)| = 0$, where $\bar{G}^T(V)$ the conjugate transpose of $G(V)$. The last equation can be easily solved because the computation of the left-hand-side runs very quickly even on an ordinary computer. Usually, there are several roots corresponding to different types of waves propagating with different phase velocities. Concerning the biosensor, the root corresponding to a shear wave, i.e. only $u_2 \neq 0$, is to be chosen.

Figure 4a shows the computed phase velocity contour for surface acoustic waves exited in the biosensor structure using the excitation frequency of 96 MHz. Figure 4b presents the slowness contour scaled by $10^3$. It is seen that the slowness contour is not convex (see remark at the end of Section 2).

<figure><img src="image_6.png" /><figcaption>Fig. 4. Characteristic contours for surface acoustic waves: (a) Phase velocity contour; (b) Slowness contour.</figcaption></figure>

# 4 Description of wave propagation

This section addresses the question how to describe the propagation of waves if the velocity surface (contour in the case of surface waves) is known. Let us first recall the classification of surfaces related to the wave propagation. Then, the applicability of Fermat's minimum time principle and the direct usage of eikonal equations will be discussed.
## 4.1 Characteristic surfaces

In acoustics, three characteristic surfaces are used to characterize the wave propagation (see [8] and [9]).

**Wave surface.** The wave surface (or the group velocity surface) describes the propagation of acoustic energy. This surface is involved in the formulation of minimum time principles, e.g. Fermat's law. The wave surface is the locus of points traced by the energy velocity vector $V_e$, drawn from a fixed point O, as the propagation direction varies. The propagation direction n (the normalized wave vector) is orthogonal to the wave surface (see Fig. 5a). It should be noticed that $V_e$ is as a rule not collinear to the wave vector in the case of anisotropic media.

**Phase velocity surface.** The phase velocity surface (see Figures 1a and 4a) describes the propagation of wave fronts. It defines the Hamiltonian of the eikonal equation. The phase velocity surface is obtained from the wave surface by projecting the vector $V_e$ onto the wave propagation direction n (see Fig. 5b) so that the phase velocity vector V is given by $V = (V_e \cdot n)n$. It should be noticed that the phase velocity surface can be constructed independently on $V_e$, e.g. as shown in Section 2 and Subsection 3.3, and the wave surface can then be defined through the phase velocity surface.

**Slowness surface.** The slowness surface (see Figures 1b and 4b) indicates the local convexity/concavity properties of the Hamiltonian of the eikonal equation. The slowness surface is related to the phase velocity surface by the inversion through the origin (see Fig. 5c), i.e. $m = n/|V|$. The energy velocity, $V_e$, is normal to the slowness surface at all points. Local concavities on the slowness surface can cause formation of cusps ("swallow tails") on the wave surface as it is shown in Fig. 5c: The arc (acb) is mapped into a "swallow tail" on the wave surface. This points out to the intersection of characteristics of the eikonal equation.

## 4.2 Fermat's principle

The Fermat principle describes how a ray, trajectory orthogonal to the wave front at all time instants, propagates from point A to point B. The principle says that the propagation time should be minimal. To express this, consider the minimization problem

$$T = \int_{A}^{B} dt = \int_{A}^{B} \frac{ds}{|V_e|} \rightarrow \min,$$

where the integrals are computed along rays. Let $x(\tau)$ be the parametrization of rays. Accounting for the relation $ds = |\dot{x}|d\tau$ yields

$$T = \int_{\tau_0}^{\tau_1} \frac{|\dot{x}|}{|V_e(x, \dot{x}/|\dot{x}|)|} d\tau := \int_{\tau_0}^{\tau_1} L(x, \dot{x}) d\tau \rightarrow \min.$$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Schematic explanation to characteristic surfaces: (a) Wave surface; (b) Phase velocity surface; (c) Slowness surface.</figcaption></figure>

Thus, feasible rays are solutions of the Euler equation

$$L_x - \frac{d}{d\tau} L_{\dot{x}} = 0.$$

This approach works well if the wave velocity $V_e(x, n)$ is well-defined for all directions $n$. It holds if the slowness surface is convex, which is as a rule violated in the case of anisotropic media. The next subsection discusses the method of direct solving eikonal equations.

## 4.3 Eikonal and Hamilton-Jacobi equations

Let $c(x, n)$ be the phase velocity depending on the spatial position $x$ and the propagation direction $n$, $|n| = 1$, Let $S(t, x)$ be the phase function that shows the phase of the wave at the time instant $t$ and at the point $x$. According to the results of Section 2, cf. equation (6), the eikonal equation reads

$$S_t - |\nabla S|c\left(x, \frac{\nabla S}{|\nabla S|}\right) = 0. \qquad (14)$$

If the Hamiltonian $c(x, p/|p|)|p|$ is convex in $p$, the method of characteristics can be used for solving equation (14) in the case of convex wave emitter. Under these conditions, the characteristics representing the rays do not intersect each other. If the convexity property is violated, equation (14) may not have classical solutions. Nevertheless, it is always uniquely solvable in the sense of viscous solutions (see [10] and [11]). Moreover, a unique viscous solution of (14) is the valid phase function. Thus, we arrive at the idea to use numerical methods of finding viscosity solutions of Hamilton-Jacobi equations. It should be noticed that common Lax-Friedrichs methods (see e.g. [12]) are not applicable in this case because they smooth solutions very strong. On the other hand, the authors
have developed numerical methods that do not contain any smoothing (see e.g. [13] and [14]). These methods assume that the Hamilton-Jacobi equations arise from conflict control problems (differential games). Therefore, their application requires solving the following problem: Given an eikonal equation, it is required to construct a differential game whose Hamiltonian coincides (up to the sign) with that of the eikonal equation. The next section shows how to do that using an unification technique proposed in [15].

## 5 Usage of differential games

Assume that a set $M \subset R^d$, $d = 2$ or 3, represents the shape of the acoustic wave emitter. For example, $M$ is a ball in the case of bulk crystal, and $M$ is the two dimensional area of the input electrodes in the case of biosensor. Let the game dynamics be described by the following system of ordinary differential equations:

$$ \dot{x} = f(x, u, v), \quad x \in R^d, t \in (-\infty, 0], \quad u \in P \subset R^a, v \in Q \subset R^b, \quad (15) $$

where $u$ and $v$ are control parameters of the first and second players, respectively. Introduce the signed distance, $\sigma$ to the set $M$ as follows: $\sigma(x) = \text{dist}(x, M)$ if $x \notin M$, and $\sigma(x) = -\text{dist}(x, R^d \setminus M)$ if $x \in M$. Consider the objective functional, $\gamma$, defined on the trajectories of (15) as follows:

$$ \gamma(x(\cdot)) = \min_{\tau \in [t, 0]} \sigma(x(\tau)). \quad (16) $$

The game is formalized using the concept of feedback strategies (see [16]). The value function is defined by the relation

$$ \Psi(t, x) = \max_{\mathcal{V}} \min_{x(\cdot) \in X(t, x, \mathcal{V})} \gamma(x(\cdot)), $$

where $\mathcal{V}$ is a feedback strategy of the second player, and the set $X(t, x, \mathcal{V})$ expresses the actions of the first player. This set consists of all limits of Euler trajectories of (15) which are obtained when the second player chooses $v \equiv \mathcal{V}(t_i, x(t_i))$ on each interval $[t_i, t_{i+1})$ of partitions of $[t, 0]$, and the first player uses admissible controls $u(\xi)$, $\xi \in [t, 0]$. In doing that, all possible partitions whose diameter tends to zero and all admissible controls of the first player are exhausted. All Euler trajectories start at $t$ from the initial state $x$.

The value function is locally bounded and Lipschitzian (see e.g. [11]).

Define the Hamiltonian

$$ H(x, p) = \max_{v \in Q} \min_{u \in P} \langle p, f(x, u, v) \rangle, \quad p \in R^d, \quad (17) $$

and consider the Hamilton-Jacobi-Bellman-Isaacs equation

$$ \Psi_t + H(x, \Psi_x) = 0, \quad \Psi(0, x) = \sigma(x). \quad (18) $$
It is proven in [13] that the value function of the game (15) with the objective functional (16) is a viscosity solution of (18). Therefore, the following proposition holds:

**Proposition 1.** Let $c(x, p/|p|) |p|$ be the Hamiltonian of the eikonal equation. If

$$H(x, p) = -c(x, p/|p|) |p|, \quad p \in \mathbb{R}^d,$$

then the wave front at any time instant $t \ge 0$ is given by the relation

$$\{x : \Psi(-t, x) = 0\}.$$

This proposition opens the way to use numerical methods of finding viscosity solutions of Hamilton-Jacobi equations. The only question consists in constructing an appropriate differential game whose Hamiltonian satisfies the condition of Proposition 1. The next subsection discusses this task.

## 5.1 Unification

Denote $E(x, p) = -c(x, p/|p|) |p|$, $p \in \mathbb{R}^d$. To find a differential game whose Hamiltonian coincides with $E$, the technique of unification (see [15]) can be used. Consider the following conflict control system

$$\dot{x} = E(x, p)p + q, \quad x, p, q \in \mathbb{R}^d, \quad |p| = 1, \quad |q| = \lambda, \quad \langle p, q \rangle \ge 0. \quad (19)$$

Here, $q$ is the control parameter of the first player who strives to minimize the objective functional (16), whereas $p$ is the control parameter of the second player who maximizes the objective functional. The parameter $\lambda$ is a constant which is greater than the Lipschitz constant of the function $E$ in $p$.

**Proposition 2.** If $|E(x, p_1) - E(x, p_2)| < \lambda|p_1 - p_2|$, $p_1, p_2 \in \mathbb{R}^d$, $|p_1| = 1$, $|p_2| = 1$, then the Hamiltonian of the game (19) satisfies the relation

$$H_{(19)}(x, s) := \max_{\substack{|p|=1 \\ \langle p, q \rangle \ge 0}} \min_{\substack{|q|=\lambda, \\ (p, q) \ge 0}} \langle E(x, p)p + q, s \rangle = E(x, s),$$

and, therefore, (19) is the required differential game.

It should be noticed that the proof of this proposition essentially uses the positive homogeneity of the function $p \to E(\cdot, p)$. Therefore, the unification procedure cannot be applied in the case of absence of positive homogeneity.

Assume now that $|p| = 1$, then $E(x, p)p = -c(x, p)p$. Taking into account that the points $c(x, p)p$ exhaust the phase velocity surface $V_{\text{surf}}(x)$ when varying $p$, the game (19) can be rewritten as

$$\dot{x} = -p + q, \quad p \in V_{\text{surf}}(x), \quad |q| = \lambda, \quad \langle p, q \rangle \ge 0.$$

Assuming that only the velocity magnitude depends on the spatial position yields the game

$$\dot{x} = -a(x)p + q, \quad p \in V_{\text{surf}}, \quad |q| = \lambda, \quad \langle p, q \rangle \ge 0, \quad (20)$$
where the phase velocity surface $V_{surf}$ does not depend on $x$. Notice that the coefficient $a(x)$ is necessary to take into account the damping area of the biosensor, the phase velocity is strongly reduced there. Therefore, $a(x) \equiv 1$ outside the damping area, and $a(x) \to 0$ towards the outer boundary.

Numerical methods developed by the authors (see [13] and [14]) provide an effective tool for finding the value function of the game (20) with the objective functional (16). The next section uses these methods and demonstrates a good reconstruction of wave patterns.

# 6 Simulation of wave propagation using the unified differential game

Now, the differential game (20) with the objective functional (16) is used for the computation of wave fronts. First, consider waves in $LiTaO_3$ piezoelectric crystals. Assume that the wave emitter, the set $M$, is a ball of radius 0.1, and $a(x) \equiv 1$, i.e. there is no damping area. The phase velocity surface is shown in Figure 1a. Thus, all data required for the formulation of the differential game (16) and (20) are available. Application of a finite difference upwind scheme described in [13] yields an approximation, $\tilde{\Psi}$, of the value function of the game. According to Proposition 1, the wave front at any time instant $t \ge 0$ is given by the relation $\{x \in \mathbb{R}^3 : \tilde{\Psi}(-t, x) = 0\}$. Figure 6 shows the wave front at $t = 0.25$, 1, and 2 ms.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Wave fronts in lithium tantalate $LiTaO_3$ piezoelectric crystal (quasi shear wave): (a) $t = 0.25$ ms; (b) $t = 1$ ms; (c) $t = 2$ ms.</figcaption></figure>

Consider now wave fronts for surface shear waves in the biosensor structure (see Section 3). The phase velocity contour is shown in Figure 4a. Figure 7 shows the position of the wave front for a time sequence with a small sample time. The wave emitter and the damping area are easily recognizable in this figure.

Figures 8 and 9 show the case of two wave emitters. Figure 10 shows the wave propagation in the presence of a hole. The hole is interpreted as an obstacle such
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Surface wave fronts in the biosensor structure. The wave emitter is an unclosed ring.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Surface wave fronts in the biosensor structure with two emitters: a half-ring and a ring.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Surface wave fronts in the biosensor structure with two half-rings as the emitters.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Surface wave fronts in the biosensor structure in the presence of a round hole. The wave emitter is a half-ring.</figcaption></figure>

that the trajectories of the game (20) can not penetrate therein. This case is numerically processed using a method for finding value functions in games with state constraints (see [14]).

## Conclusion

The technique presented in this paper is also appropriate for the numerical treatment of arbitrary Hamilton-Jacobi equations with positive homogeneous Hamiltonians. As it was seen, the unification procedure described in Section 5.1 does not use any specific features of the Hamiltonian with the exception of the positive homogeneity which is necessary for the proof of Proposition 2. This opens new possibilities of investigation of physical processes related to optimality principles involving non-convex Lagrangians and Hamiltonians.

## References

1. Melikyan, A.A.: Generalized Solutions of First Order PDEs. Birkhäuser, Boston (1998)

2. Barles, G.: Remarks on a Flame Propagation Model. Technical Report 464, INRIA (1985)

3. Hoffmann, K.-H., Botkin, N.D., Starovoitov, V.N.: Homogenization of Interfaces Between Rapidly Oscillating Fine Elastic Structures and Fluids. SIAM J. Appl. Math. 65(3), 983–1005 (2005)

4. Botkin, N.D., Hoffmann, K.-H., Pykhteev, O.A., Turova, V.L.: Dispersion Relations for Acoustic Waves in Heterogeneous Multi-Layered Structures Contacting with Fluids. Journal of the Franklin Institute 34(5), 520–534 (2007)

5. Botkin, N.D., Turova, V.L.: Mathematical Models of a Biosensor. Applied Mathematical Modelling 28(6), 573–589 (2004)
6. Botkin, N.D., Hoffmann, K.-H., Pykhteev, O.A., Starovoitova, B.N., Turova, V.L.: Two Complementary Approaches in Modelling a Biosensor. In: Hamza. M.N. (ed.) 15th IASTED International Conference on Applied Simulation and Modelling, pp. 525–530 ACTA Press, Anaheim-Calgary-Zurich (2006)

7. Botkin, N.D., Hoffmann, K.-H., Pykhteev, O.A., Turova, V.L.: Numerical Computation of Dispersion Relations for Multi-Layered Anisotropic Structures. In: Technical Proceedings of 2004 Nanotechnology Conference and Trade Show, Vol. 2, pp. 411–414. NSTI, Boston (2004)

8. Auld, B.A.: Acoustic Fields and Waves in Solids. I. Krieger Publishing Company, Malabar (1972)

9. Royer, D., Dieulesaint, E.: Elastic Waves in Solids I: Free and Guided Propagation. Springer-Verlag, Berlin Heidelberg (2000)

10. Crandall, M.G., Lions, P.L.: Viscosity Solutions of Hamilton-Jacobi Equations. Trans. Amer. Math. Soc. 277, 1–47 (1983)

11. Subbotin, A. I.: Generalized Solutions of First Order PDEs: The Dynamical Optimization Perspective. Birkhäuser, Boston (1995)

12. Mitchell, I.: Application of Level Set Methods to Control and Reachability Problems in Continuous and Hybrid Systems. PhD Thesis, Stanford University (2002)

13. Botkin, N.D., Hoffmann, K.-H., Turova, V.L.: Stable Numerical Schemes for Solving Hamilton-Jacobi-Bellman-Isaacs Equations. SIAM Journal on Scientific Computing 33(2), 992–1007 (2011)

14. Botkin, N.D., Hoffmann, K.-H., Mayer, N., Turova, V.L.: Approximation Schemes for Solving Disturbed Control Problems with Non-Terminal Time and State Constraints. Analysis 31, 355–379 (2011)

15. Krasovskii, N.N.: On the Problem of the Unification of Differential Games. Dokl. Akad. Nauk SSSR 226(6), 1260 – 1263 (1976)

16. Krasovskii, N. N., Subbotin, A. I.: Game-Theoretical Control Problems. Springer, New York (1988)