# Beyond Training: Optimizing Reinforcement Learning Based Job Shop Scheduling
Through Adaptive Action Sampling

Constantin Waubert de Puiseau, Christian Dörpelkus, Jannik Peters, Hasan
Tercan, Tobias Meisen

###### Abstract

Learned construction heuristics for scheduling problems have become
increasingly competitive with established solvers and heuristics in recent
years. In particular, significant improvements have been observed in solution
approaches using deep reinforcement learning (DRL). While much attention has
been paid to the design of network architectures and training algorithms to
achieve state-of-the-art results, little research has investigated the optimal
use of trained DRL agents during inference. Our work is based on the
hypothesis that, similar to search algorithms, the utilization of trained DRL
agents should be dependent on the acceptable computational budget. We propose
a simple yet effective parameterization, called $\delta$-sampling that
manipulates the trained action vector to bias agent behavior towards
exploration or exploitation during solution construction. By following this
approach, we can achieve a more comprehensive coverage of the search space
while still generating an acceptable number of solutions. In addition, we
propose an algorithm for obtaining the optimal parameterization for such a
given number of solutions and any given trained agent. Experiments extending
existing training protocols for job shop scheduling problems with our
inference method validate our hypothesis and result in the expected
improvements of the generated solutions.

## Introduction

Research on self-learning algorithms for scheduling problems has significantly
increased in recent years, driven primarily by advancements in algorithms and
the availability of more affordable and powerful computing units in both
industry and academia. In academia, the standardized job shop scheduling
problem (JSSP) (Pinedo 2016) is the dominant object of study for its
straightforward logical constraints and comparability thanks to public
benchmarks and libraries (van Hoorn 2018; Waubert de Puiseau et al. 2023).
While in some successful DRL-based solution methods agents steer existing
search algorithms (Zhang et al. 2022; Ni et al. 2021) and solvers (Tassel,
Gebser, and Schekotihin 2023), in a majority, agents autonomously generate
solution schedules imitating a construction heuristic (Zhang et al. 2020; Park
et al. 2021; van Ekeris, Meyes, and Meisen 2021; Tassel et al. 2022).

In these DRL-based construction heuristics, agents make iterative decisions to
schedule the next unscheduled operation until the entire schedule is complete.
For each decision, an action is sampled from the agent’s action vector that
consists of logits representing the current operation prioritization.
Depending on how actions are sampled from the vector, different solutions for
a JSSP are obtained. For our purposes, we differentiate between the following
sampling strategies:

  * •

Deterministic sampling: by greedily applying argmax policy logits we obtain
one deterministic solution per problem instance.

  * •

Stochastic sampling: sampling from the given distribution, where the
likelihood of each action is given by the corresponding logit.

  * •

Exploitative/Explorative sampling: during exploitative sampling, the
likelihoods of likely actions increase even more; during explorative sampling,
the likelihoods of all actions are assimilated. We achieve both through the
$\delta$-sampling which we present in this paper.

In their extremes, exploitative sampling converges towards deterministic
sampling whereas explorative sampling converges towards sampling from a
uniform distribution. Stochastic sampling falls in between. In the context of
scheduling, it has only very recently been shown that stochastic sampling can
produce better-than-greedy schedules even with a relatively small number of
sampled solutions (Iklassov et al. 2022). Motivated by this finding, our goal
is to identify the optimal sampling method for any given trained agent and the
desirable number of samples to be evaluated. During inference, the objective
is to sample in such a way that among all sampled solutions, the completion
times $C\textsubscript{max}$, also known as makespans, the shortest is as
short as possible. To avoid confusion with the common subscript $max$, and for
brevity, we refer to this minimum value as $C^{*}$ throughout this paper.
According to the proposed nomenclature, we formulate our working hypothesis as
follows:

  1. 1.

The expected $C^{*}$ is a function of the sampling strategy and number of
sampled solutions, i.e., the sample size.

  2. 2.

For a fixed (non-deterministic) sampling strategy, $C^{*}$ remains equal or
decreases with increasing sample size.

  3. 3.

For a small sample size, strongly exploitative sampling strategies
statistically lead to lower $C^{*}$ compared to very explorative strategies.
In other words, it is unlikely to find a local and better-than-greedy minimum
when deviating more strongly from the greedy solution.

  4. 4.

Conversely, for a large sample size, strongly explorative sampling strategies
statistically lead to better $C^{*}$ compared to exploitative strategies,
since a wider solution space is covered.

Figure 1 summarizes our hypothesis. Moreover, if the hypothesis is valid,
there should exist a sampling strategy that minimizes the expected $C^{*}$ for
a given agent and sample size.

Figure 1: Expected minimal makespans $C^{*}$ over sampling size for different
sampling strategies

Building upon our hypothesis, we introduce a new sampling method and validate
its effectiveness on agents trained with two different state-of-the-art
training protocols, testing them on JSSP instances of different sizes and with
different sample sizes. Our study leads to the following valuable
contributions:

  1. 1.

A new simple but effective tree-search-like solution generation algorithm for
JSSP based on action sampling

  2. 2.

An algorithm for the identification of the most effective parameterizations of
action sampling methods given a trained agent and a number of samples

  3. 3.

A validation of our algorithm on varios JSSP instances of different sizes as
well as on benchmark instances of Taillard (1993).

The remainder of this paper is structured as follows: After discussing related
work, we describe our methods and experiments. We then present and discuss the
experimental results, before drawing conclusions and giving an outlook on
future work.

## Related Work

DRL-based construction heuristics. Early successful attempts to solve
combinatorial optimization problems with DRL in 2016 (Bello et al. 2017) were
followed by a growing wave of research tackling the JSSP (Panzer, Bender, and
Gronau 2022). Recent research focuses on effective problem size agnostic
neural network architectures, including graph embeddings (Zhang et al. 2020;
Park et al. 2021) and recurrent networks (Iklassov et al. 2022). Another trend
is the use of curriculum learning (Iklassov et al. 2022; Waubert de Puiseau,
Tercan, and Meisen 2023). Interestingly, only Iklassov et al. (2022) sample
multiple solutions stochastically, although the possibility to obtain better
solutions through sampling in combinatorial problems was already established
by Kool, van Hoof, and Welling (2019). However, to the best of our knowledge,
our study is the first to focus on improving sampling methods for generating
JSSP solutions using trained agents.

Policy-guided tree search. Stochastic sampling from a learned policy is a
common concept to guide tree search algorithms (Orseau et al. 2018). The most
common search algorithm for this integration is Monte-Carlo tree search
(MCTS). Famously, this combination was leveraged in AlphaGo (Silver et al.
2016) and its successors. Therein, the expansion of nodes is determined by a
linear combination of learned Q-values or action predictions and the visit
count of states. Lately, the approach has been successfully transferred to
scheduling as well (Oren et al. 2021; Kumar and Dimitrakopoulos 2021).
However, MCTS comes at a high computational expense that is accepted in many
cases because it provides good expected values even in stochastic environments
that are infeasible to approximate analytically. For the deterministic JSSP,
however, the computational budget may be used more efficiently by sampling
multiple whole trajectories depth-wise from start to finish, as done in this
study.

Solution space coverage. The question of how closely to adhere to a learned
policy is captured in the exploration-exploitation tradeoff. However, this
tradeoff is usually exclusively considered during training. Therefore, methods
like epsilon-greedy sampling (Sutton and Barto 2018) or added loss terms on
the action confidence through entropy or curiosity terms (Schulman et al.
2017; Deepak Pathak et al. 2017) do not transfer to inference. A noteworthy
technique to fine-tune exploration during inference in MCTS is to use an
exponent across the action vector (Wang, Hao, and Cao 2021; Shen et al. 2018).
The primary rationale is a de- or increased effect of action priors compared
to visit counts, but a secondary effect is a (dis-)alignment of values in the
action vector itself. Our method is inspired by this idea and leverages the
secondary effect.

For combinatorial optimization problems, depth-wise search with DRL is mostly
done by stochastic policy sampling (Iklassov et al. 2022; Kool, van Hoof, and
Welling 2019; Orseau et al. 2018). One exception is beam-search, which ensures
a fixed tree width $k$ during search (Vinyals, Fortunato, and Jaitly 2015).
Known extensions make use of problem symmetries (Kwon et al. 2020) and entropy
in the action vector for the solution space expansion (Tassel et al. 2022).
Iklassov et al. (2022) compared stochastic sampling and beam-search
variations, finding that stochastic sampling outperforms the other methods.
However, they only reported results on a single sample size. We expand this
study with a new sampling method which fine-tunes sampling for each
combination of trained agent and desired sample size.

## Method and Experiments

###  $\delta$-Sampling

Our proposed sampling strategy offers the possibility to adjust the balance
between exploitative and explorative sampling. To this end, we manipulate the
action vector $a$ predicted by a trained agent to obtain $a_{sample}$ by
exponentiating the predicted logits with the discrimination exponent $\delta$
and then re-scale the vector as in equation 1.

$\displaystyle a_{sample}=\frac{1}{{\sum_{i}a_{i}^{\delta}}}\cdot a^{\delta}$
(1)

The effect is that for $\delta>1$, relative differences between logits are
emphasized favoring exploitation, whereas for $\delta<1$ differences are
decreased favoring exploration. For $\delta=1$, logits remain unchanged,
defaulting to stochastic sampling. We call this procedure $\delta$-sampling.
It easily scales to any problem size and allows to precisely adjust the
sampling behavior.

Given a base model and sample size, we are looking for $\delta^{*}$ that
minimizes $C^{*}$. We expect the functions of $C^{*}(\delta^{*})$ to be smooth
and have a single minimum based on our hypothesis. Hence, our algorithm starts
with a grid search that iteratively decreases the distance between the two
most promising candidate values for $\delta$ to neighboring candidates by 50%.
An example is illustrated in Figure 2. In the first iteration (i=0), we sample
with each candidate $\delta$ from a list of candidates and the desired sample
size on validation JSSP instances i$\delta$. The initial list of candidates is
a hyperparameter that has to be set. Sampling returns the average of $C^{*}$
across all i$\delta$ per candidate. Then, additional candidates are added
between the two most promising and their adjacent candidates, rounding to two
decimal places. If the most promising candidates lie at minimum or maximum
$\delta$ values, the initial interval is extended by half of the largest
distance between adjacent candidates. The function repeats until a specified
number of iterations or minimal distance between adjacent candidates is
reached.

Figure 2: Three iterations of searching the optimal $\delta$ value

### Experimental Setup

JSSP formalization. This work deals with the common JSSP formalization in
which $J$ jobs, consisting of $O$ operations are processed on $M$ machines.
Each job visits each machine exactly once in a predefined order, such that
$O=M$. Additionally, each operation has a specified processing time and may
not be interrupted once started (no preemption). Machines can only process one
job at once (no overlap) and setup-times are ignored. The notation of problem
sizes follows the structure $J$x$M$, e.g., a 20x15 instance consists of 20
jobs á 15 operations on 15 machines.

Base Models. We perform experiments on two base model designs that are adapted
from Zhang et al. (2020) and Iklassov et al. (2022), two competitive and
often-referenced representatives of DRL-based construction heuristics. We
refer to the base models as L2D (Zhang et al. 2020) and L2G (Iklassov et al.
2022) in reference to their publication titles. From L2D, we can extract the
openly accessible pre-trained graph neural network weights for each respective
problem size 6x6, 15x15 and 20x20, except for 100x20, where we resort to the
30x20 model as done in the original publication. The second base model, L2G,
features a neural network architecture that uses a set2set recurrent layer
(Vinyals, Bengio, and Kudlur 2016) to be problem size agnostic. Since no
trained models are publicly available, we retrain our base models according to
the reported protocol, one for each problem size, with the following
curricula: for the 15x15 and 6x6 JSSP models, we choose adaptive curriculum
learning and train on 6x6 and 15x15 JSSP instances. For the 20x20 and 100x20
base model, we use the RASCL method (Iklassov et al. 2022) and train on 6x6,
15x15, and 20x20 JSSP instances.

JSSP Instances. We test our sampling method on four problem sizes. 15x15,
20x20, and 100x20 are common benchmark sizes in the operations research
domain. 6x6 is also common in the literature since it is useful for the visual
interpretations of results in Gantt charts. To avoid overfitting on test data,
we create separate JSSP instances for training and for the search of suitable
sampling parameterizations. These Taillard-like instances are created with
separate random seeds, assigning random perturbations of machine sequences to
each job and sampling processing times uniformly from the interval [1, 99].

## Results

### Hypothesis Validation

To empirically test our hypotheses illustrated in Figure 1, we create 10,000
solutions ($C^{*}_{10k}$) on ten 6x6 JSSP instances with our L2G base model
and choose four sampling parameterizations:

  1. 1.

random sampling, which draws actions from a uniform distribution;

  2. 2.

$\delta=1$ as baseline of how sampling is usually performed;

  3. 3.

$\delta=0.05$ as extreme example of explorative sampling;

  4. 4.

$\delta=10$ as extreme example of exploitative sampling.

Figure 3 displays the averaged minimal makespan achieved over different sample
sizes (logarithmic scale). To make up for statistical outliers in small sample
sizes, we average the minimal makespans achieved over multiple sampling runs.
Using Python array slicing notation, this can be formulated as in Equation 2,

$\displaystyle\frac{1}{10,000/s}\sum_{i=0}^{10,000/s}min\left(C^{*}_{10k}[i:s(i+1)]\right),$
(2)

where s is the sample size.

Figure 3: Minimal makespans over sampling size for different sampling
strategies

In Figure 3, we divide the plot into four zones at interesting line
intersections. In zone A, which in this example only includes a sample size of
one, strongly exploitative sampling yields the best solutions. In zones B and
C, stochastic sampling with $\delta=1$ yields better makespans than strongly
exploitative or explorative sampling. For sampling sizes $<50$, strongly
explorative sampling with $\delta=0.05$ performs better than stochastic
sampling. The switch from zone B to C indicates, where explorative sampling
outperforms exploitative sampling. The results in Figure 3 thereby verify our
hypothesis.

base | problem size | sample size
---|---|---
|  | 32 |  128 |  512
L2G | 6x6 | 0.30 |  0.10 |  0.11
15x15 | 0.70 |  0.86 |  0.43
20x20 | 1.25 |  1.46 |  0.80
100x20 | 0.58 |  0.43 |  n.a.
L2D | 6x6 | 2.06 |  1.37 |  0.93
15x15 | 12.25 |  6.02 |  3.32
20x20 | 10.35 |  10.38 |  10.05
100x20 | 14.0 |  10.88 |  n.a.
Table 1: Best found delta values for different base models, problem sizes and
sample sizes Figure 4: Examples of results in the $\delta$ value search
algorithm with found minima

### Optimal $\delta$-Values

In the following we present the results of finding and using the optimal
parameterizations of our action sampling approach, adhering to the sample size
of 128 used in (Iklassov et al. 2022) for comparison and adding 32 and 512 as
significantly smaller and larger examplary sample sizes. We omitted
experiments on 100x20 JSSP with 512 samples, since we expect little additional
insights compared to the considerable required computational times for solving
large instances. Table 1 shows the results obtained from our parameter search
algorithm per base model, problem size, and sample size.

The results show two general trends. Firstly, in accordance with the above
hypothesis, larger $\delta$-values are favorable for smaller sample sizes.
Secondly, more exploitative sampling seems favorable with increasing problem
size. However, we also observe exceptions to these general trends. For L2G and
problem sizes 15x15 and 20x20, there is a small increase of found
$\delta$-values from sample size 32 to 128, before they decrease again for
512. Also for L2G, the trend of increasing $\delta$-values with increasing
problem sizes discontinues for problem size 100x20. Overall, the found
$\delta$-values of L2G are smaller in comparison to L2D. These discontinuing
trends, exceptions, and differences between base models lead us to believe
that analytical approaches to finding $\delta$-values are intractable and they
underscore the necessity to determine the sampling parameterization
empirically and individually for any trained model and sample size.

For a deeper analysis, we examine the function $C^{*}(\delta)$ which was
tested on the evaluation instances during $\delta$-value search. Noteworthy
representative examples of this function are presented in Figure 4. Across all
examples, the functions satisfy the assumed smoothness and the hypothesis of a
single minimum is valid, although some noise exists. Our iterative algorithm
to find these minima converges towards reasonable minima, which are marked by
the green triangles. Note that the small local minima, e.g. in Figure 4a) at
$\delta\approx 0.8$, could render gradient-based methods ineffective. However,
as Figure 4a) also exemplifies, $\delta$ can have a minimal overall effect,
and in cases such as Figure 4b), the function plateaus across large $\delta$
ranges. In light of these observations it is likely that the rare exceptions
to decreasing $\delta$-values with increasing sample size in L2G in Table 1,
which were discussed above, are relics of statistical noise.

sample size | base | problem size | ours | stochastic | deterministic | improvement
---|---|---|---|---|---|---
32 | L2G | 6x6  | 496.0 | 505.5 | 517.0 | 1.9% |
15x15 | 1310.4 | 1311.4 | 1364.5 | 0.1% |
20x20 | 1782.9 | 1782.3 | 1866.2 | 0.0% |
100x20  | 5759.6 | 5764.1 | 5867.3 | 0.1% |
L2D | 6x6  | 519.3 | 516.1 | 577.1 | -0.6% |
15x15 | 1406.9 | 1434.6 | 1534.7 | 1.9% |
20x20 | 1891.4 | 1957.4 | 2012.5 | 3.4% |
100x20  | 5865.8 | 5928.4 | 5996.7 | 1.1% |
128 | L2G | 6x6  | 491.7 | 503.6 | 517.0 | 2.4% |
15x15 | 1300.4 | 1300.9 | 1364.5 | 0.0% |
20x20 | 1769.3 | 1767.6 | 1866.2 | -0.1% |
100x20 | 5733.6 | 5745.1 | 5867.3 | 0.2% |
L2D | 6x6  | 504.8 | 506.2 | 577.1 | 0.3% |
15x15 | 1385.7 | 1410.1 | 1534.7 | 1.7% |
20x20 | 1865.8 | 1927.6 | 2012.5 | 3.2% |
100x20  | 5869.8 | 5899.7 | 5996.7 | 0.5% |
512 | L2G | 6x6  | 489.9 | 501.6 | 517.0 | 2.3% |
15x15 | 1281.6 | 1291.7 | 1364.5 | 0.8% |
20x20  | 1752.0 | 1750.6 | 1866.2 | -0.1% |
L2D | 6x6  | 498.0 | 498.7 | 577.1 | 0.1% |
15x15 | 1369.1 | 1386.9 | 1534.7 | 1.3% |
20x20  | 1842.2 | 1879.6 | 2012.5 | 2.0% |
Table 2: Average $C^{*}$ values on 100 generated test instances per problem size s.-size | base | p.-size | ours | stochastic | deterministic | improvement | optimal
---|---|---|---|---|---|---|---
|  |  | C* (opt. gap) | C* (opt. gap) | C* (opt. gap) | ours vs. stoch. | C
32 | L2G | 15x15 | 1345.9 (9.5%) | 1345.1 (9.5%) | 1417.4 (15.3%) | -0.1% |  | 1228.9
20x20 | 1858.9 (14.9%) | 1858.7 (14.9%) | 1935.3 (19.7%) | 0.0% |  | 1617.3
100x20 | 5759.7 (7.3%) | 5784.6 (7.8%) | 5878.4 (9.6%) | 0.4% |  | 5365.7
L2D | 15x15 | 1452.5 (18.2%) | 1465.9 (19.3%) | 1530.5 (24.5%) | 0.9% |  | 1228.9
20x20 | 1953.2 (20.8%) | 2027.1 (25.3%) | 2081.8 (28.7%) | 3.6% |  | 1617.3
100x20 | 5928.6 (10.5%) | 5983.1 (11.5%) | 6089 (13.5%) | 0.9% |  | 5365.7
128 | L2G | 15x15 | 1340.3 (9.1%) | 1340.9 (9.1%) | 1417.4 (15.3%) | 0.0% |  | 1228.9
20x20 | 1846.4 (14.2%) | 1837.4 (13.6%) | 1935.3 (19.7%) | -0.5% |  | 1617.3
100x20 | 5751.1 (7.2%) | 5766.5 (7.5%) | 5878.4 (9.6%) | 0.3% |  | 5365.7
L2D | 15x15 | 1416.2 (15.2%) | 1442 (17.3%) | 1530.5 (24.5%) | 1.8% |  | 1228.9
20x20 | 1943.4 (20.2%) | 2002.7 (23.8%) | 2081.8 (28.7%) | 3.0% |  | 1617.3
100x20 | 5888.2 (9.7%) | 5930.4 (10.5%) | 6089 (13.5%) | 0.7% |  | 5365.7
512 | L2G | 15x15 | 1323.1 (7.7%) | 1333.7 (8.5%) | 1417.4 (15.3%) | 0.8% |  | 1228.9
20x20 | 1825.5 (12.9%) | 1828.9 (13.1%) | 1935.3 (19.7%) | 0.2% |  | 1617.3
L2D | 15x15 | 1409.0 (14.7%) | 1427.6 (16.2%) | 1530.5 (24.5%) | 1.3% |  | 1228.9
20x20 | 1895.3 (17.2%) | 1946.2 (20.3%) | 2081.8 (28.7%) | 2.6% |  | 1617.3
Table 3: Results on Taillard Instances

### Performance Improvements

Test Instances. The effectiveness of the obtained parameterization is tested
on 100 test instances. The average $C^{*}$ values are presented in Table 2. We
report the results for $\delta$-sampling (ours), using stochastic sampling as
in Iklassov et al. (2022) (stochastic) as a baseline, both for three sample
sizes, and the deterministic solution (deterministic). Relative improvements
between ours and stochastic sampling ($(stochastic-ours)/stochastic$) are
depicted in the last column.

Both sampling methods find significantly better solutions than the
deterministic solution, confirming previous results (Iklassov et al. 2022;
Kool, van Hoof, and Welling 2019; Orseau et al. 2018). Our parameterized
method outperforms stochastic sampling in 17 out of 22 cases with up to 3.2%
improvement. Notably, improvements are achieved across all tested sample
sizes. The most notable exception is $\delta$ sampling for the L2D base model
on the 6x6 JSSP, which we partially attribute to statistical uncertainty given
the small sample size of 32. An exception that persists across sample sizes is
L2G on the 20x20 JSSP, where $\delta$-sampling performs slightly worse than
the stochastic baseline. Overall, we observe larger performance differences
where $\delta$ deviates more from one, i.e. stochastic sampling. This implies
that our parameter search is effective in finding stronger manipulations where
they are useful. On the other hand, constellations where $\delta$-values are
close to one and reside within plateaus of the function $C^{*}(\delta)$, as
shown in Figure 4 c) and d), only minor performance differences can be
obtained, as observed for L2G on the 20x20 JSSP.

Figure 5: Comparison of sampling methods over sample sizes

To understand those cases where very similar average $C^{*}$ are generated, we
analyze the resulting schedule differences. Within the schedules of L2G on the
15x15 JSSP with 128 samples, the most striking example, we counted in how many
cases the same schedules, measured by the same $C^{*}$, were found per
instance. Surprisingly, $C^{*}$ was equal in only 4/100 test instances. This
shows that even though the average values of $C^{*}$ were less than $0.00\%$
apart, different solutions were found, such that $\delta$ still had a
significant impact. In contrast, for 6x6 JSSPs of L2G with a sample size of
128, we obtained the same $C^{*}$ in roughly one third of the instances. This
is most likely an effect of the smaller total solution space compared to the
other problem sizes. The conclusion is that even if only a few different
solutions are found, they can still have a significant impact, as the $1.9\%$
improvement demonstrates. Therefore, our sampling method is effective for
small problem sizes with smaller solution spaces, too.

Benchmark Instances. We further report the results on benchmark instances of
Taillard (1993), which we abbreviate as TA instances, in Table 3. In
comparison, there are only ten instances per size, but the optimal solutions
are known, enabling us to additionally provide the relative gap (opt. gap) for
reference. The effects of sampling methods are fundamentally similar. As
before, both sampling approaches outperform the deterministic solution. While
the results on L2D indicate even greater improvements through our sampling
method, the results on L2G are more ambiguous. We attribute the differences of
the above results to the fact that only ten TA instances exist per problem
size, introducing a larger statistical error. Note that overall we find
slightly larger makespans than those reported in (Iklassov et al. 2022), as we
could not reproduce the reported performance on benchmarks with our trained
model. However, our contribution is a relative improvement per model and
sample size applicable to any learned construction heuristic, not the creation
of a new learning strategy or network architecture. For details on the most
recent absolute results and a comprehensive comparison of learning-based,
priority-rule-based or heuristic methods, we refer the interested reader to
(Iklassov et al. 2022), (Corsini et al. 2024) and (Falkner et al. 2022).

In Figure 5, relative improvements on TA instances through our adjusted
sampling method can be compared to improvements achievable by increasing
sample sizes. In many cases, we find that the achieved $C^{*}$ through
$\delta$-sampling is similar or better than the $C^{*}$ achieved via
stochastic sampling with a four times larger sample size. In practical terms,
the same results can be achieved with a quarter of the computational budget
through better sampling in these cases, rendering $\delta$-sampling four times
more efficient.

## Conclusion and Outlook

In this paper, we proposed a sampling method for the more effective usage of
trained DRL-based construction heuristics for the JSSP. We showed that an
optimal balance between exploration and exploitation during inference exists,
which is specific to the trained model and the chosen sample size. Based on
this insight, a method was developed that allows the tradeoff between
exploration and exploitation to be easily and smoothly parameterized during
solution generation of learned construction heuristics. The effectiveness of
this method was experimentally evaluated and showed promising results across
different base models and computational budgets on common JSSP sizes. We
believe that our approach is neither limited to the studied base models nor
the application to the JSSP, but can be used to improve any other learned
construction heuristic on both the JSSP and other related combinatorial
optimization problems.

In the future, we would like to transfer the insights gained into training
protocols that steer agents to find the best possible solutions given a
previously set computational budget to enforce a learned exploration-
exploitation tradeoff. In addition, we will combine our sampling method with
Monte Carlo Tree Search (MCTS) for a more effective policy-guided search
during inference.

## References

  * Bello et al. (2017) Bello, I.; Pham, H.; Le V, Q.; Norouzi, M.; and Bengio, S. 2017.  Neural Combinatorial Optimization with Reinforcement Learning.  _ICLR Workshop_.  http://arxiv.org/pdf/1611.09940v3.
  * Corsini et al. (2024) Corsini, A.; Porrello, A.; Calderara, S.; and Dell’Amico, M. 2024.  Self-Labeling the Job Shop Scheduling Problem.  arXiv:2401.11849.
  * Deepak Pathak et al. (2017) Deepak Pathak; Pulkit Agrawal; Alexei A. Efros; and Trevor Darrell. 2017.  Curiosity-driven Exploration by Self-supervised Prediction.  _International Conference on Machine Learning_ , 34(70): 2778–2787.
  * Falkner et al. (2022) Falkner, J. K.; Thyssens, D.; Bdeir, A.; and Schmidt-Thieme, L. 2022.  Learning to Control Local Search for Combinatorial Optimization.  _Falkner, Jonas K., et al. Learning to Control Local Search for Combinatorial Optimization. Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , 13717: 361–376.
  * Iklassov et al. (2022) Iklassov, Z.; Medvedev, D.; Solozabal, R.; and Takac, M. 2022.  Learning to generalize Dispatching rules on the Job Shop Scheduling.  _Advances in Neural Information Processing Systems_ , 33: 1621–1632.
  * Kool, van Hoof, and Welling (2019) Kool, W.; van Hoof, H.; and Welling, M. 2019.  Attention, Learn to Solve Routing Problems!  _International Conference on Learning Representations_ , 7.
  * Kumar and Dimitrakopoulos (2021) Kumar, A.; and Dimitrakopoulos, R. 2021.  Production scheduling in industrial mining complexes with incoming new information using tree search and deep reinforcement learning.  _Applied Soft Computing_ , 110: 107644.
  * Kwon et al. (2020) Kwon, Y.-D.; Choo, J.; Kim, B.; Yoon, I.; Gwon, Y.; and Min, S. 2020.  Pomo: Policy optimization with multiple optima for reinforcement learning.  _Advances in Neural Information Processing Systems_ , 34: 21188–21198.  https://arxiv.org/pdf/2010.16011.pdf.
  * Ni et al. (2021) Ni, F.; Hao, J.; Lu, J.; Tong, X.; Yuan, M.; Duan, J.; Ma, Y.; and He, K. 2021.  A Multi-Graph Attributed Reinforcement Learning based Optimization Algorithm for Large-scale Hybrid Flow Shop Scheduling Problem.  In Zhu, F., ed., _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, ACM Digital Library, 3441–3451. New York,NY,United States: Association for Computing Machinery.  ISBN 9781450383325.
  * Oren et al. (2021) Oren, J.; Ross, C.; Lefarov, M.; Richter, F.; Taitler, A.; Feldman, Z.; Daniel, C.; and Di Castro, D. 2021.  SOLO: Search Online, Learn Offline for Combinatorial Optimization Problems.  _Proceedings of the Fourteenth International Symposium on Combinatorial Search (SoCS 2021)_.
  * Orseau et al. (2018) Orseau, L.; Lelis, L. H. S.; Lattimore, T.; and Weber, T. 2018.  Single-Agent Policy Tree Search With Guarantees.  _32nd Conference on Neural Information Processing Systems (NIPS_ , 32: 3201–3211.
  * Panzer, Bender, and Gronau (2022) Panzer, M.; Bender, B.; and Gronau, N. 2022.  Neural agent-based production planning and control: An architectural review.  _Journal of Manufacturing Systems_ , 65: 743–766.
  * Park et al. (2021) Park, J.; Chun, J.; Kim, S. H.; Kim, Y.; and Park, J. 2021.  Learning to schedule job-shop problems: representation and policy learning using graph neural network and reinforcement learning.  _International journal of production research_ , 59(11): 3360–3377.
  * Pinedo (2016) Pinedo, M. 2016.  _Scheduling: Theory, algorithms, and systems_.  Springer International Publishing, fifth edition edition.  ISBN 978-3-319-26578-0.
  * Schulman et al. (2017) Schulman, J.; Wolski, F.; Dhariwal, P.; Radford, A.; and Klimov, O. 2017.  Proximal Policy Optimization Algorithms.  http://arxiv.org/pdf/1707.06347v2.
  * Shen et al. (2018) Shen, Y.; Chen, J.; Huang, P.-S.; Guo, Y.; and Gao, J. 2018.  M-Walk: Learning to Walk over Graphs using Monte Carlo Tree Search.  https://arxiv.org/pdf/1802.04394.pdf.
  * Silver et al. (2016) Silver, D.; Huang, A.; Maddison, C. J.; Guez, A.; Sifre, L.; van den Driessche, G.; Schrittwieser, J.; Antonoglou, I.; Panneershelvam, V.; Lanctot, M.; Dieleman, S.; Grewe, D.; Nham, J.; Kalchbrenner, N.; Sutskever, I.; Lillicrap, T.; Leach, M.; Kavukcuoglu, K.; Graepel, T.; and Hassabis, D. 2016.  Mastering the game of Go with deep neural networks and tree search.  _Nature_ , 529(7587): 484–489.
  * Sutton and Barto (2018) Sutton, R. S.; and Barto, A. 2018.  _Reinforcement learning: An introduction_.  Adaptive computation and machine learning. Cambridge, Massachusetts, London, England: The MIT Press, second edition edition.  ISBN 9780262039246.
  * Taillard (1993) Taillard, E. 1993.  Benchmarks for basic scheduling problems.  _European Journal of Operational Research_ , 64(2): 278–285.
  * Tassel, Gebser, and Schekotihin (2023) Tassel, P.; Gebser, M.; and Schekotihin, K. 2023.  An End-to-End Reinforcement Learning Approach for Job-Shop Scheduling Problems Based on Constraint Programming.  _Proceedings of the Thirty-Third International Conference on Automated Planning and Scheduling (ICAPS 2023)_ , 614–622.
  * Tassel et al. (2022) Tassel, P.; Kovács, B.; Gebser, M.; Schekotihin, K.; Kohlenbrein, W.; and Schrott-Kostwein, P. 2022.  Reinforcement Learning of Dispatching Strategies for Large-Scale Industrial Scheduling.  _Proceedings of the International Conference on Automated Planning and Scheduling_ , 32: 638–646.
  * van Ekeris, Meyes, and Meisen (2021) van Ekeris, T.; Meyes, R.; and Meisen, T. 2021.  Discovering Heuristics And Metaheuristics For Job Shop Scheduling From Scratch Via Deep Reinforcement Learning.  _Proceedings of the Conference on Production Systems and Logistics : CPSL 2021_ , 709–718.
  * van Hoorn (2018) van Hoorn, J. J. 2018.  The Current state of bounds on benchmark instances of the job-shop scheduling problem.  _Journal of Scheduling_ , 21(1): 127–128.
  * Vinyals, Bengio, and Kudlur (2016) Vinyals, O.; Bengio, S.; and Kudlur, M. 2016.  Order Matters: Sequence to sequence for sets.  _4th International Conference on Learning Representations 2016_.
  * Vinyals, Fortunato, and Jaitly (2015) Vinyals, O.; Fortunato, M.; and Jaitly, N. 2015.  Pointer Networks.  _Advances in Neural Information Processing Systems_ , 29: 2692–2700.
  * Wang, Hao, and Cao (2021) Wang, Q.; Hao, Y.; and Cao, J. 2021.  Learning to traverse over graphs with a Monte Carlo tree search-based self-play framework.  _Engineering Applications of Artificial Intelligence_ , 105: 104422.
  * Waubert de Puiseau et al. (2023) Waubert de Puiseau, C.; Peters, J.; Dörpelkus, C.; Tercan, H.; and Meisen, T. 2023.  schlably: A Python framework for deep reinforcement learning based scheduling experiments.  _SoftwareX_ , 22: 101383.
  * Waubert de Puiseau, Tercan, and Meisen (2023) Waubert de Puiseau, C.; Tercan, H.; and Meisen, T. 2023.  Curriculum Learning in Job Shop Scheduling using Reinforcement Learning.  _Proceedings of the Conference on Production Systems and Logistics: CPSL 2023_ , 34–43.
  * Zhang et al. (2020) Zhang, C.; Song, W.; Cao, Z.; Zhang, J.; Tan, P. S.; and Xu, C. 2020.  Learning to Dispatch for Job Shop Scheduling via Deep Reinforcement Learning.  _Advances in Neural Information Processing Systems_ , 33: 1621–1632.
  * Zhang et al. (2022) Zhang, C.; Song, W.; Cao, Z.; Zhang, J.; Tan, P. S.; and Xu, C. 2022.  Learning to Search for Job Shop Scheduling via Deep Reinforcement Learning.  https://arxiv.org/pdf/2211.10936.


# An analytical solution for supersonic flow over a circular cylinder using an
optimized shock shape

###### Abstract

An analytical solution for high supersonic flow over a circular cylinder based
on Schneider’s inverse method [1] has been presented. In the inverse method, a
shock shape is assumed and the corresponding flow field and the shape of the
body producing the shock are found by integrating the equations of motion
using the stream function. A shock shape theorised by Moeckel [2] has been
assumed and it is optimized by minimising the error between the shape of the
body obtained using Schneider’s method and the actual shape of the body. A
further improvement in the shock shape is also found by using the Moeckel’s
shock shape in a small series expansion. With this shock shape, the whole flow
field in the shock layer has been calculated using Schneider’s method by
integrating the equations of motion. This solution is compared against a fifth
order accurate numerical solution using the discontinuous Galerkin method
(DGM) and the maximum error in density is found to be of the order of
$10^{-3}$ which demonstrates the accuracy of the method used for both plane
and axisymmetric flows.

Keywords: cylinder, supersonic flow, shock wave, shock shape, inverse method,
discontinuous Galerkin method

##  1 Introduction

Analytical methods for hypersonic flows past blunt bodies in plane and
axisymmetric flows have been investigated quite extensively in the past. These
analyses can be found in papers like those of Chester [3], [4] and Freeman [5]
and in the books by Hayes and Probstein [6] and Rasmussen [7]. There are two
different problems that are examined by them. The first one is the so called
direct problem where the shape of the body is given and the shock shape and
the flow field are to be found. The second one is the inverse problem where
the shape of the shock wave is given, and it is required to find the
corresponding flow field and the shape of the body producing the shock. We are
interested in the inverse problem and particularly the method devised by
Schneider in [1] for flow past blunt bodies to solve the problem. We use the
shock shape theorised by Moeckel [2] for blunt bodies along with the method of
Schneider [1] to find an accurate analytical solution for plane flow past a
circular cylinder for high supersonic velocities. In a way, this procedure is
similar to the shock fitting technique [8] used by numerical solvers where the
shock shape is assumed and improved further based on the solution.

The analysis developed by Schneider allows a very elegant treatment of the
inviscid hypersonic blunt body problem. The fundamental advantage of this
method is its applicability in the flow field from the stagnation region up to
large distances from the nose of the body. This method for the solution in
shock layer is based on two main assumptions. First, it is assumed that the
density immediately behind a strong shock is much larger than in front of the
shock. The second assumption is based on the pressure along a streamline. With
these two assumptions, Schneider obtained the solution by integrating the
equations of motion using the stream function. The flow properties immediately
behind the shock are obtained using the Rankine-Hugoniot equations. As a
result of these assumptions and the subsequent solution, it is not necessary
to have a thin shock layer for this solution to hold. Due to these advantages,
this method has been used extensively in various aerospace applications like
in [9], [10], and [11]. Recently, this method has also been used for
astrophysical calculations in [12] and [13]. This method was also extended to
three-dimensional flows by Schwarze in [14]. Though this method was developed
for hypersonic flows, we show that this works well even in the high supersonic
flow regime.

Our method of solution is as follows. We assume a shock shape in the form
given by Moeckel in [2] as a single parameter family of hyperbolas. We
determine that parameter by minimising the error between the shape of the body
obtained using Schneider’s method and the actual shape of the body. This gives
a better shock shape $f$. Now, we improve this shock shape further by writing
it as a series expansion ($af+bf^{2}$) with unknown coefficients. We determine
these coefficients using the same optimization and get a much better
approximation for the shock shape. We use this shock shape to find the
solution of the flow in the shock layer and compare the solution obtained with
a numerical solution calculated using a fifth order accurate discontinuous
Galerkin method (DGM) and overset grids where the shock has been captured
within a grid line [15]. We find that the maximum error in density is of the
order of $10^{-3}$ (for the shock shape $af+bf^{2}$) demonstrating the
accuracy of the solution method. We also found that using more terms in the
series expansion for the shock shape (like $af+bf^{2}+cf^{3}$,
$af+bf^{2}+cf^{3}+df^{4}$) does not improve the accuracy of the solution.

The paper is organized as follows. We describe the formulation of the
Schneider’s method used for all our results in Section 2, approximation of the
shock shape is described in Section 3, the results obtained are given in
Section 4 and we conclude the paper in Section 5. We also give the formulation
and validation of the numerical method used (discontinuous Galerkin method) in
Appendix A.

##  2 Formulation

We now describe the formulation of Schneider’s method in brief. We have mostly
followed the notation used in [12] for the description of this method for the
plane or axisymmetric supersonic flow past a blunt body. The shape of the bow
shock wave (and hence the shock angle $\beta$) is given, and it is required to
find the corresponding flow field and the shape of the body producing the
shock. Let the flow be inviscid and without heat conduction, and let the gas
be perfect. $z$ and $r$ are the Cartesian coordinates for plane flow, and also
the cylindrical coordinates for axisymmetric flow. We also use a shock
oriented curvilinear coordinate system, where $x$ is the distance along the
shock surface in the plane formed by the shock normal and the direction of the
uniform fluid flow, and $y$ is the distance normal to the shock surface as
shown in Figure 1. Consider an arbitrary point Q in the shock layer. A
streamline $\psi=c$ passes through this point. The streamline intersects the
shock at the point S. The point N on the shock is determined such that a
normal to the shock at N also intersects the streamline, $\psi=c$, at the
point Q. We can now determine the $x$ and $y$ coordinates of Q and the
corresponding velocity components $u$ and $v$. The $z$-axis is taken to be
parallel to the direction of the incident flow. We emphasize that this is only
assumed for convenience, but is not a general restriction. It can be deduced
from Figure 1 that

$z=\hat{z}+y\sin\hat{\beta}$ (1)

$r=\hat{r}-y\cos\hat{\beta}$ (2)

where $\hat{\beta}$ is the shock inclination angle at the point N ($\hat{z}$,
$r$). The flow quantities immediately behind the shock at the point N are
denoted by a hat ($\wedge$), and at the point S by an asterisk (*).
Undisturbed flow quantities far upstream are denoted by the subscript
$\infty$.

O$x$ShockdxBodyQ$\hat{\beta}$N$\beta_{*}$S$y$$u$$v$$\psi=c$$r$$\hat{r}$$\hat{z}$$z$$r_{*}$$U_{\infty}$

Figure 1: The shock, body and a streamline are shown with a shock oriented
coordinate system. Here, $z$ and $r$ are the Cartesian coordinates for plane
flow with origin O, $x$ and $y$ are the distances along the shock surface and
normal to it and Q is an arbitrary point in the shock layer. $u$ and $v$ are
the velocities in $x$ and $y$ directions. $U_{\infty}$ is the freestream
velocity. $\beta$ is the shock angle. $\psi$ is the stream function. Flow
quantities immediately behind the shock at the point N (point where a normal
to the shock intersects Q) are denoted by a hat ($\wedge$), and at the point S
(intersection of the streamline through Q with the shock) by an asterisk (*).
This diagram is adapted from Schneider’s paper [1].

Since the functions $\hat{z}(x)$ and $\hat{r}(x)$ are known for a given shock
shape, equations (1) and (2) can be used to calculate the coordinates $z$ and
$r$ of a point Q from its coordinates $x$ and $y$. The curvature of the shock
contour at the point N is denoted by $\hat{\kappa}(x)$, defined as positive
when the surface is concave on the side of positive $y$.
Now, the continuity equation can be written as

$\frac{\partial(r^{j}\rho u)}{\partial
x}+\frac{\partial[(1-\hat{\kappa}y)r^{j}\rho v]}{\partial y}=0$ (3)

where $\rho$ is the fluid density. The parameter j is 0 for plane flow and 1
for axisymmetric flow. We now define a stream function $\psi$ which satisfies
(3) by

$\frac{\partial\psi}{\partial x}=(1-\hat{\kappa}y)r^{j}\rho v,$ (4)

$\frac{\partial\psi}{\partial y}=-r^{j}\rho u$ (5)

The stream function defined by (4) and (5) contains a constant of integration.
This constant is chosen such that $\psi=0$ is the body stream line. Then the
stream function represents the mass flow between the streamline $\psi$ =
constant. and the surface of the projectile, per unit depth for plane flows,
and per unit azimuthal angle (in radians) for axisymmetric flows. At the point
N, the stream function therefore is

$\hat{\psi}_{N}=\rho_{\infty}U_{\infty}\frac{\hat{r}^{1+j}}{1+j}$ (6)

Also, $\psi$ is connected to the coordinate $r_{*}$ of the point S by

$\psi_{*}=\rho_{\infty}U_{\infty}\frac{r_{*}^{1+j}}{1+j}$ (7)

We now introduce a new coordinate system with $\psi$ and $\bar{x}=x$ as
variables. With these new coordinates, the equations of energy, entropy, and
momentum conservation can written as:

$u^{2}+v^{2}+2h=u_{*}^{2}+v_{*}^{2}+2h_{*}=\text{constant}.,$ (8)

$\frac{\partial S}{\partial\bar{x}}=0\quad\text{or}\quad S=S_{*}(\psi)$ (9)

$u\frac{\partial u}{\partial\bar{x}}+v\frac{\partial
v}{\partial\bar{x}}+\frac{1}{\rho}\frac{\partial P}{\partial\bar{x}}=0$ (10)

$(1-\hat{\kappa}y)(1-j\frac{y}{\hat{r}}\cos\hat{\beta})\hat{r}^{j}\frac{\partial
P}{\partial\psi}=\hat{\kappa}u+\frac{\partial v}{\partial\bar{x}}$ (11)

where $P$ is the fluid pressure, $S$ the entropy, and $h$ is the specific
enthalpy.
The variable $y$ is now a dependent variable and we now get the transformation
equations using (4) and (5) as

$\frac{\partial y}{\partial\bar{x}}=(1-\hat{\kappa}y)\frac{v}{u}$ (12)

$\frac{\partial
y}{\partial\psi}=-\frac{1}{[1-j(y/\hat{r})\cos\hat{\beta}]\hat{r}^{j}\rho u}$
(13)

Now, the equations (10), (11), (12), (13), and (8) (or (9)) provide a set of
five equations for the five unknown dependent variables $u$, $v$, $y$, $P$ and
$h$ (or $S$).

The flow quantities immediately behind the shock may be obtained from the
Rankine-Hugoniot jump conditions in terms of the inverse compression ratio
across the shock $\hat{\chi}=\rho_{\infty}/\hat{\rho}$. They are

$\hat{u}=U_{\infty}\cos\hat{\beta}$ (14)

$\hat{v}=U_{\infty}\hat{\chi}\sin\hat{\beta}$ (15)

$\hat{P}=P_{\infty}+\rho_{\infty}U_{\infty}^{2}(1-\hat{\chi})\sin^{2}\hat{\beta}$
(16)

$\hat{h}=h_{\infty}+\frac{1}{2}U_{\infty}^{2}(1-\hat{\chi}^{2})\sin^{2}\hat{\beta}$
(17)

These equations maintain their validity even if the hats are replaced by
asterisks.
Schneider’s method for the solution in shock layer is based on two main
assumptions. First, it is assumed that the density immediately behind a strong
shock is much larger than in front of the shock; i.e.,

$\hat{\chi}=\frac{\rho_{\infty}}{\rho}\ll
1\quad\text{and}\quad\chi_{*}=\frac{\rho_{\infty}}{\rho_{*}}=O(\hat{\chi})$
(18)

Second, the pressure $P$ at the point Q of the disturbed flow field is not
much smaller than the pressure $\hat{P}$ at the point N (see Figure 1); i.e.,

$\frac{\hat{P}}{P}=O(1)\quad(\text{on }x=\text{const., }y>0)$ (19)

The symbol $f(x)=O(g(x))$ means that $|f(x)|$ is not very large in comparison
with $|g(x)|$. Using these assumptions, Schneider in [1] obtained the pressure
at an arbitrary point Q approximately as

$P=\hat{P}-\frac{\hat{\kappa}}{\hat{r}^{j}}\int_{\psi}^{\hat{\psi}}\left[u_{*}^{2}+2\left(h_{*}-h(\hat{P},S_{*})\right)\right]^{1/2}d\psi^{\prime}$
(20)

Here, we can note that all quantities on the right hand side of equation (20)
are given by the boundary conditions at the shock or by the equation of state.
The terms excluded in making the approximation are coming from the stagnation
region, as well as the region near the stagnation region where $u\gg\hat{u}$ –
are of the order of $\hat{\chi}$ and hence they contribute only negligibly to
the integral in equation (20). Therefore, the whole state of the gas is known
in the streamline coordinate system $(\bar{x},\psi)$ after evaluating
$S=S_{*}(\psi)$ and $P$ from equation (20). The actual location of the body in
space can be determined by solving the differential equation (13). We do this
by separation of variables giving the distance from the shock surface $y$ as a
function of $\bar{x}$ and $\psi$ as

$y\left(1-\frac{j\cos\hat{\beta}}{2\hat{r}}y\right)=\frac{1}{\hat{r}^{j}}\int_{\psi}^{\hat{\psi}}\frac{d\psi^{\prime}}{\rho
u}$ (21)

Neglecting errors of $O(\chi)$ and following [6], we may replace the velocity
component $u$ by

$u^{2}=u_{*}^{2}+2[h_{*}-h(P,S_{*})]$ (22)

This is obtained from the energy equation (8). The integral in equation (21)
can now be written as

$Y=\int_{\psi}^{\hat{\psi}}\frac{d\psi^{\prime}}{\rho(P,S_{*})[u_{*}^{2}+2(h_{*}-h(P,S_{*}))]^{1/2}}$
(23)

Solving the quadratic equation in $y$ on the left hand side of equation (21),
we have to distinguish between plane ($j=0$) and axisymmetric ($j=1$) flows.
This gives

$\text{for }j=0:\quad y=Y;$ (24)

$\text{for }j=1:\quad
y=\frac{\hat{r}}{\cos\hat{\beta}}\left[1-\left(1-\frac{2Y\cos\hat{\beta}}{\hat{r}^{2}}\right)^{1/2}\right]$
(25)

Equations (23), (24), and (25) give us the required results in $\bar{x}$ and
$\psi$ coordinates. We can now use the transformations (1) and (2) to get the
required results in $z$ and $r$ coordinates.
For a perfect gas with constant specific heats, the inverse compression ratio
is given by

$\chi_{*}=\frac{\gamma-1}{\gamma+1}+\frac{2}{(\gamma+1)M_{\infty}^{2}\sin^{2}\beta_{*}}$
(26)

We point out that an analogous relation is valid for $\hat{\chi}$, if all
asterisks are replaced by hats in equation (26).
Using the shock conditions (14)-(17), together with (26), the two integrals
(20) and (23), which have to be evaluated, become

$P=\hat{P}-\frac{U_{\infty}\hat{\kappa}}{\hat{r}^{j}}\int_{\psi}^{\hat{\psi}}\left[\cos^{2}\beta_{*}+\left(\frac{2}{(\gamma-1)M_{\infty}^{2}}+\sin^{2}\beta_{*}\right)\left(1-\left(\frac{\sin^{2}\hat{\beta}}{\sin^{2}\beta_{*}}\right)^{\frac{\gamma-1}{\gamma}}\right)\right]^{1/2}d\psi^{\prime}$
(27)

$Y=\frac{1}{\rho_{\infty}U_{\infty}}\int_{\psi}^{\hat{\psi}}\chi_{*}\left(\frac{\hat{P}\sin^{2}\beta_{*}}{P\sin^{2}\hat{\beta}}\right)^{1/\gamma}\left[\cos^{2}\beta_{*}+\left(\frac{2}{(\gamma-1)M_{\infty}^{2}}+\sin^{2}\beta_{*}\right)\left(1-\left(\frac{P\sin^{2}\hat{\beta}}{\hat{P}\sin^{2}\beta_{*}}\right)^{\frac{\gamma-1}{\gamma}}\right)\right]^{-1/2}d\psi^{\prime}$
(28)

The curvature of a curve in space can be calculated using the formula

$\kappa=\frac{|d^{2}r/dz^{2}|}{\left[1+(dr/dz)^{2}\right]^{3/2}}$ (29)

The density is obtained from the formula

$\rho=\rho_{*}\left(\frac{P}{P_{*}}\right)^{1/\gamma}=\left(\frac{\gamma-1}{\gamma+1}+\frac{2}{(\gamma+1)M_{\infty}^{2}\sin^{2}\beta_{*}^{2}}\right)^{-1}\rho_{\infty}\left(\frac{P\sin^{2}\hat{\beta}}{\hat{P}\sin^{2}\beta_{*}}\right)^{1/\gamma}$
(30)

For simplicity we introduce dimensionless units following [12]. The normalized
stream function then becomes $\Psi=\psi/\rho_{\infty}U_{\infty}L^{j+1}$, where
$L$ is a characteristic length.

On the surface of the body we need $\psi=0$, so the pressure on the body
surface $P_{b}(x)$, as well as the shock layer thickness can be obtained by
replacing the lower limits in equations (27) and (28) by zero as:

$\frac{P_{b}}{\rho_{\infty}U_{\infty}^{2}}=\frac{1}{\gamma
M_{\infty}^{2}}+(1-\hat{\chi})\sin^{2}\hat{\beta}-\frac{\hat{\kappa}}{\hat{r}^{j}}\int_{0}^{\hat{\Psi}}\left[\cos^{2}\beta_{*}+\left(\frac{2}{(\gamma-1)M_{\infty}^{2}}+\sin^{2}\beta_{*}\right)\left(1-\left(\frac{\sin^{2}\hat{\beta}}{\sin^{2}\beta_{*}}\right)^{\frac{\gamma-1}{\gamma}}\right)\right]^{1/2}d\Psi$
(31)

where the relation

$\frac{\hat{P}}{\rho_{\infty}U_{\infty}^{2}}=\frac{P_{\infty}}{\rho_{\infty}U_{\infty}^{2}}+(1-\hat{\chi})\sin^{2}\hat{\beta}=\frac{1}{\gamma
M_{\infty}^{2}}+(1-\hat{\chi})\sin^{2}\hat{\beta}$ (32)

has been used. Also, we get

$\text{for }j=0:\quad\Delta=Y\Bigr{|}_{\psi=0}$ (33)

$\text{for
}j=1:\quad\Delta=\frac{\hat{r}}{\cos\hat{\beta}}\left[1-\left(1-\frac{2Y\Bigr{|}_{\psi=0}\cos\hat{\beta}}{\hat{r}^{2}}\right)^{1/2}\right]$
(34)

with

$Y\Bigr{|}_{\psi=0}=\int_{0}^{\hat{\Psi}}\chi_{*}\left(\frac{\hat{P}\sin^{2}\beta_{*}}{P\sin^{2}\hat{\beta}}\right)^{1/\gamma}\left[\cos^{2}\beta_{*}+\left(\frac{2}{(\gamma-1)M_{\infty}^{2}}+\sin^{2}\beta_{*}\right)\left(1-\left(\frac{P\sin^{2}\hat{\beta}}{\hat{P}\sin^{2}\beta_{*}}\right)^{\frac{\gamma-1}{\gamma}}\right)\right]^{-1/2}d\Psi$
(35)

As soon as the bow shock wave is parameterized, the required quantities can be
computed by evaluating two integrals (namely equations (31) and (33) or (34)
and (35), together with the boundary values (9) and (14)-(17)). In deriving
this result, Schneider [1] emphasizes that it is not necessary to have a thin
shock layer for these equations to hold. Also, throughout the rest of the
paper, we consider only plane flow with $j=0$.

##  3 Approximation of shock shape

We now need an appropriate shock shape to use the results of Section 2. For
this, we use the approximate relation given by Moeckel in [2]. This has the
desirable property of yielding a single relation between freestream Mach
number and distance of the shock wave from the body. Detached shock waves are
normal to the freestream direction on the axis of symmetry, and are asymptotic
to the freestream Mach lines at great distances from the axis of symmetry.
Hence it seems plausible that the shape of the wave may be approximated by a
hyperbola. Moeckel therefore postulated that the equation of the shock is

$\beta r=\sqrt{z^{2}-z_{0}^{2}}$ (36)

where $\beta=\sqrt{M_{\infty}^{2}-1}$ and $z_{0}$ is the location of the
vertex of the wave. Transforming the coordinate system so that the origin is
located on the shock, we get the shock shape to be

$\beta r=\sqrt{(z+z_{0})^{2}-z_{0}^{2}}$ (37)

Here $z_{0}$ is the only unknown and it represents the shock offset distance
(which is the distance between the stagnation point on the body and the vertex
of the shock wave). Starting from a reasonable approximation for this $z_{0}$,
we can calculate the shape of the body using equation (35) (for plane flow).
By using the error between the obtained shape of the body and the actual shape
of the body (in this case, one quarter of a circle), we optimize the parameter
$z_{0}$ so that the error is minimised using secant method. After doing this,
the shock, the body obtained from (35) (for plane flow), and the actual shape
of the body for Mach number $M=4.0$ are shown in Figure 2 where the radius of
the circular cylinder is taken to be $0.5$.

Figure 2: The shock using equation (37) with $z_{0}=17.615$ (thick solid
line), the body obtained from (35) (for plane flow) (thin dashed line), and
the actual shape of the body (thin solid line) after minimising the error
between the obtained shape of the body and the actual shape of the body for
Mach number $M=4.0$. A zoomed-in view of the shock shape at two different
locations ([0.24,0.275]$\times$[0,0.05], [0.58,0.62]$\times$[0.45,0.5]) is
also shown in order to visualize the error clearly. Radius of the circular
cylinder is $0.5$

To get a better approximation for the shock shape, we use the relation given
by Moeckel as follows. We first define

$f=\frac{1}{\beta}\sqrt{(z+z_{0})^{2}-z_{0}^{2}}$ (38)

Now, we approximate the shape of the shock again by using two more parameters
$a$ and $b$ as

$r=af+bf^{2}$ (39)

We now have three unknown parameters which are $z_{0}$, $a$ and $b$. We again
determine these three parameters by optimizing them while minimising the error
between the obtained shape of the body (using equation (35)) and the actual
shape of the body. After optimizing the three parameters $z_{0}$, $a$, and
$b$, the shock, the body obtained from (35) (for plane flow), and the actual
shape of the body for Mach number $M=4.0$ are shown in Figure 3. We can now
clearly see the improvement over the previous result as the predicted body
shape is very close to the actual shape of the body.

Figure 3: The shock using equation (39) with $z_{0}=17.615$, $a=0.998$, and
$b=-0.045$ (thick solid line), the body obtained from (35) (for plane flow)
(thin dashed line), and the actual shape of the body (thin solid line) after
minimising the error between the obtained shape of the body and the actual
shape of the body for Mach number $M=4.0$. A zoomed-in view of the shock shape
at two different locations ([0.24,0.275]$\times$[0,0.05],
[0.58,0.62]$\times$[0.45,0.5]) is also shown in order to visualize the error
clearly. Radius of the circular cylinder is $0.5$

We can continue this procedure with higher degree approximations in $f$ (i.e.,
$af+bf^{2}+cf^{3}$, $af+bf^{2}+cf^{3}+df^{4}$ and so on) by optimizing each of
the coefficients of $f$ by minimising the error between the obtained shape of
the body and the actual shape of the body. We show the results obtained for
body shape using such approximations for shock shapes in Figure 4 for Mach
number $M=4.0$. We can see from Figure 4 that using higher degree
approximations (more than 2) does not really improve the body shape that much.
Even the solution obtained using such higher degree approximations is not much
better as will be explained in Section 4. For this reason, we have used
equation (39) as the shock shape for all the remaining calculations unless
otherwise specified.

Figure 4: The body obtained from (35) (for plane flow) for shock shapes given
by $f$ (dashed line), $af+bf^{2}$ (dotted line), $af+bf^{2}+cf^{3}$ (dash and
dot line), $af+bf^{2}+cf^{3}+df^{4}$ (dash and double dot line) and the actual
shape of the body (thin solid line) after minimising the error between the
obtained shape of the body and the actual shape of the body for Mach number
$M=4.0$. Radius of the circular cylinder is $0.5$

After obtaining the shock shape, we can now get the full flow field between
the shock and the body on the streamlines by evaluating equations (27), (28),
and (30). We present all our results for flow over a circular cylinder in the
next section.

##  4 Results

###  4.1 Plane flow:

We consider the plane flow past a circular cylinder for high supersonic Mach
numbers and approximate the shape of the shock using equation (39) as
discussed in Section 3. We consider five different Mach numbers
$M=4.0,5.0,6.0,7.0,$ and $8.0$ to illustrate our results. For all these Mach
numbers, the parameters $z_{0}$, $a$, and $b$ are optimized by minimising the
error between the obtained shape of the body (using equation (35)) and the
actual shape of the body. The parameters obtained for the considered Mach
numbers are tabulated in Table 1.

Mach Number | $z_{0}$ | a | b
---|---|---|---
4.0 | 17.615 | 0.998 | -0.045
5.0 | 26.755 | 0.998 | -0.050
6.0 | 38.495 | 0.998 | -0.052
7.0 | 51.982 | 0.998 | -0.054
8.0 | 67.984 | 0.998 | -0.058
Table 1: The optimized parameters $z_{0}$, $a$, and $b$ in (39) obtained by
minimising the error between the obtained shape of the body (using equation
(35)) and the actual shape of the body.

After using this shock shape, we now calculate the full flow field between the
shock and the body on the streamlines by evaluating equations (27), (28), and
(30). We have used 200 streamlines for each Mach number. The density
variations for Mach numbers $4.0$, $5.0$, $6.0$, $7.0$, and $8.0$, obtained
with this shock shape are shown in Figures 5, 6, 7, 8, and 9 respectively. All
the plots are obtained by plotting the point data and using a delaunay
triangulation to convert them into cells to cover the whole domain using
paraview [16].

Figure 5: Density variation for Mach number $4.0$ obtained using (27), (28),
and (30) with shock shape given in (39) with parameters given in Table 1 for
200 streamlines. Figure 6: Density variation for Mach number $5.0$ obtained
using (27), (28), and (30) with shock shape given in (39) with parameters
given in Table 1 for 200 streamlines. Figure 7: Density variation for Mach
number $6.0$ obtained using (27), (28), and (30) with shock shape given in
(39) with parameters given in Table 1 for 200 streamlines. Figure 8: Density
variation for Mach number $7.0$ obtained using (27), (28), and (30) with shock
shape given in (39) with parameters given in Table 1 for 200 streamlines.
Figure 9: Density variation for Mach number $8.0$ obtained using (27), (28),
and (30) with shock shape given in (39) with parameters given in Table 1 for
200 streamlines.

We now validate our procedure by using a numerically calculated solution. We
have used a fifth order accurate (formally called $\mathbf{P}^{4}$ based)
discontinuous Galerkin method (DGM) along with overset grids with 9976
elements for accurate shock capturing [15] to solve the two-dimensional Euler
equations over a circular cylinder and used this solution to calculate the
density error in our solution. The formulation and validation of our solver
using DGM, along with a brief outline of the shock capturing procedure are
given in Appendix A. The density errors obtained using this fifth order
accurate numerical solution for Mach number $4.0$, $5.0$, $6.0$, $7.0$, and
$8.0$, are shown in Figures 5, 11, 12, 13, and 14 respectively. We don’t show
the solution obtained by the discontinuous Galerkin method as it is quite
similar to the solution obtained using our procedure and shown in Figures 5 \-
9.

Figure 10: Density error obtained for Mach number $4.0$ using a fifth order
accurate (formally called $\mathbf{P}^{4}$ based) discontinuous Galerkin
method along with overset grids using 9976 elements for accurate shock
capturing [15]. Figure 11: Density error obtained for Mach number $5.0$ using
a fifth order accurate (formally called $\mathbf{P}^{4}$ based) discontinuous
Galerkin method along with overset grids using 9976 elements for accurate
shock capturing [15]. Figure 12: Density error obtained for Mach number $6.0$
using a fifth order accurate (formally called $\mathbf{P}^{4}$ based)
discontinuous Galerkin method along with overset grids using 9976 elements for
accurate shock capturing [15]. Figure 13: Density error obtained for Mach
number $7.0$ using a fifth order accurate (formally called $\mathbf{P}^{4}$
based) discontinuous Galerkin method along with overset grids using 9976
elements for accurate shock capturing [15]. Figure 14: Density error obtained
for Mach number $8.0$ using a fifth order accurate (formally called
$\mathbf{P}^{4}$ based) discontinuous Galerkin method along with overset grids
using 9976 elements for accurate shock capturing [15].

We now explain the choice of the selected Mach numbers. We have calculated the
maximum density error obtained using the procedure in Section 2 and tabulated
them in Table 2. We also show the value of $\rho_{\infty}/\rho$ across a
normal shock in Table 2 for reference. We can clearly see that the error is
increasing quite fast when we go below Mach $4.0$. In fact, it has increased
by two orders of magnitude when we go from Mach $4.0$ to Mach $3.5$. The
reason for this can be the first assumption for this procedure given by
equation (18) where we have assumed that the density immediately behind a
strong shock is much larger than in front of the shock. The values of
$\rho_{\infty}/\rho$ (across a normal shock) given in Table 2 decrease when we
increase the Mach number and we probably have to cross a threshold for this
procedure to work well which can be between Mach numbers 3.5 and 4.0 as seen
from Table 2. For this reason, we have selected Mach numbers $\geq 4$ to
illustrate our solution.

Mach Number | Maximum Density Error |  $\rho_{\infty}/\rho$ across a normal shock
---|---|---
3.5 | 1.1E-01 | 0.2347
3.8 | 8.9E-02 | 0.2244
3.9 | 2.1E-02 | 0.2214
4.0 | 1.9E-03 | 0.2188
5.0 | 2.0E-03 | 0.2
6.0 | 3.8E-03 | 0.1898
7.0 | 5.1E-03 | 0.1837
8.0 | 5.9E-03 | 0.1797
Table 2: Maximum density error obtained using the procedure given in Section 2
for various Mach numbers. For these Mach numbers, the density ratio
($\rho_{\infty}/\rho$) across a normal shock is also shown. The error is
calculated using a fifth order accurate discontinuous Galerkin method along
with overset grids using 9976 elements for accurate shock capturing [15]

To further validate the solution obtained, we calculate the functional -
integrated surface pressure on the body of the cylinder and compare it with
the numerical solution obtained. The error in integrated surface pressure is
tabulated in Table 3. From the values obtained, we can see that the error in
the integrated surface pressure is smaller than the maximum density error.
This shows that the procedure is quite accurate near the body (except at some
locations as will be discussed later) in predicting the surface pressure which
is an important quantity.

Mach Number | Error in integrated surface pressure
---|---
4.0 | 3.45E-04
5.0 | 2.75E-04
6.0 | 2.15E-04
7.0 | 2.03E-04
8.0 | 1.88E-04
Table 3: Error in the functional - integrated surface pressure for various
Mach numbers using the shock shape given by $af+bf^{2}$. The error is
calculated using a fifth order accurate discontinuous Galerkin method along
with overset grids using 9976 elements for accurate shock capturing [15]

Discussion of results: We can see that the results obtained for each Mach
number are very accurate and the maximum error in comparison with a numerical
solution for each case is only about the order of $10^{-3}$ (using a shock
shape given by $af+bf^{2}$). For all the Mach numbers, the error pattern is
similar and the maximum error occurs near the upper portion of the body. To
explain this, we show the pressure solution for Mach $4.0$ in Figure 15. We
also remember the assumption that the pressure at the point Q (Figure 1) of
the disturbed flow field is not much smaller than the pressure at point N (see
Figure 1). This assumption is given by equation (19). From Figure 15, we can
see that for the points that are near the upper portion of the body, the
pressure ratio $\hat{P}/P$ is larger compared to the other points. Though we
showed the pressure variation for Mach 4.0 only, this happens to be true for
Mach numbers 5.0, 6.0, 7.0, and 8.0 as well. This can be why the error
obtained is the largest in that region as we are slightly violating an
assumption with which we derived the solution.

Figure 15: Pressure variation for Mach number $4.0$ with shock shape given in
(39) with parameters given in Table 1 for 200 streamlines.

We have also changed the shock shape further by looking at equations of the
form $af+bf^{2}+cf^{3},af+bf^{2}+cf^{3}+df^{4}$ with $f$ defined by (38). We
have optimized the parameters $z_{0}$, $a$, $b$, $c$, and $d$ as mentioned in
Section 2. We tabulate the maximum density errors for shock shapes with
polynomial expansions in $f$ till fourth degree in Table 3. From this table,
we can see that the solution accuracy increases only till the second degree
approximation in $f$ and doesn’t improve much for higher degree
approximations.

As a final shock shape, we consider the accurate numerical solution obtained
in [15] where the shock has been captured aligned to a grid line. We use that
solution and obtain a cubic spline for the shock shape. Using the cubic spline
as the shock shape, we obtain the body and the full solution in the shock
layer using the method given in Section 2. The shape of the body obtained
along with the shock and the exact shape of the body are shown in Figure 16
for Mach $4.0$ flow. We can see that the shape of the body obtained in this
fashion is quite accurate. The density error for this solution is shown in
Figure 17. We can see that the maximum density error is about $1.4\times
10^{-3}$. We also tabulate the maximum density error obtained in this fashion
for various Mach numbers in Table 4.

Figure 16: The shock using the cubic spline from the numerical solution (thick solid line), the body obtained from (35) (for plane flow) (thin dashed line), and the actual shape of the body (thin solid line) for Mach number $M=4.0$. Radius of the circular cylinder is $0.5$ Figure 17: Density error obtained for Mach number $4.0$ using the cubic spline shock shape from the numerical solution using a fifth order accurate (formally called $\mathbf{P}^{4}$ based) discontinuous Galerkin method along with overset grids using 9976 elements for accurate shock capturing [15]. | Maximum density error with shock shape using Schneider’s method for plane flow:
---|---
Mach number | $f$ |  $af+bf^{2}$ (shock shape used for all our calculations) | $af+bf^{2}+cf^{3}$ | $af+bf^{2}+cf^{3}+df^{4}$ |  cubic spline obtained from numerical solution
4.0 | 1.8E-02 | 1.9E-03 | 1.8E-03 | 1.8E-03 | 1.4E-03
5.0 | 1.9E-02 | 2.0E-03 | 1.7E-03 | 1.7E-03 | 1.1E-03
6.0 | 2.1E-02 | 3.8E-03 | 3.5E-03 | 3.4E-03 | 2.4E-03
7.0 | 2.2E-02 | 5.1E-03 | 4.9E-03 | 4.8E-03 | 3.6E-03
8.0 | 2.8E-02 | 5.9E-03 | 5.8E-03 | 5.5E-03 | 4.3E-03
Table 4: Plane flow: Maximum density error obtained for various Mach numbers
with various shock shapes obtained using an expansion in terms of $f$ and also
using a cubic spline obtained from the numerical solution. The error is
calculated using a fifth order accurate discontinuous Galerkin method along
with overset grids using 9976 elements for accurate shock capturing [15]

From these results, we can conclude that using this kind of a shock shape,
this might be the best accuracy we can get using this analytical method with
the stated approximations and assumptions.

To see how well the shock shapes assumed work in a typical shock fitting
solver, we use the shock shapes $f$, $af+bf^{2}$, $af+bf^{2}+cf^{3}$,
$af+bf^{2}+cf^{3}+df^{4}$ and optimise the various parameters ($z_{0},a,b,c,d$
as required) based on the boundary condition at the circular cylinder (the no-
penetration condition). We get different values of the optimised parameters to
those obtained using Schneider’s method (as expected) and we tabulate the
maximum density error for each of the shock shapes in Table 5. Here unlike the
case for Schneider’s method, we see that the maximum density error decreases
reasonably as we expand the shock shape in series expansion. The decrease is
quite prominent (two orders of magnitude) when we go from $f$ to $af+bf^{2}$.
The final column of the table is given just to show that the shock shape
obtained from numerical solution as a cubic spline solves the problem very
accurately. All these results show that the choice of shock shape used for our
solution is quite good.

| Maximum density error with shock shape with shock fitting using DGM for
plane flow:
---|---
Mach number | $f$ | $af+bf^{2}$ | $af+bf^{2}+cf^{3}$ | $af+bf^{2}+cf^{3}+df^{4}$ |  cubic spline obtained from numerical solution
4.0 | 3.3E-06 | 5.5E-08 | 3.7E-09 | 4.1E-10 | 1.8E-16
5.0 | 2.4E-06 | 2.3E-08 | 4.4E-09 | 3.6E-10 | 2.1E-16
6.0 | 4.1E-06 | 4.2E-08 | 3.2E-09 | 3.9E-10 | 1.3E-16
7.0 | 5.2E-06 | 6.3E-08 | 8.3E-09 | 6.7E-10 | 2.2E-16
8.0 | 4.8E-06 | 7.7E-08 | 6.9E-09 | 7.8E-10 | 3.2E-16
Table 5: Plane flow: Maximum density error obtained for various Mach numbers
using shock fitting with DGM with various shock shapes obtained using an
expansion in terms of $f$ and also using a cubic spline obtained from the
numerical solution. The error is calculated using a fifth order accurate
discontinuous Galerkin method along with overset grids using 9976 elements for
accurate shock capturing [15]

###  4.2 Axisymmetric Flow:

We have also tested Schneider’s method given in Section 2 for the axisymmetric
case. To obtain an accurate numerical solution for this case, we have
converted our plane flow solver into an axisymmetric flow solver using the
method given by Yu [17]. This gives us a fifth order accurate discontinuous
Galerkin solution for axisymmetric flow. We have used this numerical solution
with 7500 elements to calculate the errors in the analytical solution. We used
the same shock shape given by Moeckel with

$r=f$ (40)

where $f$ is given by equation (38). Here, $r$ is the radial coordinate and
$z$ is the longitudinal coordinate.

We now use the shock shapes $f$, $af+bf^{2}$, $af+bf^{2}+cf^{3}$,
$af+bf^{2}+cf^{3}+df^{4}$ and optimise the various parameters ($z_{0},a,b,c,d$
as required) again using the error between the exact body shape and the body
shape obtained using the method given in Section 2. We get different values of
the optimised parameters to those obtained for plane flow (as expected) and we
tabulate the maximum density error for each of the shock shapes in Table 6.
Here, we can see that the pattern for maximum density error is similar to that
of plane flow decreasing from shock shapes $f$ to $af+bf^{2}$ and remaining
almost the same for further series expansions. We also show the error obtained
using a shock shock obtained from numerical solution in the last column of
Table 6. From these results, we can again conclude that using this kind of a
shock shape, this might be the best accuracy we can get for axisymmetric flow
using the analytical method given in Section 2.

| Maximum density error with shock shape using Schneider’s method for
axisymmetric flow:
---|---
Mach number | $f$ | $af+bf^{2}$ | $af+bf^{2}+cf^{3}$ | $af+bf^{2}+cf^{3}+df^{4}$ |  cubic spline obtained from numerical solution
4.0 | 8.5E-02 | 4.3E-03 | 4.2E-03 | 3.8E-03 | 1.6E-03
5.0 | 8.3E-02 | 3.4E-03 | 3.1E-03 | 2.8E-03 | 1.3E-03
6.0 | 7.4E-02 | 3.2E-03 | 2.8E-03 | 2.6E-03 | 1.1E-03
7.0 | 7.3E-02 | 2.8E-03 | 2.5E-03 | 2.3E-03 | 1.2E-03
8.0 | 6.6E-02 | 2.6E-03 | 2.3E-03 | 2.2E-03 | 1.7E-03
Table 6: Axisymmetric Flow: Maximum density error obtained for various Mach
numbers with various shock shapes obtained using an expansion in terms of $f$
and also using a cubic spline obtained from the numerical solution. The error
is calculated using a fifth order accurate discontinuous Galerkin solver for
axisymmetric flow with 7500 elements.

##  5 Conclusion

We have calculated the high supersonic flow over a circular cylinder using
Schneider’s inverse method [1]. We assumed the shock shape of a hyperbola in
the form given by Moeckel in [2] and optimized the shock shape by minimising
the error between the shape of the body obtained using Schneider’s method and
the actual shape of the body. We have improved this shock shape further by
writing it as a series expansion ($af+bf^{2}$, where $f$ gives the shock shape
of Moeckel) with unknown coefficients. We have determined these coefficients
using the same optimization and obtained a much better approximation for the
shock shape. We have used this shock shape to find the solution of the flow in
the shock layer using Schneider’s method by integrating the equations of
motion using the stream function. We have compared the solution obtained with
a numerical solution calculated using a fifth order accurate discontinuous
Galerkin method which captures the shock very accurately using overset grids
[15]. We have found that the maximum error in density is only of the order of
$10^{-3}$ (with the shock shape $af+bf^{2}$) demonstrating the accuracy of the
solution method. We also found that using more terms in the series expansion
for the shock shape (like $af+bf^{2}+cf^{3}$, $af+bf^{2}+cf^{3}+df^{4}$) does
not improve the accuracy of the solution. We also used a cubic spline obtained
from the numerical solution to find the solution and the maximum error in
density is again of the order of $10^{-3}$. This suggests that this might be
the best accuracy we can get using this analytical method. We have also used
the shock shape given by Moeckel and used it in a series expansion for
axisymmetric flow case and optimised the parameters. The density errors again
follow the same pattern of plane flow and the maximum density error is again
of the order of $10^{-3}$.

##  Appendix A Formulation and Validation of discontinuous Galerkin Method

Consider the Euler equations in conservative form as given by

$\frac{\partial\textbf{Q}}{\partial t}+\frac{\partial\textbf{F(Q)}}{\partial
x}+\frac{\partial\textbf{G(Q)}}{\partial y}=0\quad\text{in the
domain}\quad\Omega$ (41)

where $\textbf{Q}=(\rho,\rho u,\rho v,E)^{T}$,
$\textbf{F(Q)}=u\textbf{Q}+(0,P,0,Pu)^{T}$ and
$\textbf{G(Q)}=v\textbf{Q}+(0,0,P,Pv)^{T}$ with
$P=(\gamma-1)(E-\frac{1}{2}\rho(u^{2}+v^{2}))$ and $\gamma=1.4$. Here, $\rho$
is the density, $(u,v)$ is the velocity, $E$ is the total energy and $P$ is
the pressure. We approximate the domain $\Omega$ by $K$ non overlapping
elements given by $\Omega_{k}$.

We look at solving (41) using the discontinuous Galerkin method (DGM). We
approximate the local solution in an element $\Omega_{k}$, where $k$ is the
element number, as a polynomial of order $N$ which is given by:

$Q_{h}^{k}(r,s)=\sum_{i=0}^{N_{p}-1}Q_{i}^{k}\psi_{i}(r,s)$ (42)

where $N_{p}=(N+1)(N+1)$ and $r$ and $s$ are the local coordinates. Here, the
subscript $i$ represents the particular degree of freedom, $h$ represents the
grid size, and the superscript $k$ is the element number. The polynomial basis
used ($\psi_{i}(r,s)$) is the tensor product orthonormalized Legendre
polynomials of degree $N$. The number of degrees of freedom are given by
$N_{p}=(N+1)(N+1)$. Now, using $\psi_{j}(r,s)$ as the test function, the weak
form of the equation (41) is obtained as

$\sum_{i=0}^{N_{p}-1}\frac{\partial Q_{i}^{k}}{\partial
t}\int_{\Omega_{k}}\psi_{i}\psi_{j}d\Omega+\int_{\partial\Omega_{k}}\hat{F}\psi_{j}ds-\int_{\Omega_{k}}\vec{F}.\nabla\psi_{j}d\Omega=0\quad
j=0,\ldots,N_{p}-1$ (43)

where $\partial\Omega_{k}$ is the boundary of $\Omega_{k}$,
$\vec{F}=(\textbf{F(Q)},\textbf{G(Q)})$ and $\hat{F}=\bar{F^{*}}.\hat{n}$
where $\bar{F^{*}}$ is the monotone numerical flux at the interface which is
calculated using an exact or approximate Riemann solver and $\hat{n}$ is the
unit outward normal. We have used the Lax-Friedrichs numerical flux for all
our calculations unless otherwise specified. This is termed to be
$\mathbf{P}^{N}$ based discontinuous Galerkin method.

Equation (43) is integrated using an appropriate Gauss Legendre quadrature and
is discretized in time by using the fifth order Runge-Kutta time
discretization given in [18] unless otherwise specified. To control spurious
oscillations which occur near discontinuities, a limiter is used with a
troubled cell indicator (a shock detector). We have used the KXRCF troubled
cell indicator [19] and the compact subcell WENO (CSWENO) limiter proposed in
[20] for all our calculations. To capture the shock accurately, we have used
overset grids as proposed in [15]. The communication between overset grids
occurs using the procedure given in [21]. The shock capturing procedure with
overset grids is briefly given below.

Shock capturing procedure with overset grids: The step by step procedure is as
follows:

Step 1: Run the solver on a coarse grid with a given troubled cell indicator
and limiter to steady state and obtain the solution.

Step 2: Look at the troubled cells (with a reliable shock detector [19]) to
locate the discontinuities (shocks) that occur in the solution.The troubled
cells give us a good idea of the location of the discontinuities.

Step 3: Construct an overset grid conforming to the computational domain which
is refined in a direction perpendicular to the discontinuities such that they
are approximately parallel to a grid line.

Step 4: Using this overset grid, we rerun the solver with the coarse grid
solution as the initial condition. While running the solver, we use the
troubled cell indicator and the limiter only on the overset grid. We also use
a high resolution numerical flux on the overset grid to capture the shock
accurately. We have used the SLAU2 [22] numerical flux in the overset grid and
the less expensive Lax-Friedrichs flux elsewhere. Using this procedure, we
obtain a more accurate solution with the discontinuities approximately aligned
to a grid line.

Validation of the solver using Isentropic Vortex Problem: Consider the two-
dimensional Euler equations given by equation (41) in the domain
$[0,10]\times[-5,5]$ for the Isentropic Euler Vortex problem suggested in [23]
as a test case. The analytical solution is given by:
$\rho=\left(1-\left(\frac{\gamma-1}{16\gamma\pi^{2}}\right)\beta^{2}e^{2(1-r^{2})}\right)^{\frac{1}{\gamma-1}}$,
$u=1-\beta e^{(1-r^{2})}\frac{y-y_{0}}{2\pi}$, $v=\beta
e^{(1-r^{2})}\frac{x-x_{0}-t}{2\pi}$, and $p=\rho^{\gamma}$, where
$r=\sqrt{(x-x_{0}-t)^{2}+(y-y_{0})^{2}}$, $x_{0}=5$, $y_{0}=0$, $\beta=5$ and
$\gamma=1.4$. Initialization is done with the exact solution at $t=0$ and
periodic boundary conditions are used at the edges of the domain. The solution
is obtained for various orders making sure that the spatial (DGM) and temporal
(Runge-Kutta) schemes are of the same order. The errors in density and
numerical orders of accuracy are calculated at $t=2$ for various grid sizes
and are presented in Table 7. We can clearly see that the theoretical order of
accuracy of the scheme is obtained by the DG solver.

| $K\times K$ |  $L_{2}$ error | Order |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order
---|---|---|---|---|---|---|---
$P^{1}$ | 20$\times$20 | 1.58E-03 |  | 2.41E-03 |  | 1.77E-01 |
40$\times$40 | 3.66E-04 | 2.11 | 5.54E-04 | 2.12 | 4.39E-02 | 2.01
80$\times$80 | 7.28E-05 | 2.33 | 1.11E-04 | 2.32 | 1.11E-02 | 1.99
160$\times$160 | 1.39E-05 | 2.39 | 2.15E-05 | 2.37 | 2.80E-03 | 1.98
$P^{2}$ | 20$\times$20 | 5.19E-04 |  | 6.33E-04 |  | 9.42E-02 |
40$\times$40 | 6.27E-05 | 3.05 | 7.70E-05 | 3.04 | 1.15E-02 | 3.03
80$\times$80 | 7.11E-06 | 3.14 | 8.67E-06 | 3.15 | 1.30E-03 | 3.15
160$\times$160 | 8.35E-07 | 3.09 | 1.01E-06 | 3.10 | 1.62E-04 | 3.00
$P^{3}$ | 20$\times$20 | 3.30E-06 |  | 4.66E-06 |  | 8.21E-05 |
40$\times$40 | 2.06E-07 | 4.00 | 2.93E-07 | 3.99 | 5.17E-06 | 3.99
80$\times$80 | 1.27E-08 | 4.02 | 1.85E-08 | 3.99 | 3.30E-07 | 3.97
160$\times$160 | 7.89E-10 | 4.01 | 1.12E-09 | 4.04 | 2.09E-08 | 3.98
$P^{4}$ | 20$\times$20 | 1.76E-07 |  | 3.16E-07 |  | 9.83E-06 |
40$\times$40 | 5.54E-09 | 4.99 | 1.01E-08 | 4.97 | 3.32E-07 | 4.89
80$\times$80 | 1.78E-10 | 4.96 | 3.24E-10 | 4.96 | 1.10E-08 | 4.92
160$\times$160 | 5.68E-12 | 4.97 | 1.05E-11 | 4.95 | 3.57E-10 | 4.94

Table 7: 2D Euler equations for the Isentropic Vortex problem with periodic
boundary conditions, $t=2$, uniform mesh with $K\times K$ elements, $L_{1}$,
$L_{2}$ and $L_{\infty}$ errors for $P^{1}$, $P^{2}$, $P^{3}$ and $P^{4}$
based DG

## References

  * [1] W. Schneider, “A uniformly valid solution for the hypersonic flow past blunted bodies.,” Journal of Fluid Mechanics, vol. 31, pp. 397–415, 1968\.
  * [2] W. Moeckel, “Approximate Method for Predicting Form and Location of Detached Shock Waves Ahead of Plane or Axially Symmetric Bodies.,” NACA Technical Note, No.1921, 1949.
  * [3] W. Chester, “Supersonic flow past a bluff body with a detached shock. Part I, Two-dimensional body.,” Journal of Fluid Mechanics, vol. 1, pp. 353–365, 1956.
  * [4] W. Chester, “Supersonic flow past a bluff body with a detached shock. Part II, Axisymmetrical body.,” Journal of Fluid Mechanics, vol. 1, pp. 490–496, 1956.
  * [5] N. Freeman, “On the theory of hypersonic flow past plane and axially symmetric bluff bodies.,” Journal of Fluid Mechanics, vol. 1, pp. 366–387, 1956.
  * [6] W. Hayes and R. Probstein, Hypersonic flow theory.  Academic Press, New York, 1966.
  * [7] M. Rasmussen, Hypersonic flow.  Wiley, 1994.
  * [8] M. D. Salas, A Shock-Fitting Primer.  CRC Applied Methematics and Nonlinear Science, Chapman and Hall, 1st edition, 2009.
  * [9] Y. Chou, “Radiative Coupled Viscous Flow with Massive Blowing.,” NASA CR-2236, 1973.
  * [10] S. Maslen, “Inviscid Flow about Blunted Cones of Large Opening Angle at Angle of Attack.,” NASA CR-132652, 1975.
  * [11] K. Nagaraja, “On the Shock Stand-off Distance in an Inviscid Hypersonic Source Flow Past Two-Dimensional Bluff Bodies.,” ARL 71-0303, 1971.
  * [12] M. Schulreich and D. Breitschwerdt, “Astrophysical bow shocks:an analytical solution for the hypersonic blunt body problem in the intergalactic medium.,” Astronomy and Astrophysics, vol. 531, A13, 2011.
  * [13] K. Scherer, F. Fichtner, J. Kleimann, T. Wiengarten, D. Bomans, and K. Weis, “Shock structures of astrospheres.,” Astronomy and Astrophysics, vol. 586, A111, 2016.
  * [14] H. Schwarze, “A Uniformly Valid Solution for the Three-Dimensional Inviscid Supersonic Flow past Blunt Bodies with Strong Compression in the Shock Wave.,” Acta Mechanica, vol. 42, pp. 11–35, 1982.
  * [15] S. R. Siva Prasad Kochi and M. Ramakrishna, “Shock capturing with discontinuous Galerkin Method using Overset grids for two-dimensional Euler equations,” accepted for publication in Journal of Computational Physics.
  * [16] J. Ahrens, B. Geverci, and C. Law, ParaView: An End-User Tool for Large Data Visualization.  Elsevier, 2005.
  * [17] S.-T. Yu, “Convenient Method to Convert Two-Dimensional CFD Codes into Axisymmetric Ones.,” Journal of Propulsion and Power, vol. 9, pp. 493–495, 1993.
  * [18] J. C. Butcher, Numerical Methods for Ordinary Differential Equations.  John Wiley and Sons, 2016.
  * [19] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, and J. Flaherty, “Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws.,” Appl. Numer. Math., vol. 48, pp. 323–338, 2004.
  * [20] S. R. Siva Prasad Kochi and M. Ramakrishna, “A compact subcell WENO limiting strategy using immediate neighbours for Runge-Kutta discontinuous Galerkin methods.,” International Journal of Computer Mathematics, vol. 98, no. 3, pp. 608–626, 2021.
  * [21] S. R. Siva Prasad Kochi and M. Ramakrishna, “A Discontinuous Galerkin Overset Scheme Using WENO Reconstruction and Subcells for Two-Dimensional Problems,” Journal of Scientific Computing, vol. 93:35, 2022.
  * [22] K. Kitamura and E. Shima, “Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes.,” Journal of Computational Physics, vol. 245, pp. 62–83, 2013.
  * [23] C.-W. Shu, “Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.,” Lecture Notes in Mathematics, Springer, vol. 1697, pp. 325–432, 1998.