\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small Eighth Mississippi State - UAB
Conference on Differential Equations and Computational
Simulations. {\em Electronic Journal of Differential Equations},
Conf. 19 (2010), pp. 235--244.\newline ISSN: 1072-6691. URL:
http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{235}
\title[\hfilneg EJDE-2010/Conf/19/\hfil Numerical solution]
{Numerical solution to nonlinear Tricomi equation using WENO
schemes}
\author[A. Sescu, A. A. Afjeh, C. Sescu \hfil EJDE/Conf/19 \hfilneg]
{Adrian Sescu, Abdollah A. Afjeh, Carmen Sescu} % not in alphabetical order
\address{Adrian Sescu \newline
MIME Department, University of Toledo,
2801 West Bancroft Street,
Toledo, OH 43606-3390OH, USA}
\email{adrian.sescu@utoledo.edu Tel +1 419 530 8160, fax: +1 419
530 8206}
\address{Abdollah A. Afjeh \newline
MIME Department, University of Toledo,
2801 West Bancroft Street,
Toledo, OH 43606-3390OH, USA}
\email{aafjeh@utoledo.edu}
\address{Carmen Sescu \newline
MIME Department, University of Toledo,
2801 West Bancroft Street,
Toledo, OH 43606-3390OH, USA}
\email{carmen.sescu@utoledo.edu}
\thanks{Published September 25, 2010.}
\subjclass[2000]{76Q05, 35L60, 35L65, 65M22}
\keywords{Nonlinear aeroacoustics; hyperbolic conservation law;
\hfill\break\indent discretized equations}
\begin{abstract}
Nonlinear Tricomi equation is a hybrid (hyperbolic-elliptic)
second order partial differential equation, modelling the sonic
boom focusing. In this paper, the Tricomi equation is transformed
into a hyperbolic system of first order equations, in conservation
law form. On the upper boundary, a new mixed boundary condition
for the acoustic pressure is used to avoid the inclusion of the
Dirac function in the numerical solution. Weighted Essentially
Non-Oscillatory (WENO) schemes are used for the spatial
discretization, and the time marching is carried out using the
second order accurate Runge-Kutta total-variation diminishing
(TVD) scheme.
\end{abstract}
\maketitle \numberwithin{equation}{section}
\section{Nonlinear Tricomi Equation}
The weak shock wave focusing at a smooth caustic has been elucidated as
a result of the concern about the focusing of sonic boom induced
by maneuvering supersonic aircraft. Guiraud \cite{guiraud}
elaborated a consistent theory including both diffraction and
nonlinear effects at first order and leading to the nonlinear
Tricomi equation. In a review report, Coulouvrat \cite{coulouvrat}
extended the previous results to a three-dimensional heterogeneous
medium. Hunter and Keller \cite{hunter} showed that the nonlinear
Tricomi equation occurs for the general case of weakly nonlinear
wave solutions of a system of hyperbolic equations. Kandil and
Zheng \cite{kandil2} solved the nonlinear, non-conservative
Tricomi equation using a frequency-domain scheme, a time domain
scheme and a time domain with overlapping grid scheme. A
conservative form of the nonlinear Tricomi equation with the
pressure potential as the dependent variable has been developed
and solved using a time-domain scheme. The four schemes have been
applied to several incoming waves which include an $N$-wave, a
Concorde aircraft wave and symmetric and asymmetric flat-top and
ramp-top waves. Kandil and Khasdeo~\cite{kandil3} investigated the
effects of several parameteres (longitudinal size fo the
computational domain, the incoming shock footprint width and its
level of strength) on the computational solution to the nonlinear
Tricomi equation.
According to the catastrophe theory \cite{berry}, the pressure
near the fold caustic can be shown to be a function of two
independent variables only: the distance to the caustic $z$, and
the phase of the signal $t$. The corresponding dimensionless
variables are:
(1) the dimensionless delayed time:
$\bar{\tau}=[t-x(1-z/R_{\rm sec}/c_0)]/T_{\rm ac}$, where $x$ is the
coordinate tangent to the caustic, $R_{\rm sec}$ is the radius of
curvature of the intersection of the caustics with the $Oxz$
plane, $c_0$ is the ambient sound speed, and $T_{\rm ac}$ is the
characteristic duration of the acoustical signal;
(2) the dimensionless distance to the caustic: $\bar{z}=z/\delta$,
where $\delta$ is the thickness of the diffraction boundary; and
(3) the dimensionless pressure: $\bar{p}=(p-p_0)/p_{\rm ac}$, where
$p_0$ is the ambient pressure, and $p_{\rm ac}$ is the maximal
pressure.
Using these dimensionless variables, the Euler equations are
simplified to get the nonlinear Tricomi equation \cite{guiraud}:
\begin{equation} \label{1}
\frac{\partial^2{\bar{p}}}{\partial{\bar{z}^2}} -
\bar{z}\frac{\partial^2{\bar{p}}}{\partial{\bar{\tau}^2}} +
\frac{\mu}{2}
\frac{\partial^2{\bar{p^2}}}{\partial{\bar{\tau}^2}} = 0,
\end{equation}
where the coefficient $\mu=2\beta M_{\rm ac}[R_{\rm cau}/(2c_0 T_{\rm ac})
]^{2/3}$ measures the nonlinear effects relative to the
diffraction effects (here $R_{\rm cau}=1/(1/R_{\rm sec} + 1/R_{\rm ray})$,
and $R_{\rm ray}$ is the radius of curvature of the projection of
the acoustical ray on the $Oxz$ plane), $\beta$ is the nonlinearity
parameter of the medium, and $M_{\rm ac}$ is the acoustical Mach number.
In figure \ref{f1}, a schematic of the focusing of a $N$-wave i
s shown. The domain consists of two zones: the hyperbolic
zone which is located above the $\bar{z}=0$ axis, and the elliptic
zone which is located below the $\bar{z}=0$ axis. The amplification
of the $N$-wave occurs on the caustic, in the vicinity of the
$\bar{\tau}=0$, $\bar{z}=0$.
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig1}
\end{center}
\caption{ $N$-wave focusing}
\label{f1}
\end{figure}
The associated boundary conditions are:
(1) In time, for a transient signal, the medium is not
perturbed before or after the acoustical wave has passed:
$\bar{p}(\bar{z},\bar{\tau} \to \pm \infty) = 0$
or, for a periodic signal (with a period T), one simply gets
$\bar{p}(\bar{z},\bar{\tau}+ T) = \bar{p}(\bar{z},\bar{\tau})$;
(2) Away from the caustic in the shadow zone, the
acoustic pressure decays exponentially:
$\bar{p}(\bar{z} \to -\infty,\bar{\tau}) = 0$;
(3) Away from the caustic on the illuminated side, the
field matches the geometrical acoustics approximation:
\begin{equation} \label{2}
\bar{p}(\bar{z} \to +\infty,\bar{\tau}) =
\bar{z}^{-1/4}
[F(\bar{\tau}+\frac{2}{3}\bar{z}^{3/2}) +
G(\bar{\tau}-\frac{2}{3}\bar{z}^{3/2})]
\end{equation}
The function $F$ represents the incoming wave, and is supposed
to be known. On the contrary, the function $G$ is the time
waveform along the outgoing ray. Unlike the incoming signal $F$,
the outgoing signal $G$ has undergone the diffraction effects
after having tangented the caustic, and is unknown. To eliminate
this unknown function, the matching boundary condition \eqref{2}
can be written as a ``radiation condition", by a combination of
its derivatives with respect to $\bar{z}$ and $\bar{\tau}$:
\begin{equation} \label{3}
\bar{z}^{ 1/4}\frac{\partial{\bar{p}}}{\partial{\bar{\tau}}} +
\bar{z}^{-1/4}\frac{\partial{\bar{p}}}{\partial{\bar{z}}}
+ z^{-5/4}\frac{\bar{p}}{4} \to
2\frac{dF}{d\tau}
\big(\bar{\tau}+\frac{2}{3}\bar{z}^{3/2}\big),\quad
\text{as } \bar{z} \to +\infty
\end{equation}
\subsection{Conservation Law Form}
In what follows the bar notation for both the dependent and
independent variables will be dropped. Upon introducing the
pseudo-time derivatives and the new dependent variable $q$, the
conservation law form of the nonlinear Tricomi equation can be
written as
\begin{equation} \label{4}
\frac{\partial{\mathbf{u}}}{\partial{t}} +
\operatorname{div}\mathbf{f}(\mathbf{u}) = 0, \quad
\text{in }\Omega \times (0,\theta)
\end{equation}
where $\Omega$ is a real domain in $R^2$, $(0,\theta)$
is a time interval, and
\begin{equation} \label{5}
\mathbf{u} = \begin{pmatrix}
p \\
q \end{pmatrix}; \mathbf{f}(\mathbf{u}) =
\mathbf{f}_1(\mathbf{u}) \vec{i} + \mathbf{f}_2(\mathbf{u})
\vec{j} = \begin{pmatrix}
zp - \frac{\mu}{2}p^2 \\
q \end{pmatrix} \vec{i} +
\begin{pmatrix}
q \\
p \end{pmatrix} \vec{j}
\end{equation}
\noindent\textbf{Definition.}
The conservation law \eqref{4} is
hyperbolic if any real combination of the Jacobians
$\sum_{i=1}^{d}\psi_i \partial{\mathbf{f}_i}/\partial{\mathbf{u}}$
has 2 real eigenvalues and a complete set of eigenvectors.
The real combination of the Jacobians is given by
\begin{equation} \label{6}
J(\mathbf{u}) =
\begin{bmatrix}
\psi_1(z-\mu p) & \psi_2 \\
\psi_2 & \psi_1 \end{bmatrix}
\end{equation}
and the eigenvalues are
$
\lambda_{1,2} = 1/2
[\psi_1(z-\mu p+1) \pm \sqrt{\psi_1^2(z-\mu p-1)^2+4\psi_2^2}]
$,
which are real for any combination of $\psi_1$ and $\psi_2$.
The corresponding linearly independent right eigenvectors are:
\begin{equation} \label{7}
\mathbf{r}_1 = C_1
\begin{pmatrix}
A + \sqrt{A^2 + 4\psi_2^2} \\
2\psi_2 \end{pmatrix}; \quad
\mathbf{r}_2 = C_2
\begin{pmatrix}
2\psi_2 \\
-A - \sqrt{A^2 + 4\psi_2^2} \end{pmatrix};
\end{equation}
where $A = \psi_1(z-\mu p-1)$.
As a result, the conservation law \eqref{4} is hyperbolic.
The pseudo-time derivative in equation \eqref{4} has been
introduced to change the character of the equation from
mixed (elliptic/hyperbolic) to hyperbolic. The iterations
in time have the scope to drop the time derivatives to zero
(the residuals go to zero), so that the numerical solution
satisfies the original equation \eqref{1}.
\subsection{Modified Boundary Conditions}
For $z \to +\infty$ the boundary condition is a little
problematic: equation \eqref{3} contains the derivative of $F$
with respect to $\tau$ and, in the case of an incoming shock wave,
this leads to a boundary condition with a sharp singularity (delta
Dirac function), which may be problematic from the numerical point
of view.
In this work, a different treatment of the boundary condition for
$z \to +\infty$ is adopted. Function $F(x)$ which represents the
incoming wave is assumed to be zero except in a finite interval in
the vicinity of $x=0$. This is the case of an N wave, for example,
modelling the sonic boom focusing. When equation \eqref{4} is
solved numerically, the domain is truncated and the boundary
conditions at infinity become conditions for large $\tau$ or $z$.
If the upper boundary ($z \to +\infty$) is considered far enough
from 0, then the incoming and the outgoing waves are separated by
the $z$-axis. This suggests the division of the upper boundary in
two regions: $\tau < 0$ and $\tau \ge 0$. Thus the boundary
condition on the upper side can be written as
\begin{equation} \label{8}
\Phi(-\tau)p + \Phi(\tau)
\Big(z^{ 1/4}\frac{\partial{p}}{\partial{\tau}} +
z^{-1/4}\frac{\partial{p}}{\partial{z}}
+ z^{-5/4}\frac{\bar{p}}{4}\Big)
\to
\Phi(-\tau)F\Big(\tau+\frac{2}{3}z^{3/2}\Big)
\end{equation}
as $z \to +\infty$, where $\Phi(\tau)$ is the Heaviside
step function.
\subsection{Linear Case: Analytical Solution}
An analytical solution~\cite{auger} to the linear Tricomi equation
is presented in this section, which will be used to validate
the numerical solution. The linear Tricomi equation is obtained
by setting $\mu=0$:
\begin{equation} \label{9}
z\frac{\partial^2{p}}{\partial{\tau^2}}
- \frac{\partial^2{p}}{\partial{z^2}} = 0
\end{equation}
with boundary conditions defined as before.
An analytical solution can be found using the Fourier analysis
\cite{auger} in the form
\begin{equation} \label{13}
p(\tau,z) =
TF^{-1} [
\sqrt{2}(1+i \operatorname{sgn}(\omega))|\omega|^{\frac{1}{6}}
Ai(-|\omega|^{2/3}z) TF(F)]
\end{equation}
where $TF$ and $TF^{-1}$ stand for the Fourier transform and its
inverse:
\begin{equation} \label{10}
\hat{p}(z,\omega) =
\int\limits_{-\infty}^{+\infty}p(\tau,z)e^{-i\omega \tau} d\tau; \quad
p(\tau,z) = \frac{1}{2\pi}
\int\limits_{-\infty}^{+\infty}\hat{p}(z,\omega)e^{i\omega \tau}
d\omega
\end{equation}
The analytical solution at $z=0$, for an incoming $N$-wave in
the form
\begin{equation}\label{14}
F(\tau) = \begin{cases}
-\tau/T, &\text{if } |\tau|0$, the geometrical acoustics
approximation \eqref{2} can be used, where function $G$
can be explicitly determined, \cite{auger2}
\begin{equation} \label{17}
G(\tau) =
-\frac{T}{\pi}
[2 + \frac{\tau}{T}
\ln \frac{|T-\tau|}{|T+\tau|}]
\end{equation}
\section{WENO Schemes}
WENO schemes have been designed in recent years as a class of
high order finite volume or finite difference schemes to solve
hyperbolic conservation laws, with the property of maintaining
both uniform high order accuracy and an essentially non-oscillatory
solution in the vicinity of discontinuities or large gradients.
The first WENO scheme is constructed in~\cite{liu} for a
third-order finite volume version in one space dimension.
In~\cite{jiang}, third and fifth-order finite difference WENO
schemes in multi space dimensions are constructed, with a general
framework for the design of the smoothness indicators and
nonlinear weights. WENO schemes are designed based on the
successful ENO schemes in~\cite{harten,shu1}. Both ENO and WENO
schemes use the idea of adaptive stencils in the reconstruction
procedure based on the local smoothness of the numerical solution
to automatically achieve high order accuracy and a non-oscillatory
property near discontinuities. ENO uses just one (optimal in some
sense) out of many candidate stencils when doing the
reconstruction, while WENO uses a convex combination of all the
candidate stencils, each being assigned a nonlinear weight which
depends on the local smoothness of the numerical solution based on
that stencil. For a detailed review of ENO and WENO schemes, we
refer to the lecture notes~\cite{shu2}.
\subsection{Upwind WENO Schemes}
Upwind WENO schemes of third and fifth order of accuracy are
applied in this work, and therefore they are briefly introduced
here. The schemes are written as
\begin{equation} \label{18}
\hat{f}_{i+1/2}^{\pm} = \sum_{r=0}^{k-1}\omega_{r}^{\pm}\hat{f}_{i+1/2,r}^{\pm}
\end{equation}
where $\omega_{r}^{\pm}$ and $\alpha_{r}^{\pm}$ are given by
\begin{equation} \label{23}
\omega_{r} = \frac{\alpha_r}{\sum_{l=0}^{k-1}\alpha_{l}},
\quad\text{where }
\alpha_{r} = \frac{C_r}{(\epsilon + IS_r)^2},
\end{equation}
where $C_r$ are the ideal weights and the parameter $\epsilon$
is a small number ($10^{-10}$ to $10^{-6}$) introduced to avoid
the cancellation of the denominator, and $IS_r$ is a measure of
the smoothness.
\subsubsection{Third Order ($r=2$)}
The third order (WENO3) scheme uses a three-point stencil
$[x_{i-1},x_{i+1}]$ or $[x_{i},x_{i+2}]$, corresponding to
$\hat{f}_{i+1/2}^{+}$ or $\hat{f}_{i+1/2}^{-}$, respectively.
There are two ENO candidate stencils (for $\hat{f}_{i+1/2}^{+}$
or $\hat{f}_{i+1/2}^{-}$), given by
\begin{equation} \label{19}
\begin{gathered}
\hat{f}_{i+1/2,0}^{+} =
- \frac{1}{2} \hat{f}_{i-1}^{+}
+ \frac{3}{2} \hat{f}_{i}^{+} \\
\hat{f}_{i+1/2,1}^{+} =
+ \frac{1}{2} \hat{f}_{i}^{+}
+ \frac{1}{2} \hat{f}_{i+1}^{+}
\end{gathered}
\end{equation}
with the corresponding smoothness indicators,
\begin{equation} \label{20}
IS_0^{+} = \Big(\hat{f}_{i+1}^{+} - \hat{f}_{i}^{+}\Big)^2,\quad
IS_1^{+} = \Big(\hat{f}_{i}^{+}- \hat{f}_{i-1}^{+}\Big)^2
\end{equation}
The optimal weights are $C_0 = 1/3$ and $C_1 = 2/3$.
\subsubsection{Fifth Order ($r=3$)}
The fifth order (WENO5) scheme uses a five-point stencil
$[x_{i-2},x_{i+2}]$ or $[x_{i-1},x_{i+3}]$, corresponding
to $\hat{f}_{i+1/2}^{+}$ or $\hat{f}_{i+1/2}^{-}$, respectively.
There are three ENO candidate stencils (for $\hat{f}_{i+1/2}^{+}$
or $\hat{f}_{i+1/2}^{-}$), given by
\begin{equation} \label{21}
\begin{gathered}
\hat{f}_{i+1/2,0}^{+} =
\frac{2}{6} \hat{f}_{i-2}^{+}
- \frac{7}{6} \hat{f}_{i-1}^{+}
+ \frac{11}{6} \hat{f}_{i}^{+} \\
\hat{f}_{i+1/2,1}^{+} =
- \frac{1}{6} \hat{f}_{i-1}^{+}
+ \frac{5}{6} \hat{f}_{i}^{+}
+ \frac{2}{6} \hat{f}_{i+1}^{+} \\
\hat{f}_{i+1/2,2}^{+} =
\frac{2}{6} \hat{f}_{i}^{+}
+ \frac{5}{6} \hat{f}_{i+1}^{+}
- \frac{1}{6} \hat{f}_{i+2}^{+}
\end{gathered}
\end{equation}
with the corresponding smoothness indicators
\begin{equation} \label{22}
\begin{gathered}
IS_0^{+} = \frac{13}{12}\Big( \hat{f}_{i-2}^{+}
- 2 \hat{f}_{i-1}^{+}
+ \hat{f}_{i}^{+}\Big)^2
+ \frac{1}{4}\Big(
\hat{f}_{i-2}^{+}
- 4 \hat{f}_{i-1}^{+}
+ 3 \hat{f}_{i}^{+}
\Big)^2 \\
IS_1^{+} =
\frac{13}{12}\Big( \hat{f}_{i-1}^{+}
- 2 \hat{f}_{i}^{+}
+ \hat{f}_{i+1}^{+}\Big)^2
+ \frac{1}{4}\Big( \hat{f}_{i-1}^{+}
- \hat{f}_{i+1}^{+}\Big)^2 \\
IS_2^{+} =
\frac{13}{12}\Big( \hat{f}_{i}^{+}
- 2 \hat{f}_{i+1}^{+}
+ \hat{f}_{i+2}^{+}\Big)^2
+ \frac{1}{4}\Big( 3 \hat{f}_{i}^{+}
- 4 \hat{f}_{i+1}^{+}
+ \hat{f}_{i+2}^{+}\Big)^2
\end{gathered}
\end{equation}
The optimal weights are $C_0 = 1/10$, $C_1 = 6/10$ and
$C_2 = 3/10$.
\subsection{WENO Scheme for Nonlinear Tricomi Equation}
Consider the nonlinear Tricomi Equation written as a system
of conservation law (equation~\eqref{4}). The procedure of
calculating the numerical fluxes at a cell face using WENO
schemes proposed by Jiang and Shu~\cite{jiang} can be
summarized as follows:
\begin{itemize}
\item[1.] Projection to the characteristic space.
\item[2.] Flux splitting.
\item[3.] WENO interpolation.
\item[4.] Reverse projection.
\end{itemize}
The left eigenvector matrix $\mathbf{R}^{-1}_{i+1/2}$ is used
for the projection of the conservative variables and inviscid fluxes
to the characteristic space,
\begin{equation} \label{23b}
\mathbf{U} = \mathbf{R}^{-1}_{i+1/2} \mathbf{u}, \quad
\mathbf{F}_{i+1/2} = \mathbf{R}^{-1}_{i+1/2} \mathbf{f}_{i+1/2},
\end{equation}
The projected flux is split by flux splitting methods such as
the Lax-Friedrichs method:
\begin{equation} \label{24}
F^{(l)\pm} = \frac{1}{2}(F^{(l)} \pm \lambda_{\rm max}^{(l)} U{(l)})
\end{equation}
where the superscript $l$ denotes each characteristic field,
the signs $+$ and $-$ mean the positive and negative flux parts,
the scalar variables $F$ and $U$ are the elements of the vectors
$\mathbf{F}$ and $\mathbf{U}$, respectively, and the eigenvalue
$\lambda_{\rm max}^{(l)}$ is defined over a grid line:
$\lambda_{\rm max}^{(l)} = \max |\lambda_{i}^{(l)}|$,$
i=1,\dots ,N_x$. Then, the numerical fluxes obtained by the WENO
interpolation are given by
\begin{equation} \label{25}
\hat{\mathbf{F}}_{i+1/2} =
\hat{\mathbf{F}}_{i+1/2}^{+} + \hat{\mathbf{F}}_{i+1/2}^{-}
\end{equation}
where $\hat{\mathbf{F}}_{i+1/2}^{\pm}$ are the interpolated fluxes
at the cell face $i + 1/2$. Finally, the fluxes are obtained by
the reverse projection as
\begin{equation} \label{26}
\hat{\mathbf{f}}_{i+1/2} =
\mathbf{R}_{i+1/2}\hat{\mathbf{F}}_{i+1/2}
\end{equation}
After the spatial reconstruction is obtained with WENO schemes,
a set of ordinary differential equations is obtained:
\begin{equation} \label{27}
\frac{du}{dt} = L(u)
\end{equation}
where $L$ is the spatial discrete operator. This set of ordinary
differential equations is discretized here using the second
order TVD Runge-Kutta method~\cite{liu}:
\begin{gather*} %\label{21}
u^{(0)} = u^n \\
u^{(1)} = u^{(0)}+\Delta t L(u^{(0)}) \\
u^{n+1} = \frac{1}{2}u^{(0)}+\frac{1}{2}u^{(1)}+\frac{1}{2}\Delta t L(u^{(1)})
\end{gather*}
\section{$N$-wave Focusing}
An $N$-wave is considered as the incoming signal (function $F$),
with an amplitude of $z_{up}^{-1/4}$, where $z_{up}$ is the upper
limit of the domain, along $z$ direction. In order to be able to
apply the boundary conditions for large stencils
(such as WENO schemes), a number of additional ghost cells are
considered outside the domain. The values of the dependent
variables in the ghost points are set to very large numbers
(the order of $10^{4}$).
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig2a}
\includegraphics[width=0.47\textwidth]{fig2b}
\end{center}
\caption{Acoustic pressure versus $\tau$ for $z=2$ (linear solution):
WENO3 (left); WENO5 (right)}
\label{f2}
\end{figure}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig3a}
\includegraphics[width=0.47\textwidth]{fig3b}
\end{center}
\caption{Acoustic pressure versus $\tau$ for $z=1.2$ (linear solution):
WENO3 (left); WENO5 (right)}
\label{f3}
\end{figure}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig4a}
\includegraphics[width=0.47\textwidth]{fig4b}
\end{center}
\caption{Acoustic pressure versus $\tau$ for $z=0$ (linear solution):
WENO3 (left); WENO5 (right)}
\label{f4}
\end{figure}
Figures~\ref{f2} and~\ref{f3} show comparisons between the numerical
solution and the analytical solution for $z=1.2$ and $z=2$,
respectively; the analytical solution in the radiation zone is
obtained using equation \eqref{17}, for both cases ($z=1.2$ and $z=2$).
Also, the comparison between the numerical and the analytical
solutions in linear case, obtained using equation \eqref{15},
for $z=0$, is shown in figure~\ref{f4}. The nonlinear solver
was used to determine the numerical linear solution, by setting
a very small value to the coefficient
$\mu=2\beta M_{\rm ac} [R_{\rm cau}/(2c_0 T_{\rm ac})]^{2/3}$ defined in
section 1. WENO5 schemes are more accurate than WENO3 schemes,
as expected: the $L^2$-error is in the order of $10^{-2}$ for
WENO3 schemes, and in the order of $10^{-4}$ for WENO5 schemes.
Figure~\ref{f5} shows contours of acoustic pressure in nonlinear
case, for third (left) and fifth (right) order WENO schemes.
The $N$-wave is entering the domain on the left side, and is
amplified on the caustic (in the vicinity of $(z=0, \tau=0)$);
the amplitude becomes three to four times larger. Figure~\ref{f6}
shows the acoustic pressure as a function of $\tau$ for $z=0.207$,
corresponding to the maximum pressure.
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig5a}
\includegraphics[width=0.47\textwidth]{fig5b}
\end{center}
\caption{ Acoustic pressure contours (nonlinear solution): \hfill
WENO3 (left); WENO5 (right)}
\label{f5}
\end{figure}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig6a}
\includegraphics[width=0.47\textwidth]{fig6b}
\end{center}
\caption{ Acoustic pressure versus $\tau$ for $z=0.207$
(nonlinear solution): WENO3 (left); WENO5 (right)}
\label{f6}
\end{figure}
\subsection*{Concluding Remarks}
Nonlinear Tricomi equation, modelling the sonic boom focusing, has
been transformed in hyperbolic conservation law form, by
introducing a new dependent variable and using pseudo-time
derivatives. The boundary condition on the illuminated side has
been modified; it was split into two parts: on the left side,
Dirichlet boundary condition has been imposed for pressure; on the
right side, a radiation condition has been imposed. The equation
was solved using WENO schemes of third and fifth order of
accuracy, and the time marching was carried out using a second
order TVD Runge-Kutta method. The numerical solution for the
linear case was compared to the analytical solution found via
Fourier analysis, and the agreement was good.
\begin{thebibliography}{00}
\bibitem{abramowitz} Abramowitz, M. and Stegun, I.A.:
Handbook of Mathematical Functions.
National Bureau of Standards, Washington (1965).
\bibitem{auger2} Auger, T.:
Mod\'elisation et simulation num\'erique de la focalisation
d'ondes de chocs acoustiques en milieu en mouvement. Application
\`a la focalisation du bang sonique en acc\'el\`eration. PhD
thesis, Universit\'e Paris VI, (2001).
\bibitem{auger} Auger, T. and Coulouvrat, F.:
Focusing of shock waves at smooth caustics, 15th International
Symposium on Nonlinear Acoustics, Gottingen, Germany,
September 1-4, (1999).
\bibitem{berry} Berry, M. V.:
{Waves and Thoma's theorem,} Adv. Phys. 25, 1-26 (1976).
\bibitem{coulouvrat} Coulouvrat, F.:
Synth\'ese bibliographique sur la focalisation de la d\'etonation
balistique, Rapport final de contrat A\'erospatiale 96/065,
Universit\`e Pierre et Marie Curie, Paris, (1997).
\bibitem{guiraud} Guiraud, J.-P.:
{Acoustique Geometrique, Bruit Balistique des Avions
Supersoniques et Focalisation,} J. Mec. 4, 215-267 (1965).
\bibitem{harten} Harten, A., Engquist, B., Osher, S. and Chakravathy,
S.: {Uniformly high order accurate essentially non-oscillatory
schemes, III}, Journal of Computational Physics 71, 231-303
(1987).
\bibitem{hunter} Hunter, J. K. and Keller, J. B.:
Caustics of nonlinear waves. Wave Motion, 9(5):429-443, (1987).
\bibitem{jiang} Jiang, G. and Shu, C.-W.:
{Efficient implementation of weighted ENO schemes},
Journal of Computational Physics 126, 202-228 (1996).
\bibitem{kandil2} Kandil, O.A. and Zheng, X:
Computational solution of Nonlinear Tricomi Equation for sonic
Boom Focusing and Applications, Paper Number 2005-43,
International Sonic Boom Forum, Penn State University, PA, (2005).
\bibitem{kandil3} Kandil, O. and Khasdeo, N.:
{Parametric Investigation of Sonic Boom Focusing using Solution
of Nonlinear Tricomi Equation,} AIAA Paper 2006-0415 (2006).
\bibitem{liu} Liu. X., Osher, S. and Chan, T.:
{Weighted essentially non-oscillatory schemes}, Journal of
Computational Physics 115, 200-212 (1994).
\bibitem{shu1} Shu, C.-W. and Osher, S.:
{Efficient implementation of essentially non-oscillatory
shock-capturing schemes}, Journal of Computational Physics
77 439-471 (1988).
\bibitem{shu2} Shu C.-W.: {Essentially non-oscillatory and
weighted essentially non-oscillatory schemes for
hyperbolic conservation laws}, in:
{\it Advanced Numerical Approximation of Nonlinear Hyperbolic
Equations}, Lecture Notes in Mathematics, vol. 1697,
Springer-Verlag, Berlin, pp. 325-432 (1998).
\end{thebibliography}
\end{document}