The definitions and physical units of the variables used in the PDE’s are given below. In the following table, \(L\,\ \)stands for length, \(T\) for time, \(M\) for mass and \(N\) for force.
| term | meaning |
dimensions | SI units |
| \(p\) | acoustic air pressure in medium |
\(\frac {N}{L^{2}}\) or \(\frac {ML}{T^{2}}\frac {1}{L^{2}}\) or \(\frac {M}{LT^{2}}\) | \(N/Meter^{2}\) |
| \(u\) | acoustic perturbation velocity |
\(L/T\) | \(Meter/Second\) |
| \(c\) | speed of sound in medium |
\(L/T\) | \(Meter/Second\) |
| \(K\) | bulk modulus or modulus of bulk elasticity for gas10 | \(\frac {M}{T^{2}L}\) | \(kg\) per \(meter\) per \(second^{2}\) |
| \(\rho \) | air density |
\(M/L^{3}\) | \(kg/meter^{3}\) |
To show that the system is hyperbolic, the PDE’s are written in matrix form\begin {align*} p_{t}+Ku_{x} & =0\\ \rho u_{t}+p_{x} & =0 \end {align*}
Therefore\begin {align*} \begin {pmatrix} p\\ u \end {pmatrix} _{t}+\overset {A}{\overbrace {\begin {pmatrix} 0 & K\\ 1/\rho & 0 \end {pmatrix} }}\begin {pmatrix} p\\ u \end {pmatrix} _{x} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \\ \mathbf {q}_{t}+A\mathbf {q}_{x} & =\mathbf {0} \end {align*}
If the eigenvalues of \(A\) are real and distinct, implying the existence of linearly independent eigenvectors for \(A\), then the system is called strictly hyperbolic11. The eigenvalues of \(A\) are found by solving the following equation\begin {align*} Det\left ( A-\lambda \mathbf {I}\right ) & =0\\ \left ( -\lambda \right ) \left ( -\lambda \right ) -\left ( k\right ) \left ( 1/\rho \right ) & =0\\ \lambda ^{2} & =\frac {k}{\rho }\\ \lambda _{1,2} & =\pm \sqrt {\frac {k}{\rho }} \end {align*}
The quantity \(\frac {k}{\rho }\)is positive and real because \(\rho \) is density (which is a real positive number) and \(k\) is bulk modulus of compressibility which is also real positive number.
Therefore both eigenvalues of \(A\) are real and distinct. Hence the system is strictly hyperbolic. The system is diagonalizable as well, since the transpose of \(A\) is a diagonal matrix, but this property was not needed to show the system is hyperbolic. The speed of sound in the medium is given by \(\sqrt {\frac {k}{\rho }}\). Hence a sound wave will travel in one direction at speed\(\sqrt {\frac {k}{\rho }}\) and another sound wave will travel in the same speed but in the opposite direction.
The following diagram illustrates the grid numbering used in the numerical solution
The Lax-Wendroff scheme for the linear system \(\mathbf {q}_{t}+A\mathbf {q}_{x}=\mathbf {0}\) is given by
\[ \mathbf {q}_{j}^{n}+1=\mathbf {q}_{j}^{n}-\frac {\Delta t}{2h}A\left ( \mathbf {q}_{j+1}^{n}-\mathbf {q}_{j-1}^{n}\right ) +\frac {\Delta t^{2}}{2h^{2}}A^{2}\left ( \mathbf {q}_{j-1}^{n}-2\mathbf {q}_{j}^{n}+\mathbf {q}_{j+1}^{n}\right ) \]
Where \(A=\begin {pmatrix} 0 & K\\ 1/\rho & 0 \end {pmatrix} \) is a constant matrix.
In this problem, the solution at time \(n\) is \[ \mathbf {q}_{j}^{n}=\begin {pmatrix} p\\ u \end {pmatrix} _{j}^{n}=\mathbf {q}_{j}^{n}=\begin {pmatrix} p_{j}^{n}\\ u_{j}^{n}\end {pmatrix} \] The following are the boundary conditions used \[ \mathbf {q}_{0}^{n}=\begin {pmatrix} p\\ u \end {pmatrix} _{0}^{n}=\begin {pmatrix} p_{1}^{n}\\ u_{1}^{n}\end {pmatrix} \]\[ \mathbf {q}_{N+1}^{n}=\begin {pmatrix} p\\ u \end {pmatrix} _{N+1}^{n}=\frac {1}{2}\begin {pmatrix} p_{N}^{n}+u_{N}^{n}\sqrt {k\rho }\\ \frac {p_{N}^{n}}{\sqrt {k\rho }}+u_{N}^{n}\end {pmatrix} \]
To find the time step \(\Delta t\), Courant number \(r=0.8\) was used12, and \(\Delta t\) found by using the CFL condition\[ r=\left \vert \frac {\Delta t}{h}\lambda \right \vert \] Solving for \(\Delta t\) gives \[ \Delta t=\frac {rh}{\left \vert \lambda \right \vert }\] The solution was implemented in Matlab and the result is given below. For each run, a number of plots are shown to illustrate the solution at different time instances. The following table describes the simulations done. Three different initial conditions are used with two different runs for each initial condition. The first run used the boundary conditions given in this problem, and the second run used different boundary conditions which caused the sound wave to reflect when it reached both the left and the right boundaries, and not just the left boundary. Therefore a total of 6 simulations were made, the first three used the following boundary conditions\begin {align*} p_{0}^{n} & =p_{1}^{n}\\ u_{0}^{n} & =-u_{1}^{n}\\ p_{N+1}^{n} & =\frac {1}{2}\left ( p_{N}^{n}+u_{N}^{n}\sqrt {k\rho }\right ) \\ u_{N+1}^{n} & =\frac {1}{2}\left ( \frac {p_{N}^{n}}{\sqrt {k\rho }}+u_{N}^{n}\right ) \end {align*}
And the second three simulations used the following boundary conditions \begin {align*} p_{0}^{n} & =p_{1}^{n}\\ u_{0}^{n} & =-u_{1}^{n}\\ p_{N+1}^{n} & =p_{N}^{n}\\ u_{N+1}^{n} & =-u_{N}^{n} \end {align*}
The images below show the three initial conditions for the pressure \(p\left ( x,0\right ) \). The initial velocity \(u\left ( x,0\right ) \) was set to zero for all simulations. The following section shows the simulation plots for each one of the 6 simulations. All snapshots were taken at the same time for each run in order to compare the results. All runs were made with the following parameters: \begin {align*} h & =0.005\ meter\\ \Delta t & =0.1278\ ms\\ Courant\ number\ & =0.8\\ maximum\ run\ time & =0.005\ \sec \end {align*}
Animations of these runs are available above (in HTML version only).
| \(\sin \left ( 10\pi x\right ) \text { from }x=0.4\text { to }x=0.6\) | \(\sin \left (20\pi x\right ) \text { from }x=0.4\text { to }x=0.6\) | triangle function |
|
|
|
|
Simulation using first initial data and reflect from left end only
This simulation used \(p\left ( x,0\right ) =\sin \left ( 10\pi x\right ) \) from \(x=0.4\) to \(x=0.6\). The pressure wave starts in the middle, and immediately starts to split into two smaller waves, each one became half the amplitude of the original wave. Each smaller wave traveled in opposite directions. The wave that reached the left boundary was reflected back while the wave that reached the right boundary was absorbed into the boundary. After the left wave reflected back and eventually reached the right boundary, it was also absorbed. This resulted in the original wave disappearing. As the left wave reflected from the left end, it also flipped upside down, such that the leading half of the wave remained with positive amplitude and the trailing half remained with the negative amplitude.
Simulation using second initial data and reflect from left end only
These images show the simulation result using \(p\left ( x,0\right ) =\sin \left ( 20\pi x\right ) \) from \(x=0.4\) to \(x=0.6\). Each frame is taken at the same time as the first simulation. The same result can be seen as described in the first simulation.
Simulation using third initial data and reflect from left end only
This simulation uses the triangle pulse as the initial data. Each frame is taken at the same time as the first simulation. The same result can be seen as was described in the first simulation.
Simulation using first initial data and reflecting from both ends
The following 3 simulations are a repeat of the first 3, but using boundary conditions that caused the pressure wave to reflect from both the left and the right boundaries. This resulted in the wave reflecting back and forth all the time. When both waves met again at the middle, the original wave form was reconstructed for a very short time but in an upside down form compared to its original form, and then the whole cycle was repeated. When the waves met again for the second time in the middle, the original wave was reconstructed again, but this time with the same shape it was at the initial time. This process continued again. Since there was no diffusion term present in the PDE, this cycle repeated for the duration of the simulation and no energy was lost. The times of each frame is the same as was used in all the previous simulations.
Simulation using second initial data and reflect from both ends
This simulation used \(p\left ( x,0\right ) =\sin \left ( 20\pi x\right ) \) from \(x=0.4\) to \(x=0.6\), but using boundary conditions that caused the pressure wave to reflect from both the left and the right boundaries. The same observation can be made as with the previous simulation.
Simulation using third initial data and reflect from both ends
This simulation used a triangle pressure wave as its initial data but using boundary conditions that caused the pressure wave to reflect from both the left and the right boundaries. The same observation can be made as with the previous simulation.
The boundary conditions given in the problem are\begin {align*} p_{0}^{n} & =p_{1}^{n}\\ u_{0}^{n} & =-u_{1}^{n}\\ p_{N+1}^{n} & =\frac {1}{2}\left ( p_{N}^{n}+u_{N}^{n}\sqrt {k\rho }\right ) \\ u_{N+1}^{n} & =\frac {1}{2}\left ( \frac {p_{N}^{n}}{\sqrt {k\rho }}+u_{N}^{n}\right ) \end {align*}
At the left most cell (cell 0), the acoustic perturbation velocity \(u\) is negative its value on the inside cell, therefore the average value of \(u\) right at the left edge (start of the physical domain) will be zero, as shown by the following diagram
Physically, this represent a barrier or a wall where perturbation velocity is zero at the wall resulting in deflection. Having zero velocity at the edge means that the momentum of the wave is zero at the left boundary. Since momentum is conserved, then it must have a direction which is opposite to what it was in the previous time step. This is similar to a ball hitting a perfectly elastic wall. For the pressure boundary conditions, having the acoustic pressure in the left most cell and the ghost cell being the same means that the pressure drop or gradient is zero between these two cells. Therefore, no sound will be transmitted through the boundary since sound is transmitted only due to presence of a pressure gradient between adjacent spatial points in the medium.
On the right side, when taking the average between the right-most cell and the ghost cell at the right results in\begin {align*} u_{right\_edge} & =\frac {3}{4}u_{N}^{n}+\frac {1}{4}\frac {p_{N}^{n}}{\sqrt {k\rho }}\\ p_{right\_edge} & =\frac {3}{4}p_{N}^{n}+\frac {1}{4}u_{N}^{n}\sqrt {k\rho } \end {align*}
Therefore, the perturbation velocity \(u\) at the right edge is no longer zero, but it has the same sign as the velocity at the right most cell. Physically this means the acoustic wave will continue to have momentum in the same direction and will not reflect. For the pressure, there exists now a pressure gradient, therefore sound will travel across the right boundary. Physically, this boundary can be thought of as a sound absorbing wall. (For example, a wall treated with special paint or covering).
Given a sequence \(u_{j}^{0}\) which is monotone in \(j\), we need to show that when a TVD scheme is applied to this sequence, the resulting sequence \(u_{j}^{n}\) is also monotone at any \(n\). This is the same as saying that a TVD is monotone preserving.
We are given that the sequence \(u_{j}^{n}\) has the fixed boundary conditions at \(j=\pm \infty \) for any \(n.\)
A monotone sequence can be either monotone increasing or monotone decreasing but not both. A monotone increasing sequence \(u\) is one where \(u_{j}\leq u_{j+k}\) for any \(j\) and for any \(k>j\). A monotone decreasing sequence is one where \(u_{j}\geq u_{j+k}\) for any \(j\) and for any \(k>j\). In the following discussion, a monotone sequence is taken to mean either an increasing or a decreasing sequence.
The following diagram illustrates this point. In this diagram the scheme is viewed as a system or an operator which transforms a sequence to a new sequence. We need to show that this transformation is monotone preserving when the operator is the TVD scheme.
Since \(u_{j}^{0}\) (the initial sequence) can be assumed to be monotone, then the total variation of \(u_{j}^{0}\) is known, which is\[ TV\left ( u_{j}^{0}\right ) =\left \vert U_{+\infty }-U_{-\infty }\right \vert \] The total variation is defined as the sum of the total amount the sequence change (in absolute values). In other words, the TV of the initial sequence is
\begin {align*} TV\left ( u^{0}\right ) & ={\displaystyle \sum \limits _{j}} \left \vert u_{j}^{0}-u_{j-1}^{0}\right \vert \\ & =\left \vert U_{+\infty }-U_{-\infty }\right \vert \end {align*}
\(TV\left ( u_{j}^{0}\right ) =\left \vert U_{+\infty }-U_{-\infty }\right \vert \) is valid since \(u^{0}\) is monotone. We could not have said this if \(u^{0}\) was not monotone. The following diagram helps illustrate why this is the case, showing a monotone sequence, and showing that adding all the differences between successive values in the sequence is the same as the difference between the left-most value and the right-most values (in absolute values).
The above is similar to walking up a staircase. If we are told that each step could only go up (or remain flat), then the total height of the overall staircase is the total variation, which is the sum of the height difference between each 2 successive steps.
We know that a TVD scheme, by definition, is one in which satisfies the following relation for any \(n\)\begin {equation} TV\left ( u^{n}\right ) \leq TV\left ( u^{0}\right ) \tag {1} \end {equation} We now need to show, that when \(u^{0}\) is monotone, then \(u^{n}\) will also be monotone when applied to a TVD scheme.
The proof will be by contradiction. The idea is to assume that the scheme is TVD, hence Eq. (1) is true, and then to assume that the scheme, when applied to an initial monotone sequence \(u^{0}\) has resulted in a sequence \(u^{n}\) which is no longer monotone. Then we show that this result is a contradiction to the assumption, meaning that \(u^{n}\) must be monotone.
The following proof below is for a monotone increasing sequence \(u^{0}\), but the same idea of the proof can be used for a monotone decreasing sequence.
Proof
Let the scheme be TVD, therefore \(TV\left ( u^{0}\right ) \leq TV\left ( u^{n}\right ) \), and let a monotone increasing sequence be \(u_{j}^{0}\) with a total variation \(TV\left ( u^{0}\right ) =\Delta \), where \(\Delta \) is some constant that does not change with \(n\). In this problem this constant is given as \(\left \vert U_{+\infty }-U_{-\infty }\right \vert \).
Let result of applying the TVD scheme to \(u_{j}^{0}\) be the sequence \(u_{j}^{n}.\) Now, assume that \(u_{j}^{n}\) is no longer a monotone increasing sequence. Since \(u_{j}^{n}\) is not monotone sequence, it must contain at least one local minimum and/or one local maximum. To illustrate this in a diagram, assume \(u_{j}^{n}\) had one local minimum. The same idea would apply if we assumed a local maximum.
Since \(u^{n}\ \)has a local minimum, then the total variation of \(u_{j}^{n}\) is now larger than the total variation of what it would had been if it did not have this local minimum. In the above diagram, \(u_{j}^{n}\) is shown as being monotone increasing, except for the one local minimum which appeared as a result of applying the TVD scheme.
Due to the presence of this local minimum, the total variation has become larger than \(\left \vert U_{+\infty }-U_{-\infty }\right \vert \). The extra amount added to \(TV\left ( u^{0}\right ) \) is seen as \(2\left \vert u_{j}^{n}-u_{j-1}^{n}\right \vert \), as this is the distance needed to be traversed in going down the local minimum and climbing back up the same level before meeting this local minimum.
Therefore, having a local minimum (or a local maximum) in a sequence increases its total variation. Therefore \(TV\left ( u^{n}\right ) >TV\left ( u^{0}\right ) \). However, we started by assuming that the scheme is TVD, which means that \(TV\left ( u^{n}\right ) \leq TV\left ( u^{0}\right ) \), so this result is a contradiction to our assumption.
Therefore \(u_{j}^{n}\) can not be a non monotone sequence, hence it must be a monotone sequence. This completes the proof.
The PDE \[ u_{t}+au_{x}=0 \] was solved using finite volume method using the 7 flux limiter functions listed in the problem statement above. The following tables summarize the observations made after running the simulations using each of these limiter functions. Each method was given a letter grade based on how close it was to the exact solution and how well the numerical solution appeared. Numerical solutions that showed ripples around the region of discontinuous data (corners) or showed more spatial lag relative to the exact solution, or had large amount of diffusion were graded lower than those which did not show any of these result.
| method | comment | grade | ||
| Upwinding | Very large diffusion seen at wave crest and trough, but no shift (lag).
| F | ||
| Lax-Wendroff |
| B- | ||
| Beam-Warming | Similar to Lax-Wendroff, but shift was to the right relative to exact solution. | B | ||
| minmod | Diffusion was present at wave crest and trough, but no shifting. | C | ||
| superbee | No shifting and very small amount of diffusion at crest and trough.
| B+ | ||
| MC | Similar to superbee, but a little more diffusion at crest and trough.
| B | ||
| Van Leer | Similar to MC limited, but much more diffusion at crest and trough.
| B- | ||
Among the high resolution limiter functions, superbee had the best numerical result.
| method | comment | grade | ||
| Upwinding | No shifting, but large amount of diffusion at crest and trough of the wave.
| C | ||
| Lax-Wendroff | No shifting and no diffusion.
| A | ||
| Beam-Warming | Very similar to Lax-Wendroff.
| A | ||
| minmod | No shifting, but small amount of diffusion was present near crest and trough. | B | ||
| superbee |
| A- | ||
| MC |
| A | ||
| Van Leer | No diffusion and no shifting
| A | ||
Among the high resolution methods, MC and Van Lee had the best results. Among the non high resolutions method, Lax-Wendroff and Beam-Warming were the best.
| method | comment | grade | ||
| Upwinding |
| C | ||
| Lax-Wendroff |
| C | ||
| Beam-Warming | The ripples are larger and have a larger extent than Lax-Wendroff. | C- | ||
| minmod | No ripples and little diffusion. An improved version of upwinding. | C+ | ||
| superbee |
| A- | ||
| MC | Similar to superbee, but more diffusion.
| B | ||
| Van Leer | Similar to MC limited.
| B+ | ||
Among the high resolution methods, superbee was the best. Among the non high resolutions method, Lax-Wendroff and Beam-Warming are best.
Among the high resolutions methods, I would choose superbee. It handled discontinues data the best and did well for smooth data, even though MC and Van Leer did a little better on the low frequency data, superbee had less diffusion in the wave packet data. So, overall, and in particular since it handled discontinues data better than any other flux limiter function, it is the method I would choose in practice.
Among the non high resolution methods, Lax-Wendroff and Beam-Warming were very similar. Upwinding did not do well. All the non high resolution methods did relatively worst in the step function test compared to the high resolution methods, as they were not able to handle solution near the discontinues regions as well as the high resolution methods did.
Numerical solutions using all the above methods have been animated and available to run at my course web page. All the animations run for 5 seconds each.
10The term \(k\) behaves like the stiffness of a spring. The larger the stiffness on the matrial, the faster a sound wave will travel in it. Hence sound speed in a perfectly incompressible fluid (which is not possible) will have an infinite speed.
11Another method to show that the system is hyperbolic, is to show that \(A\) is real and symmetric, because this implies that \(A\) is diagonalizable. In this case, the system is called symmetric hyperbolic.
12For stability, the Courant number must be less than 1.
