Page 331 - Computational Fluid Dynamics for Engineers

10.13 Model Problem for the Implicit Method: Unsteady Shock Tube 321

10.13 Model Problem for the Implicit Method:
Unsteady Shock Tube

The model problem is the shock tube problem of Section 10.10 and the initial
conditions (subsection 10.10.1) and boundary conditions (subsection 10.10.2)
are identical.

10.13.1 Solution Procedure and Sample Calculations

The computer program is given in Appendix B and, since the shock tube
problem was solved using the explicit algorithm of MacCormack in subsection
10.10.3, the numerical solution of the problem is constructed by borrowing from
the subroutines already developed. First, a 1-D grid with equidistant spacing is
generated (subroutine generate grid). The flux variables are then initialized
according to the specified boundary condition of subsection 10.9.1 for a diatomic
gas (subroutine i n i t i a l conditions). Then the algorithm marches in time,
s
with increments satisfying the stability restrictions (subroutine time tep). The
implicit steps are computed (subroutine implicit) to yield the updated flux
values. The program ends when the time has reached 250. The program uses the
same subroutine as in subsection 10.10.3 to compute the Euler fluxes F given
the primitive variables (subroutine flux).
The Jacobian matrix A is first constructed at all points for the implicit
steps and, since it operates on a 1st-derivative central difference operator, it is
computed for the lower and upper band of the block tridiagonal matrix. The
main band is formed by the identity matrix. The implicit artificial dissipation,
which is operated on by a 2nd-derivative central difference operator, is added
on the diagonal of each matrix: lower, upper and main diagonal. The explicit
artificial dissipation construction is straightforward, and is added to the right
hand side vector.
Zero-order extrapolation is applied to the smallest and largest vector element
of the right hand side as numerical boundary scheme, which effectively let the
waves escape from the computational domain. Once the left side tridiagonal
matrix formed, a standard block-tridiagonal subroutine is called to update the
primitive quantities.
The numerical results and the theoretical values are presented in Fig. 10.17
for several CFL number (0.5, 0.9, 2.0, 5.0) with the explicit artificial dissipation
coefficient held constant at 0.1 and the implicit artificial dissipation coefficient
set 2.5 times the explicit one. It can be seen that there are no dispersion er-
rors, but dissipative error remains as expected from the truncation error terms
operating on second derivatives. The dissipation increases as the CFL number
increased, and the implicit procedures allows CFL number greater than 1, which
is consistent with the stability analysis.