\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}