# Convolutions on Spherical Images

Marc Eder Jan-Michael Frahm
University of North Carolina at Chapel Hill
{meder<EMAIL_ADDRESS>
###### Abstract

Applying convolutional neural networks to spherical images requires particular
considerations. We look to the millennia of work on cartographic map
projections to provide the tools to define an optimal representation of
spherical images for the convolution operation. We propose a representation
for deep spherical image inference based on the icosahedral Snyder equal-area
(ISEA) projection, a projection onto a geodesic grid, and show that it vastly
exceeds the state-of-the-art for convolution on spherical images, improving
semantic segmentation results by $12.6\%$.

##  1 Introduction and Related Work

Omnidirectional imaging is becoming increasingly popular thanks to the
proliferation of consumer-grade $360^{\circ}$ cameras and the benefit of the
wide field-of-view to a number of applications. At the same time, deep
learning using convolutional neural networks (CNNs) has never been a more
widely-used tool for computer vision tasks. Yet, applying CNNs to spherical
images is not quite as trivial as simply training a vanilla CNN on spherical
projections onto an image plane. Any planar representation of a sphere
necessarily contains some degree of content deformation. This violates
convolution’s translation equivariance, which requires the convolution kernel
and the signal to be uniformly discretized. To address the unique problems of
convolution on spherical images, several methods have been proposed.

Figure 1: We show that convolving on a planar geodesic approximation to a
sphere drastically improves results for the dense prediction task of semantic
segmentation.

Su and Grauman [13] train a CNN to transfer center-perspective-projection-
trained models to the equirectangular domain using a learnable adaptive
convolutional kernel to match the outputs. Observation that certain regions of
cube-map faces contain less distortion than others, Xiong and Grauman [16]
develop a content-aware rotation of the spherical data that “snaps-to” a view
where the relevant information has minimal distortion. Zioulis _et al_. [17]
explore the use of rectangular filter-banks, rather than square filters, to
horizontally expand the network’s receptive field to account for distortion in
equirectangular images. The spherical convolution derived by Cohen _et al_.
[3, 4] uses a generalized formulation of convolution to filter rotations of
the feature maps rather than translations. Finally, work from Coors _et al_.
[5] and Tateno _et al_. [14] develop a location-dependent transformation of
the planar convolutional kernel using a gnomonic (rectilinear) projection to
adapt to the spatial distortion in spherical images.

We take a slightly different approach than these previous works. Rather than
changing the tools of convolution, we look to find a better representation of
the data. We first analyze the convolution operation to identify the important
properties that must be considered when applying convolutions to spherical
images and use this to determine a more optimal way to represent the data. We
resolve that a planar geodesic approximation to the sphere is a more optimal
representation and show that its use yields state-of-the-art performance for
dense semantic segmentation.

##  2 Cartographic Map Projections

First, we briefly review cartographic map projection properties and explain
their relevance to this problem.

(a) Gnomonic

(b) Equirectangular

(c) Mercator

(d) Gall-Peters

(e) ISEA (Icosphere)

Figure 2: Different map projections with superimposed Tissot’s indicatrices to
illustrate the distortion. Ellipse eccentricity shows conformality, size shows
areal equality, and spacing shows equidistance. Note that the Gnomonic
projection displayed only depicts a $90^{\circ}\times 90^{\circ}$ segment of
the Earth, as would a face on a cube-map. Earth image from [1].

###  2.1 Distortion in map projections

Gauss proved with his Theorema Egregium that it is impossible to represent a
sphere on a plane without some degree of distortion. The nature of this
distortion is dependent on the projection model used to map the sphere to the
plane. There are numerous cartographic projections that satisfy useful
properties for specific applications, like navigation, but generally, there
are three major projection properties to consider:

  * •

Equal area: Preserves relative scale of objects

  * •

Conformal: Preserves local angles

  * •

Equidistant: Preserves distance between locations on the sphere

Resulting from Gauss’s theorem, it is impossible for a projection to be both
conformal and equal area, and similarly a map cannot be equidistant everywhere
[10]. However, certain projections may provide one or another of these
properties. A typical method for visualizing the distortion properties of
different map projections is to use Tissot’s indicatrices [10], which
illustrate the local distortions at points on the map. A Tissot’s indicatrix
represents the projection onto a map of a circle with infinitesimal radius on
the surface of the sphere. The result is an ellipse where the major and minor
axes relate the scale and angular distortions. By placing the circles with
regular spacing on the sphere, the indicatrices also depict the distortion of
distances in the projection.

Figures 2(a)-2(d) show Earth represented by the gnomonic (rectilinear)
projection used for faces of cube-maps, the equirectangular (plate carée)
projection, the common Mercator projection, and the Gall-Peters projection,
respectively, with Tissot’s indicatrices overlaid. The gnomonic projection is
neither conformal, equal areal, nor equidistant, demonstrated by the varied
eccentricity and shape of the ellipses and the spacing between them. The
equirectangular projection is an equidistant projection, preserving the
distance between parallels of latitude, but is neither equal-area nor
conformal, so the indicatrices are neither circular nor equal in size. The
Mercator projection is conformal but not equal-area, and thus relative scale
varies throughout the map. Note how Antarctica appears significantly larger
than Africa, despite actually having less than 50% of Africa’s area. All
indicatrices remain circular as it is conformal, but they vary in size
depending on location. Finally, the Gall-Peters projection is an equal-area
projection, which addresses this scale imbalance by preserving the relative
size of objects in the map. This method is not conformal, though, and
therefore local angles are not preserved, which can be observed by the varying
eccentricities of the ellipses.

Figure 3: Simple encoder-decoder network architecture used for our
experiments. Each layer but the last is followed by an exponential linear unit
activation function [2] without batch normalization. The layers are described
as (input_channel, output_channel, kernel_size, up/downsampling) and ‘same’
padding is used for all layers.

###  2.2 Relevance to convolution

The 1D discrete convolution of a filter $f$ of size $K$ and a signal $g$ is
given in Eq. (1):

$(f*g)[n]=\sum_{m=-\lfloor\frac{K}{2}\rfloor}^{\lfloor\frac{K}{2}\rfloor}f[m]g[n-m].$
(1)

The operation has two components: a sampling and a weighted summation, or more
explicitly:

$(f*g)[n]=\sum_{m=-\lfloor\frac{K}{2}\rfloor}^{\lfloor\frac{K}{2}\rfloor}f[m]\left(\sum_{l=-\infty}^{\infty}g[l]\delta[l-n+m]\right)$
(2)

where $\delta[\cdot]$ is the Dirac delta function. The Dirac delta function
can be expressed as the limit of a zero-centered, isotropic Gaussian
distribution as the variance goes to zero, and thus in 2D can be represented
by an infinitely small circle. Hence when convolving with a spherical image,
each sample is an infinitesimal circle on the surface of a sphere. As
discussed previously, we can examine how this sampling is affected by
distortion using Tissot’s indicatrix.

From the definition of convolution, the kernel must sample the same area at
each location in the data. Thus, there is an implicit assumption that the data
is undistorted (i.e. the indicatrices are all perfect circles of equal area).
As illustrated by Figure 2(b), this is not the case with equirectangular
projections. Neither does it hold true for cube-maps, which use the gnomonic
projection shown in Figure 2(a). These two common spherical image formats
violate this distortion-free assumption, which explains why the performance of
traditional CNNs on these images falls short when compared to center
perspective images that typically have little to no distortion. It follows
that there are two viable solutions for resolving this issue: replace the
Direct delta function in Eq. (2) with a bivariate, anisotropic Gaussian with
location-dependent variance, or resample the image to a different projection
with less distortion. The former approach would slow down convolution on
spherical images, as a large variance at a sample location could require the
accumulation of large regions of pixels, instead of just sampling at one. The
latter approach is thus preferable. We still are blocked by the theoretical
limitations of map projections, but we can select a compromise projection with
the least distortion.

###  2.3 Icosahedral grids

To re-state, we desire a map projection that is: 1) nearly conformal,
preserving the local shape of the data; 2) nearly equal-area, maintaining the
local size of the data; and 3) nearly equidistant, allowing fixed-size
convolutional filters to retain translational equivariance over the entire
image. This suggests that the best projection should yield a very close planar
approximation to the sphere itself. Stemming from R. Buckminster Fuller’s work
on discretizing the sphere [6], researchers have developed various geodesic
approximations for cartographic map projections. One such projection, the
icosahedral Snyder equal area (ISEA) projection [11], is particularly well-
suited for our task of representing spherical images. It projects the data on
the sphere’s surface onto the planar faces of a recursively subdivided
icosahedron using a triangle partitioning. This compromise projection is among
the least distorted among geodesic grids [7], which satisfies our
aforementioned criteria, and it is easily represented by a triangular 3D mesh.
Therefore, we propose that this ISEA projection is an ideal spherical image
representation for inference with CNNs.

To perform the ISEA projection, we first create a $7^{th}$ order “icosphere,”
a regular icosahedron subdivided $7$ times using Loop subdivision [9] with the
vertices subsequently normalized to the unit sphere. We choose a $7^{th}$
order icosphere because the number of vertices ($163,842$) is most similar to
the number of pixels ($131,072$) in the $256\times 512$ equirectangular images
we use for comparison. Next, we resample the image to the vertices of the
icosphere using barycentric interpolation on the faces. The result of this
projection is shown in Figure 2(e), along with the Tissot indicatrices of
deformation on the geodesic’s surface. Notice how the indicatrices are all
nearly equally-sized circles. This illustrates that the content deformation is
very limited, and the convolution’s sampling, as described in Section 2.2,
will be nearly unaffected. We use this representation for all spherical image
inputs in our experiments.

Representation | Floor | Ceiling | Wall | Door | Cabinet | Rug | Window | Curtain
---|---|---|---|---|---|---|---|---
Equirect. (gnom. kernel) [5, 14] | $0.9315$ | $\mathbf{0.9710}$ | $0.8597$ | $0.6466$ | $0.6376$ | $\mathbf{0.7284}$ | $0.7012$ | $0.4703$
Equirect. (equirect. kernel) | $0.9327$ | $0.9654$ | $0.8566$ | $0.6692$ | $0.6638$ | $0.7003$ | $0.7180$ | $0.4248$
Icosphere (ours) | $\mathbf{0.9352}$ | $0.9703$ | $\mathbf{0.8797}$ | $\mathbf{0.6890}$ | $\mathbf{0.7037}$ | $0.6970$ | $\mathbf{0.7562}$ | $\mathbf{0.5744}$
| Sofa | Partition | Bed | Chair | Table | Shelving | Chandelier | All Classes
Equirect. (gnom. kernel) [5, 14] | $0.7114$ | $0.4172$ | $0.7133$ | $0.4219$ | $0.4587$ | $0.3278$ | $0.4491$ | $0.5904$
Equirect. (equirect. kernel) | $0.7337$ | $0.4448$ | $0.7228$ | $0.4415$ | $\mathbf{0.5363}$ | $0.2843$ | $\mathbf{0.4566}$ | $0.5969$
Icosphere (ours) | $\mathbf{0.7374}$ | $\mathbf{0.4683}$ | $\mathbf{0.7776}$ | $\mathbf{0.4375}$ | $0.5018$ | $\mathbf{0.3733}$ | $0.4472$ | $\mathbf{0.6639}$
Table 1: Semantic segmentation mIOU for the different representations. Classes
ordered by descending frequency.

##  3 Semantic Segmentation Evaluation

We evaluate the icosphere representation using semantic segmentation and
demonstrate that limiting distortion in the representation spawns significant
accuracy improvements.

###  3.1 Dataset

For our experiments, we use the SUMO [15] dataset, a derivative of SunCG [12]
containing $58,631$ RGB-D cube-maps representing synthetic spherical image
captures. We divide the SUMO dataset into 80% training, 10% validation, and
10% testing splits and pre-process the data by resampling the images onto an
$7^{th}$ order icosphere. An example of this resampling is displayed with a
super-imposed convolutional kernel in Figure 1. As the dataset contains $132$
semantic classes, some with high inter-class similarity (e.g. ‘bunk bed’ and
‘baby bed’), we aggregated these similar classes to balance the data for these
experiments. There is still a very high variance in class frequency even after
aggregation, so we only evaluate performance on the $15$ most common classes.
For exact comparison with existing work, we render the predictions on the
icosphere back into an equirectangular image to assess per-pixel accuracy.

###  3.2 CNN Training

We follow the method of Coors _et al_. [5] and Tateno _et al_. [14] to sample
the spherical data using a gnomonic projection of the planar filter. However,
as our data is no longer in a 2D image, we project the planar filter onto the
icosphere’s surface and sample from the faces using barycentric interpolation.
We set the filter resolution to the mean distance between adjacent vertices on
the icosphere, to maintain analogy with the notion of resolution of the planar
filter being set to the distance between adjacent pixels. During convolution,
we apply this filter at each vertex of the icosphere and sample from the faces
using barycentric interpolation of the data at the $3$ vertices. For up- and
down-sampling, we leverage different subdivision orders, which, similar to
pixels in 2D images, vary roughly by a factor of 4 in the number of vertices.
We compare convolving with our icosphere representation to the method of [5,
14], which convolves with equirectangular images. We use the encoder-decoder
architecture given in Figure 3 for both methods. We train for 10 epochs using
Adam [8] with an initial learning rate of $10^{-4}$, reduced by half every 3
epochs. Our training criterion is the cross-entropy loss weighted by inverse
class frequency.

###  3.3 Results and Discussion

Results of our semantic segmentation experiments are given in Table 1 as the
per-class mean intersection-over-union (mIOU). While certain classes perform
marginally better in the equirectangular format, we see an $12.6\%$
improvement in the overall mIOU compared to [5, 14] simply by resampling the
image. There are two reasons for this jump. The first is that the gnomonic
projection used for the kernel in [5, 14] is not optimal for modeling the
distortion in an equirectangular image. In fact, an equirectangular projection
of the kernel would be a more appropriate, as it would correctly adapt the
sampling locations to equirectangular distortion. The difference in the
sampling patterns is slight, though, as shown in Figure 4, and the
quantitative result is only slightly better (second row in Table 1). It is
worth mentioning that the gnomonic kernel does make sense for sampling from
the icosphere, as it represents the projection of a tangent plane onto a
sphere. In our case, the planar filter is tangent to our icosahedral spherical
approximation. Nevertheless, the more important consideration is that neither
the gnomonic nor equirectangular filter projections correct for the
information imbalance caused by the distortion. The filter projection method
of [5, 14] dynamically adjusts where the kernel samples, but it only accounts
for the conformality and equidistance properties of distortion. It does not
address the areal inequality. Unequal area projections disrupt the
distribution of information throughout the image. For example, all of the
pixels along the top row in an equirectangular images are essentially
redundant, having all been sampled from the pole. That whole row carries the
same information capacity of a single pixel at the equator. For convolutions,
this means that the amount of information aggregated at each pixel varies
throughout the image. The ISEA projection, being nearly equal-area, corrects
for this, ensuring that predictions at each vertex of the icosphere are made
by accumulating equivalent quantities of information. This translates to
better inference. It is also worth noting that the objects for which the
equirectangular representation perform marginally better are typically found
at the top (ceiling, chandelier) or bottom (rug) of an image. The
equirectangular format’s pixel redundancy increases the number of training
samples for these objects, which would explain the better performance at test
time.

(a) Gnomonic

(b) Equirectangular

Figure 4: Comparing convolutional kernel sampling patterns using gnomonic and
equirectangular projections.

##  4 Conclusion

Borrowing from the field of cartography, we have identified three key
principles of local distortion that impact our ability to convolve on
spherical data. Applying these principles to our choice of data
representation, we have proposed the ISEA projection on an icosphere as an
ideal spherical image representation. We demonstrated that by significantly
limiting distortion, this representation provides a $12.6\%$ boost in semantic
segmentation performance.

## References

  * [1] Planet texture maps.  http://planetpixelemporium.com/earth8081.html.  Accessed: 2018-04-16.
  * [2] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289, 2015.
  * [3] Taco Cohen, Mario Geiger, Jonas Köhler, and Max Welling.  Convolutional networks for spherical signals.  arXiv preprint arXiv:1709.04893, 2017.
  * [4] Taco S Cohen, Mario Geiger, Jonas Köhler, and Max Welling.  Spherical cnns.  arXiv preprint arXiv:1801.10130, 2018.
  * [5] Benjamin Coors, Alexandru Paul Condurache, and Andreas Geiger.  Spherenet: Learning spherical representations for detection and classification in omnidirectional images.  In European Conference on Computer Vision, pages 525–541. Springer, 2018.
  * [6] R Buckminster Fuller.  Synergetics: explorations in the geometry of thinking.  Estate of R. Buckminster Fuller, 1982.
  * [7] Jon A Kimerling, Kevin Sahr, Denis White, and Lian Song.  Comparing geometrical properties of global grids.  Cartography and Geographic Information Science, 26(4):271–288, 1999\.
  * [8] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [9] Charles Loop.  Smooth subdivision surfaces based on triangles.  Master’s thesis, University of Utah, Department of Mathematics, 1987\.
  * [10] John Parr Snyder.  Map projections–A working manual, volume 1395.  US Government Printing Office, 1987.
  * [11] John P Snyder.  An equal-area map projection for polyhedral globes.  Cartographica: The International Journal for Geographic Information and Geovisualization, 29(1):10–21, 1992.
  * [12] Shuran Song, Fisher Yu, Andy Zeng, Angel X Chang, Manolis Savva, and Thomas Funkhouser.  Semantic scene completion from a single depth image.  IEEE Conference on Computer Vision and Pattern Recognition, 2017\.
  * [13] Yu-Chuan Su and Kristen Grauman.  Learning spherical convolution for fast features from 360 imagery.  In Advances in Neural Information Processing Systems, pages 529–539, 2017.
  * [14] Keisuke Tateno, Nassir Navab, and Federico Tombari.  Distortion-aware convolutional filters for dense prediction in panoramic images.  In European Conference on Computer Vision, pages 732–750. Springer, 2018.
  * [15] Lyne Tchapmi and Daniel Huber.  The sumo challenge.
  * [16] Bo Xiong and Kristen Grauman.  Snap angle prediction for 360° panoramas.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 3–18, 2018.
  * [17] Nikolaos Zioulis, Antonis Karakottas, Dimitrios Zarpalas, and Petros Daras.  Omnidepth: Dense depth estimation for indoors spherical panoramas.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 448–465, 2018.

Convolutions on Spherical Images: Supplementary Material

##  5 Qualitative Results

Figure 5 shows a qualitative comparison between semantic segmentation results
obtained using the equirectangular image representation and our proposed ISEA
representation on the icosphere. As this work is focused on examining the
methodology, the results are not post-processed in any way, and we note that
the architecture and loss we use are not state-of-the-art. However, the
representation of the data is network-agnostic, and we expect the outcome to
be the same regardless of the choice of network architecture and training
routine.

The quantitative results in Table 1 of the main paper illustrate that the ISEA
representation improves segmentation quality. This is borne out qualitatively
as well. In Figure 5, rows 1 and 2 provide clear examples of this effect. The
improvement comes from both more accurate pixel classification, as with the
door in the center-right of the image in row 3, and more precise segmentation
coverage, as with the stools in the center of the image in row 4. Rows 5 and 6
show examples where our representation has provided a better quantitative
result, but the actual segmentations are not noticeably better than with the
equirectangular image representation. Instead, we seem to trade some mistakes
for others. The last row of samples demonstrates that the equirectangular
image representation can achieve similar performance to the icosphere
representation. In fact, there are many examples where differences in the
qualitative results are indiscernible. Nonetheless, we conclude that there is
a clear and demonstrable benefit to using the ISEA representation, justified
by both the quantitative metrics and the qualitative comparison given here.

1)

2)

3)

4)

5)

6)

(b) Equirectangular Image (c) Icosphere (d) Ground Truth

7)

(a) Input

Figure 5: Sample semantic segmentation results. In column (b) are the results
of applying the method of [5, 14] to equirectangular image representations of
spherical data. In column (c) are the results of our method using the ISEA
representation on the icosphere. The icosphere results are resampled back to
an equirectangular projection for visualization here. Note that the
visualizations in column (c) have pixelated edges due to the use of nearest-
vertex interpolation in the resampling from the icosphere. This is an artifact
of the rendering, not a result of the predictions.


Conference Paper
2022

Shock System Dynamics of a Morphing Bump Over a Flat Plate

Ahmed A. Hamada1,2, Lubna Margha1,2, Mohamed M. AbdelRahman2, Amr Guaily3,2
1Department of Ocean Engineering, Texas A$\&$M University, College Station,
TX, 77843, USA
2Department of Aeronautical and Aerospace Engineering, Cairo University, Giza,
12613, Egypt
3Smart Engineering Systems Research Center (SESC), Nile University, Giza,
12588, Egypt

## Abstract

The shock wave boundary layer interaction (SW-BLI) phenomenon over transonic
and supersonic airfoils captured the attention of aerospace engineers, due to
its disastrous effect on the aerodynamic performance of these vehicles. Thus,
the scientific community numerically and experimentally investigated several
active and passive flow control elements to reduce the effect of the
phenomenon, such as vortex generator, cavity, and bump. They focused on
designing and optimizing the shape and location of the bump control element.
However, the transit movement of the bump from the state of a clean airfoil to
the state of an airfoil with a bump needs more investigation, especially the
dynamics of the shock system. Thus, it is preferred to start with simple
geometry, such as a flat plate, to fully understand the flow behavior with a
morphing bump. In this paper, the shock dynamics due to the movement of a bump
over a flat plate flying at supersonic speed are numerically investigated. The
bump is located at the impingement position of the shock wave and is moved at
different speeds. This study determines the suitable speed that achieves the
minimum entropy change, which is the representation parameter of the
transition period. The two-dimensional unsteady Navier-Stokes equations are
solved using OpenFOAM to simulate the flow field variables, while the motion
of the bump is tracked using the Arbitrary Lagrangian-Eulerian (ALE)
technique. The results show that a spatial lag on the shock system from the
steady-state solution occurs due to the movement of the bump. Further, the
spatial lag increases with the increase in the bump’s speed. This causes a
high increase in the flow parameters and consequently the total entropy
changes on the bump surface. Generally, it is common to move the bump over the
longest possible time to approximate a quasi-steady flow during the motion.
However, this causes a deviation in the flow parameters between the final time
of transition and the steady-state case of bump existence. Thus, it is
concluded that the optimal non-dimensional time for a morphing bump in a
supersonic flow of Mach number of 2.9 is 2, which is different than the
longest time of 10.

Keywords: Moving bump; Supersonic flow; Flat plate; Active flow control.

## 1\. Introduction

The Shock wave boundary layer interaction (SWBLI) over a high-speed wing
disastrously affects its aerodynamics performance and structural lifetime. The
phenomenon of SWBLI was experimentally observed for the first time by Ferri
[1] in 1939. Green [2] summarized the occurrence of the phenomenon in four
conditions: externally as in transonic airfoils, and near control flaps, or
internally as in high-speed inlets (scram-jets), and nozzles at off-design.
Based on strength of the shock wave, the

phenomenon varies the height of the boundary layer, the surface drag of the
body, and/or the heat transfer. The SWBLI adversely affects the structure and
geometry of the flying vehicles [3]. The flow experiences a high instantaneous
increase in pressure and thermal transfer due to the existence of SWBLI, which
leads to a reduction in the fatigue life of the structure. Consequently, that
imposes strong constrictions on choosing the material of the structure,
resulting in an expensive and heavy design. Moreover, this phenomenon happens
in transonic airfoils after a critical Mach number, causing the sudden
increase in total drag, which is called “drag divergence”.

There are several ways to eliminate or reduce the SWBLI effect which is
discussed in detail by Dennis Bushnell [4]. The focus in this paper will be on
the bump approach. In 1922, Ashill et al. [5] were the first researchers who
used a two-dimensional bump at the upper surface of airfoils to affect the
strength of incurring shock waves. Milholen et al. [6] and Patzold et al. [7]
proved that using bumped airfoils increases the lift, reduces the drag, and
postpones the buffeting, which enhances the aerodynamic performance. Fulker
[8], who is a DASA-Airbus researcher, showed that applying a bump in A-340’s
hybrid laminar wing will save fuel by $2.11\%$ at $M_{\infty}=0.84$. Several
investigations are conducted to conclude with the optimal shape and location
of the bump [9, 10, 11, 12]. Further, Eastwood and Jarrett [13] investigated
the three-dimensional bump control technique. Yun Tian et al. [14] showed that
airfoil geometry and free stream condition determine the bump’s optimal
geometric parameters which specify its shape and location over the airfoil. In
addition, varying these geometric parameters highly affects the aerodynamic
performance.

By the beginning of the 21st century, the morphing wing concept become
industrial applicable. Stanewsky et al. discussed the design aspects of the
morphing concept from the view of the civil aircraft industry [15]. From that
time, many scholars started to investigate the flow behavior over the morphing
wing at different velocity regimes [16, 17, 18, 19]. Further, Bruce and
Colliss [20] conducted a review study for the shock control bumps, showing
that the morphing bump is a future ideal solution to limit the SWBLI effect on
wings at transonic flows.

To simplify the phenomenon, consider a strong two-dimensional shock wave that
impinges with a boundary layer of a flat plate, see Figure 1. Downstream the
impinging region, separation, and relaxation zones are generated. The flow
changes from a non-equilibrium state in the separation zone to an equilibrium
state in the relaxation zone. Through these zones, the flow experiences a
quick increase in pressure and heat transfer, detaching, and reattaching the
boundary layer. Moreover, the turbulent kinetic energy increases due to the
generation of the strong shear and the adverse pressure gradient in and near
the separation zone, respectively. During the phenomenon, the flow contains
different types of waves generated from the SW-BL.I, such as; a reflected
shock wave, expansion waves, and compression waves [21].

The static condition for the bump control element is fully investigated.
Particularly, they concentrated on designing and optimizing the bump control
element’s shape and position. However, additional detailed investigations are
needed for the transition of the morphing bump from a clean airfoil to an
airfoil with a bump, focusing on the dynamics of the shock system. The aim of
this research is to numerically study the transient effect of the morphing
pump over a simple geometry (flat plate) on the unsteady flow features of the
shock dynamics. We started to study the transient phenomenon over a simple
geometry first to fully understand the transition effects. The location of the
bump is chosen to be at the impingement location of the incident shock wave.
Then, starting the unsteady motion of the bump at different speeds and a
constant velocity profile. This enabled us to determine the suitable morphing
bump’s speed that would minimize the entropy change.

Figure 1: A schematic of 2D SW-BLI separation flow [21].

## 2\. Physical Model

### 2.1 Model Description

The flow configuration of the morphing bump in a supersonic flow is shown in
Figure 2. Further, the variation of the morphing bump shape with time is shown
in Figure 3. The height of the computational inflow boundary is $1.1m$, and
the height of the outflow at the right is $1m$. The difference in height
between the two edges represents the height of the wedge, which is the shock-
source surface. The wedge is far from the inlet with a length of $0.2656m$ and
the wedge stream-wise length is $0.5173m$. The morphing bump is placed at the
center of the bottom surface, representing the origin of the domain. The
origin is $2.25m$ far away from the inlet. The final shape of the morphing
bump is a parabolic arc of height-to-length ratio, $\alpha=4\%$, and
represented by Equation (1).

$y=\sqrt{\left(\frac{\alpha^{2}+0.25}{2\alpha}\right)^{2}-x^{2}}+\frac{\alpha^{2}-0.25}{2\alpha}$
(1)

When the supersonic flow with a free-stream Mach number $M_{\infty}=2.9$ and a
free-stream Reynolds number $Re_{\infty}=6.6\times 10^{7}$ impinges the wedge,
an incident shock wave (I) is generated hitting the bottom surface. Then, a
reflected shock wave (R) appears. Further, a series of expansion fans
generates from the wedge trailing edge (E). The expansion fan interacts with
the shock system causing a curvature to them. The remaining top surface after
the wedge is considered an outlet to enable the flow to go out freely. The
morphing bump rises to its final shape with a constant velocity profile and
different speeds (different motion periods of $0.1$, $1.0$, $5.0$, $7.5$,
$10.0$). The initial condition of the problem is shown in Figure 2 (a) with
the solid lines, while the dashed lines are for the final condition. The
values of the problem’s parameters are shown in Table 1. The appearance of the
morphing bump causes a generation of two oblique shocks at its leading edge
(L) and at its trailing edge (T). Further, a series of expansion fans appear
over the morphing bump due to its geometric curvature (Series of Es).

Figure 2:  Schematics of shock system over the morphing bump. Figure 3:  Shape of the morphing bump with time. Table 1: System properties and parameters. Mach number, $M_{\infty}$ | $2.9$
---|---
Reynolds number, $Re_{\infty}$ | $6.6\times 10^{7}$
Inlet length, $l_{i}$ | $1.1m$
Wedge chord, $c_{w}$ | $0.527m$
Wedge angle, $\alpha_{w}$ | $11^{\circ}$
Slip wall length, $l_{sl}$ | $0.26555m$
Bump final length, $l_{b}$ | $0.5m$
Bump final height, $h_{b}$ | $0.01m$
Motion periods, $t^{\ast}_{f}$ |  $0.1$, $1$, $1.8$, $2$, $3$, $5$, $7.5$, $10$

### 2.2 Governing Equations

Two-dimensional unsteady compressible Navier-Stokes equations are used to
model the supersonic flow, and are expressed as:

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle+\nabla\cdot\left(\rho\textbf{{V}}\right)=0$ (2a)
$\displaystyle\rho\left(\frac{\partial\textbf{{V}}}{\partial
t}+\left(\textbf{{V}}\cdot\nabla\right)\textbf{{V}}\right)$
$\displaystyle=-\nabla p+\nabla\cdot$
$\displaystyle\left(\mu\left(\left(\nabla\textbf{{V}}\right)+\left(\nabla\textbf{{V}}\right)^{T}\right)\right)-\frac{2}{3}\nabla\left(\mu\nabla\cdot\textbf{{V}}\right)$
(2b) $\displaystyle\rho\left(\frac{\partial e}{\partial
t}+\left(\textbf{{V}}\cdot\nabla\right)e\right)$
$\displaystyle+p\left(\nabla\cdot\textbf{{V}}\right)$
$\displaystyle=-\frac{2}{3}\mu$
$\displaystyle\left(\nabla\cdot\textbf{{V}}\right)^{2}+\mu\left(\left(\nabla\textbf{{V}}\right)+\left(\nabla\textbf{{V}}\right)^{T}\right):\left(\nabla\textbf{{V}}\right)$
(2c)

where $\rho$ is the density, $V$ is the velocity vector, $p$ is the pressure,
$e$ is the specific energy per unit mass, $\mu$ is the fluid dynamic
viscosity.

## 3\. Computational Model

### 3.1 Numerical Implementation

rhoCentralDyMFoam is a toolbox flow solver in OpenFOAM®$-$v2006. It solves
transient, compressible flow problems using a density-based approach,
supporting dynamic mesh applications. The semi-discrete and upwind-central
non-staggered schemes of Kurganov and Tadmor are combined in the solver [22,
23], which uses an operator-splitting method [24]. The van Leer limiter is
implemented in the solver to balance between shock capture, oscillations-free
fields, and computational cost during the solution [25].

The letters ”DyM” in the name of the rhoCentralDyMFoam solver indicates that
the dynamic mesh simulation is applicable with a dynamic FV mesh technique. In
this study, the motion of the bump will be actively controlled. Thus, the
coded motion solver is used in this study to arbitrarily move the
computational mesh inside the domains, optimizing the cells’ shape. Hence, the
dynamic FV mesh follows the Arbitrary Lagrangian-Eulerian technique, which is
presented by Hirt et al. [26].

### 3.2 Computational Domain

A body-fitted computational domain is constructed without the existence of the
bump. The domain is clustered in the direction of the $y$-axis along the
bottom surface and also clustered in the $x$-direction at the area of the
morphing bump. The domain was divided into 10 blocks where the cells are
nearly orthogonal and quadrilateral. The grid will be updated each time step
during the motion of the bump, preserving the body-fitted mesh by applying the
Arbitrary Lagrangian-Eulerian (ALE). Four types of boundaries define the
borders of the domain: inlet, outlet (right and top after the wedge), non-slip
walls (bottom and wedge), and slip wall (top before wedge). Figure 4 shows the
boundary and initial conditions over the Schematic of the grid.

A mesh-independent study is essential to ensure the result’s accuracy. The
study was implemented over the case of without bump, starting with mesh 1 of
size $100\times 45$ cells. Then, high-density meshes (2, 3, and 4) were
generated by doubling the number of cells in $x$ and $y$ directions. The time
step was controlled with the Courant–Friedrichs–Lewy (CFL) number of $0.2$.
Figure. 5 shows the variation of the pressure distribution at height of
$y=0.5m$ for different mesh sizes. The results with mesh 3 show a converged
solution that does not depend on the mesh quality. This is because the
difference between mesh 3 and mesh 4 is negligible. Thus, Mesh 3 of size
$400\times 180$ was chosen to implement all the simulations to minimize
computational time. The minimum element size in mesh 3 is $2.8mm\times 1.2mm$,
and the time step to ensure the CFL number of $0.2$ is $111.1\mu s$.

Figure 4:  Schematics of the computational domain for mesh 1. Figure 5: Mesh
independent study for a supersonic flow over a clean flat plate, by comparing
the pressure distribution at the height, $y=0.5m$.

## 4\. Results and Discussion

### 4.1 Shock Structure over a Stationary Bump

Firstly, the benefit of using the bump in the SWBLI problem has to be shown
before investigating the transitional effect of the morphing bump. Thus,
supersonic flows over a flat plate without and with a stationary bump were
simulated to obtain the steady-state solutions of the morphing bump limiting
positions, as shown in Figure 6. The benefit of using the bump is indicated by
plotting the pressure distribution over the lower boundary for both problems;
without and with a stationary bump, as shown in Figure 7.

The existence of the bump created a series of weak expansion fans which in
total decreased the pressure dramatically to the inlet value of pressure as
shown in Figure 7 within the region of $x\approx(0,0.5)$. Then, the pressure
increased at the end of the bump to the original case (clean flat plate) at
$x=0.5$, due to the generation of the trailing oblique shock. This decrease of
the pressure in this simple geometry (flat plate) case shows that the bump is
applicable to control the SWBLI in the airfoil case, where the pressure
decreases over the upper surface, and accordingly the lift force increases.
This aligns with the results of Mazaheri et al. [12].

(a) Without a bump.

(b) With a stationary bump.

Figure 6:  Pressure contours of supersonic flow over a flat plate with and
without a bump at dimensionless time of 10. Figure 7: Comparison between the
two problems, supersonic flow over a flat plate with and without a bump using
the pressure distribution over the bottom boundary.

(a) $t^{\ast}=0$

(b) $t^{\ast}=2$

(c) $t^{\ast}=4$

(a) $t^{\ast}=6$

(b) $t^{\ast}=8$

(c) $t^{\ast}=10$

Figure 9:  Pressure contours, shown in (a), zoomed picture of the morphing
mesh, shown in (b), and entropy, temperature, pressure, and density
distributions, shown in (c) and (d), of supersonic flow over a flat plate with
a morphing bump control element with constant velocity profile along the
Bottom boundary.

### 4.2 Shock Dynamics over Morphing Bump

The transition phase of the morphing bump was investigated in detail with a
constant velocity profile. The motion of the morphing bump occurred during
different dimensionless times, $t^{\ast}_{f}=10$, $7.5$, $5$, $3$, $2$, $1.8$,
$1$, and $0.1$. This corresponds to a constant local Mach number, $M_{b}$, of
$0.004$, $0.00533$, $0.008$, $0.0133$, $0.02$, $0.0222$, $0.04$, and $0.4$ at
the mid-point of the morphing bump ($x=0$), respectively, and zero-velocity
points at the two edges of the bump. Figure 9 shows the transitional solution
of the morphing bump with the slowest constant velocity (the Mach number at
$x=0$ is $0.004$ and the final dimensionless time, $t^{\ast}_{f}$ is $10$) by
plotting the pressure contours, shown in (I), zoomed picture of the morphing
mesh, shown in (II), and entropy, temperature, pressure, and density
distributions, shown in (III) and (IV), along the bottom boundary at different
dimensionless times, $t^{\ast}$.

(a) Stationary bump and morphing bump at $t^{\ast}_{f}=10$.

(b) Morphing bump at $t^{\ast}_{f}=10$.

(c) Morphing bump at $t^{\ast}_{f}=1$.

(d) Morphing bump at $t^{\ast}_{f}=0.1$.

Figure 10: Pressure contours of supersonic flow over a flat plate with a
stationary/morphing bump at different constant velocities, to show the spatial
lag in the shock system.

There is a spatial lag in the leading and trailing oblique shocks between the
transitional solution at different values of constant velocity (different
values of the final dimensionless times, $t^{\ast}_{f}$ ), and the stationary
condition, as shown in Figure 10. This spatial lag in the shock system from
the stationary condition increases when the value of the morphing bump’s
constant velocity increases, i.e. the operational time of the morphing bump
was reduced. This was concluded also in a rotating wedge problem, made by
Margha et al. [27]. The spatial lag was clearly indicated by comparing the
pressure, and Mach number distributions at a height of $y=1$, between the
stationary condition and the morphing conditions at different $t^{\ast}_{f}$,
as shown in Figure 11. Further, the change in the spatial lag can be neglected
when the operational time of the morphing bump exceeds the final dimensionless
time, $t^{\ast}_{f}=5$.

(a) Pressure distribution.

(b) Mach number distribution.

Figure 11: Pressure, and Mach number distributions of supersonic flow over a
flat plate with a stationary/morphing bump control element at a height of
$y=1m$, to show the change of spatial lag in the shock system at different
dimensionless morphing periods.

A comparison between the start state (clean flat plate) and the final state
(at the end of the morphing motion) over the bump’s surface at different
motion periods was performed by obtaining the flow parameter distributions.
Then, the temporal change in the entropy difference distribution with the case
of a clean flat plate problem at a steady state over the lower boundary,
$\left(\Delta s_{@t=t^{\ast}_{f}}(x)-\Delta s_{@t=0}(x)\right)/c_{v}$, is
calculated, as shown in Figure 12. Further, the entropy difference
distribution, $\Delta s(x)$, is defined as the difference between the entropy
at a certain $x$ location with the inlet condition.

(a) Pressure distribution.

(b) Temporal change in the entropy difference distribution.

Figure 12:  Pressure distribution, and the temporal change in the entropy
difference distribution for the supersonic flow over a flat plate with a
morphing bump over the lower boundary.

When the leading oblique shock appears due to the morphing bump motion, the
pressure and the density over the bump increase. Thus, the temporal change in
the entropy difference rises up at the bump’s leading edge. Then, the series
of expansion, before and after the incident oblique shock, decrease the
pressure and density along the morphing bump. This may cause negative values
for the temporal change in the entropy difference over the bump at relatively
low morphing speeds. Furthermore, the location of the incident shock moves
slightly to the left during the upward motion of the bump. The trailing
oblique shock again increases the temporal change in the entropy difference.
In addition, there is a deviation between the two entropy difference over the
lower boundary; the case of the stationary bump at steady state, and the case
of the morphing bump at final states of the motion, see Figure 9 (f) and
Figure 10 (a). The deviation happens due to the effect of spatial lag on the
shock system, which results from the bump’s motion.

(a) Change in temporal entropy difference.

(b) Change in temporal entropy difference.

Figure 13:  Pressure distribution, and the change in temporal entropy
difference distribution for the supersonic flow over a flat plate with a
morphing bump over the lower boundary.

when the morphing bump is moved with a relatively high speed (small
$t^{\ast}_{f}$), we lose the gain in the pressure difference
($p_{@t=t^{ast}_{f}}(x)-p_{t=0}(x)$), that is obtained from the series of
expansion fans. This is explained as the pressure applied on the morphing bump
would increase, which can be considered as a resistant action to the motion of
the bump. Further, the compressibility effect increases with the bump’s speed,
and consequently the density increases. When the air is compressed with the
motion of the bump, the molecular distances decrease, and then the friction
between adjacent molecules increases. Thus, the temperature increases with the
bump’s motion. Thus, when the morphing bump moves faster, the flow parameters
(pressure, density, and temperature) increase, as shown in Figure 12 (a). This
is clearly shown in the temporal change in the entropy difference distribution
of Figure 12 (b). Hence, the results show that the slowest velocity for the
morphing bump, which approximates a quasi-steady flow during the motion,
decreases the lag effect.

Despite the existence of the lag effect in the shock system, the slowest
tested morphing bump’s speed ($t^{\ast}_{f}=10$) is not the optimal one, due
to the deviation in the entropy from the stationary bump at steady state, see
Figure 9 (f). This deviation represents the losses in the shock system,
resulting from the morphing speed. Thus, comparing the entropy difference for
each speed was made with the stationary bump at a steady state, as shown in
Figure 13. The results show the suitable motion period, $t^{\ast}_{f}$, for
the morphing bump at a free-stream Mach number of 2.9 is 2. The reason is that
its entropy difference is very close to that of the stationary bump at a
steady state. Furthermore, the lag effect in the system resulting from that
morphing speed is within acceptable moderate levels, as shown in Figure 11.

## 5\. Conclusion

The current research work aimed to investigate the morphing of a bump control
element over a flat plate at a free-stream Mach number of $2.9$ and a high
Reynolds number of $6.6\times 10^{7}$. Different values of constant velocity
were conducted to study the effect of the morphing bump on the shock system
with non-dimensional motion periods, $t^{\ast}_{f}$ of $0.1$, $1.0$, $1.8$,
$2.0$, $3.0$, $5.0$, $7.5$, $10.0$. Further, the steady state of supersonic
flows over a clean flat plate and a flat plate with a stationary bump were
conducted to show the beneficial effect of the bump’s existence. Furthermore,
a comparison between the dynamic and static cases was achieved. The results
showed that a spatial lag in the shock system appears due to the dynamic
motion. In addition, the lag effect remarkably increases when the local Mach
number at the tip of the morphing bump, $M_{b}$, increases higher than
$0.008$. The relatively fast morphing bump would increase the compressibility
effect in the near area which may compensate for the beneficial effect of the
bump’s existence. Besides, the relatively slow morphing speed results in a
deviation in the entropy from the stationary bump case, which is a
representation of flow momentum losses. Thus, the suitable speed to morph with
is the one that results in neither a remarkable lag effect in the shock system
nor high losses in the entropy deviation from the stationary steady-state
case. For the case of supersonic flow with $M_{\infty}=2.9$ over a flat plate,
the suitable bump’s morphing period, $t^{\ast}_{f}$, was found to be 2. For
future work, various velocity profiles with different speed values are
recommended to be tested for the morphing bump. This will determine the
suitable morphing bump’s velocity profile and speed, that achieves low time-
average entropy for the flat plate problem. Then, this configuration of motion
will be tested over a transonic airfoil to achieve a high time-averaged lift-
to-drag ratio.

## References

  * Ferri and ATTI [1939] Antonio Ferri and DG ATTI.  Experimental results with airfoils tested in the high speed tunnel at guidonia: Naca tm 946.  Technical report, Washington, DC: National Advisory Committee for Aeronautics, 1939.
  * Green [1970] JE Green.  Interactions between shock waves and turbulent boundary layers.  _Progress in Aerospace Sciences_ , 11:235–340, 1970.
  * Dolling [2001] David S Dolling.  Fifty years of shock-wave/boundary-layer interaction research: what next?  _AIAA journal_ , 39(8):1517–1531, 2001.
  * Bushnell [2004] Dennis M Bushnell.  Shock wave drag reduction.  _Annu. Rev. Fluid Mech._ , 36:81–96, 2004.
  * Ashill and Fulker [1992] PR Ashill and JL Fulker.  92-01-022 a novel technique for controlling shock strength of laminar-flow aerofoil sections.  _DGLR BERICHT_ , pages 175–175, 1992.
  * Milholen and Owens [2005] William Milholen and Lewis Owens.  On the application of contour bumps for transonic drag reduction.  In _43rd AIAA Aerospace Sciences Meeting and Exhibit_ , page 462, 2005\.
  * Pätzold et al. [2006] Martin Pätzold, Thorsten Lutz, Ewald Kramer, and Siegfried Wagner.  Numerical optimization of finite shock control bumps.  In _44th AIAA Aerospace Sciences Meeting and Exhibit_ , page 1054, 2006.
  * Fulker [1999] J Fulker.  The euroshock programme (a european programme on active and passive control of shock waves).  In _17th Applied Aerodynamics Conference_ , page 3174, 1999.
  * Sommerer et al. [2000] Andreas Sommerer, Thorsten Lutz, and Siegfried Wagner.  Numerical optimisation of adaptive transonic airfoils with variable camber.  In _Proceedings of the 22nd International Congress of the Aeronautical Sciences_ , volume 27. Optimage Ltd. Edinburgh, UK, 2000.
  * Wong et al. [2007] WS Wong, A Le Moigne, and N Qin.  Parallel adjoint-based optimisation of a blended wing body aircraft with shock control bumps.  _The Aeronautical Journal_ , 111(1117):165–174, 2007.
  * Tian et al. [2014] Yun Tian, PeiQing Liu, and Zhi Li.  Multi-objective optimization of shock control bump on a supercritical wing.  _Science China Technological Sciences_ , 57(1):192–202, 2014.
  * Mazaheri et al. [2015] K Mazaheri, KC Kiani, A Nejati, M Zeinalpour, and R Taheri.  Optimization and analysis of shock wave/boundary layer interaction for drag reduction by shock control bump.  _Aerospace Science and Technology_ , 42:196–208, 2015.
  * Eastwood and Jarrett [2012] Jeremy P Eastwood and Jerome P Jarrett.  Toward designing with three-dimensional bumps for lift/drag improvement and buffet alleviation.  _AIAA journal_ , 50(12):2882–2898, 2012.
  * Tian et al. [2011] Yun Tian, PeiQing Liu, and PeiHua Feng.  Shock control bump parametric research on supercritical airfoil.  _Science China Technological Sciences_ , 54(11):2935, 2011.
  * Monner et al. [2000] HP Monner, E Breitbach, Th Bein, and H Hanselka.  Design aspects of the adaptive wing—the elastic trailing edge and the local spoiler bump.  _The Aeronautical Journal_ , 104(1032):89–95, 2000.
  * Olivett et al. [2021] Anthony Olivett, Peter Corrao, and M Amin Karami.  Flow control and separation delay in morphing wing aircraft using traveling wave actuation.  _Smart Materials and Structures_ , 30(2):025028, 2021.
  * Elbadry et al. [2021] Yusuf T Elbadry, Amr Gamal Guaily, Mohammed A Boraey, and Mohamed M AbdelRahman.  Active morphing control of airfoil at low reynolds number using level-set method.  In _2021 3rd Novel Intelligent and Leading Emerging Sciences Conference (NILES)_ , pages 265–268. IEEE, 2021.
  * Botez et al. [2018] RM Botez, A Koreanschi, OS Gabor, Y Tondji, M Guezguez, JT Kammegne, LT Grigorie, D Sandu, Y Mebarki, M Mamou, et al.  Numerical and experimental transition results evaluation for a morphing wing and aileron system.  _The Aeronautical Journal_ , 122(1251):747–784, 2018.
  * Popov et al. [2010] Andrei V Popov, Lucian T Grigorie, Ruxandra Botez, Mahmood Mamou, and Youssef Mébarki.  Real time morphing wing optimization validation using wind-tunnel tests.  _Journal of Aircraft_ , 47(4):1346–1355, 2010\.
  * Bruce and Colliss [2015] PJK Bruce and SP Colliss.  Review of research into shock control bumps.  _Shock Waves_ , 25(5):451–471, 2015.
  * Ma et al. [2014] Li Ma, Lipeng Lu, Jian Fang, and Qiuhui Wang.  A study on turbulence transportation and modification of spalart–allmaras model for shock-wave/turbulent boundary layer interaction flow.  _Chinese Journal of Aeronautics_ , 27(2):200–209, 2014.
  * Kurganov and Tadmor [2000] Alexander Kurganov and Eitan Tadmor.  New high-resolution semi-discrete central schemes for hamilton–jacobi equations.  _Journal of Computational Physics_ , 160(2):720–742, 2000.
  * Kurganov et al. [2001] Alexander Kurganov, Sebastian Noelle, and Guergana Petrova.  Semidiscrete central-upwind schemes for hyperbolic conservation laws and hamilton–jacobi equations.  _SIAM Journal on Scientific Computing_ , 23(3):707–740, 2001.
  * Greenshields et al. [2010] Christopher J Greenshields, Henry G Weller, Luca Gasparini, and Jason M Reese.  Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows.  _International Journal for Numerical Methods in Fluids_ , 63(1):1–21, 2010.
  * Marcantoni et al. [2012] Luis F Gutiérrez Marcantoni, José P Tamagno, and Sergio A Elaskar.  High speed flow simulation using openfoam.  _Mecánica Computacional_ , 31(16):2939–2959, 2012.
  * Hirt et al. [1974] Cyrill W Hirt, Anthony A Amsden, and JL Cook.  An arbitrary lagrangian-eulerian computing method for all flow speeds.  _Journal of computational physics_ , 14(3):227–253, 1974.
  * Margha et al. [2021] Lubna Margha, Ahmed A Hamada, Doyle D Knight, and Ahmed Eltaweel.  Dynamic transition from regular to mach reflection over a moving wedge.  In _ASME International Mechanical Engineering Congress and Exposition_ , volume 85581, page V004T04A006. American Society of Mechanical Engineers, 2021.