See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/359877788

A new strategy for solving store separation problems using OpenFOAM

Article in Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering · April 2022
DOI: 10.1177/09544100221080771

CITATIONS READS

3 2,011

3 authors, including:

Saleh Abuhanieh
TAI - Turkish Aerospace Industries, Inc.
6 PUBLICATIONS 24 CITATIONS

SEE PROFILE

All content following this page was uploaded by Saleh Abuhanieh on 19 June 2022.

The user has requested enhancement of the downloaded file.

Original Article

Proc IMechE Part G:

J Aerospace Engineering
A new strategy for solving store separation 2022, Vol. 0(0) 1–15
© IMechE 2022

problems using OpenFOAM Article reuse guidelines:

sagepub.com/journals-permissions
DOI: 10.1177/09544100221080771
journals.sagepub.com/home/pig

Saleh Abuhanieh1,2 , Hasan U. Akay1 and Barış Bicer2

Abstract
The ability of OpenFOAM to solve the problem of a store separating from an air vehicle (store separation problem) has
been evaluated using a dynamic mesh (Overset/Chimera) technique for an industry-class (transonic and generic)
benchmark test case. The major limitations of the standard libraries have been determined. To tackle these challenges,
a new strategy has been proposed and implemented using only open-source libraries and tools. The strategy combines
porting, modifying, and adapting an overset library from the OpenFOAM fork platform (foam-extend) to the standard
OpenFOAM platform (ESI). Furthermore, in order to overcome the well-known weakness of the standard OpenFOAM
compressible solvers, the newly adapted overset library was integrated with an open-source, density-based, and coupled
solver (HiSA), which uses the OpenFOAM technology. Additionally, a force restrained model was developed to consider
the externally applied forces on the store by the store ejectors. The accuracy of the developed strategy has been compared
with wind tunnel tests and the solutions of two well-known commercial codes, showing good agreements with them. While
the study has focused on simulations with inviscid Euler equations (typical of the test case considered here), the viscosity
effect on the solution has also been studied with Navier–Stokes equations and compared with other results in the literature,
showing minor differences. To the best of the authors’ knowledge, this is the first work which studies and validates the
store separation problem in transonic regime with OpenFOAM.

Keywords
trajectory prediction, large mesh movement, overset method, compressible flow, parallel computing, open-source tools

Date received: 12 July 2021; revised: 11 December 2021; accepted: 25 January 2022

Introduction freedom (6-DOF) solver is used to obtain the linear dis-

placements, linear velocities, angular displacements, and
Safe separation of stores is a crucial mission for many air angular velocities. For the on-line approach, a transient
vehicles. Predicting the trajectory of the store (to decide simulation which utilizes a dynamic mesh technique is re-
the safe separation envelop) can be done using wind tunnel quired. Since the resulting deformations are normally large,
testings1 and flight tests.2 However, both methods are the suitable dynamic mesh techniques are either based on
costly and raise some safety issues. The recent develop- a deformable mesh (with re-meshing capabilities)6,7 or an
ments in computer hardware as well as parallel Compu- overset/chimera8–11 technique. At each time step, the flow
tational Fluid Dynamics (CFD) algorithms have made fields are computed, the forces and moments are calculated,
computer simulations possible and more feasible. This and using a 6-DOF solver, the displacements and velocities
reduced the need for these methods. are calculated. Finally, the mesh is moved (or deformed/
There are two CFD approaches to solve this problem, re-meshed) according to these displacements.
which are inherited from the wind tunnel testing.3 The first In general, if the required number of cases to decide the
approach is the off-line (or grid survey)4 and the second safe separation envelop is enormous for the same ge-
approach is the on-line (or captive trajectory simulation).5 ometry, using the first approach by generating large da-
In the off-line approach, an aerodynamic database of tabase will be efficient. However, if the number of cases
forces and moments are obtained by solving the flow fields
(as a steady-state case) over a static mesh for different
scenarios. The scenarios (or grid test matrix) are mixture 1
Department of Mechanical Engineering, Atilim University, Turkey
2
of different Mach numbers, altitudes, angle of attacks, side Turkish Aerospace, OpenSource CFD Group, Turkey
angles, store positions, and orientations. Finding an op-
Corresponding author:
timum set of scenarios is normally obtained using the Saleh Abuhanieh, Turkish Aerospace, OpenSource CFD Group, Ankara,
Modern Design of Experiments (MDOE) method.3 After Turkey.
creating the aerodynamic database, a six degrees of Email: [email protected]
2 Proc IMechE Part G: J Aerospace Engineering 0(0)

are not so many, and for different geometries, the second Additionally, a dynamic mesh model driven by the
approach can be more efficient. Moreover, the second 6-DOF equations of motion, described below, is to be
approach is time-accurate. In this study, the second ap- integrated with the CFD model
proach is used. X d
Although OpenFOAM1 is a very popular and suc- F¼ ðmV Þ (4)
cessful open-source platform for CFD, it has been rarely dt
used for solving the store separation problem. The work of X d
Wadibhasme12 can be mentioned, where the Mesquite M¼ ðH Þ (5)
dt
dynamic mesh library of Menon13 was used, which was
available under OpenFOAM in earlier releases. However, where F and M are the applied force and moment vectors
Wadibhasme solved only a sample projectile case using an computed at the center of gravity of the store from the
incompressible solver (pimpleDyMFoam) for demon- CFD analysis, mV is the linear momentum vector, and H is
stration. The reasons that OpenFOAM is not used widely the angular momentum vector. The linear and angular
in store separation problems are twofold. Firstly, store displacement vectors of the store are then computed by
separation analysis is normally required in transonic and integrating equations (4) and (5), respectively, at each time
supersonic regimes, whereas in OpenFOAM, there is step.
a well-known limitation in its standard compressible
solvers, such that there are no density-based coupled
The overset/chimera method
compressible solvers, which are normally used in these
regimes. Secondly, the OpenFOAM dynamic mesh li- Since the first time it was proposed to the aerospace
braries have limited capabilities in this area. community four decades ago,14 the overset method proved
Because of these limitations, OpenFOAM was not used to be a very useful technique in CFD. Its basic idea is to
in the past for store separation simulations, neither was assemble the computational domain using separate
compared before to commercial codes which are used meshes (sub-domains). It has been mainly used to solve
extensively for this kind of problems. Therefore, this work two problems. Firstly, to simplify the meshing of complex
carried out to address this gap in the existing literature. geometries. In this case, each mesh can be prepared alone,
Accordingly, the main objects of this study are as follows: and that allows for generation of high-quality meshes
much easier including the block-structured meshes. Sec-
(a) To develop effective and accurate strategy in ondly, to simulate the cases where a solid body is moving
OpenFOAM for solving the store separation problem. inside a computational domain. The latter capability has
(b) To compare the developed strategy with the standard been used widely to solve the store separation problems.
OpenFOAM. An example is shown in Figure 1, where a pitching airfoil
(c) To compare the developed strategy with commercial case is presented. The airfoil mesh (the overset mesh is in
codes. red) is prepared separately and merged/overlapped with
(d) To analyze the viscosity effect on store separation at a background mesh (in blue). This configuration allows
the transonic regime. the airfoil to pitch or to move significantly during the
simulation without any need for re-meshing.
The Eglin1 test case has been used for code validation The main task in the overset method is to connect the
and comparisons due to the availability of the wind tunnel multiple meshes in a single computational domain where
results. the discretized governing equations are solved. This
process has been commonly termed in the literature as the
Overset Grid Assembly (OGA).15–18 In OpenFOAM, the
Problem formulation same has been called as the cell-to-cell mapping.2
Governing equations. The problem can be described as The OGA process can be divided into the following steps:
solving a transient, compressible, and viscous flow over
a moving body. The governing equations can be written as

∂½ρ
þ =  ðρvÞ ¼ 0 (1)
∂t
∂ðρvÞ
þ =  fρvvg ¼ ð=pÞ þ =  fτg (2)
∂t
∂½ρe
þ =  ðρveÞ ¼ =  q þ fσg : f=vg (3)
∂t
where ρ is the density, t is the time, v is the velocity vector,
p is the pressure, τ is the viscous stress tensor, e is the
specific internal energy, q is the heat flux vector, and σ is Figure 1. The pitching airfoil case. The background mesh is
the mechanical stress tensor. in blue and the overset mesh is in red.
Abuhanieh et al. 3

Figure 3. The interpolation stencil of cell I0 consists of the

donors D0  D4.

Figure 2. The overset DCI for the pitching airfoil case. The The next step, which is normally at the linear system
overset/store mesh (a) and the background mesh (b). The solver level, is to use the interpolation stencil to evaluate
calculated cells are in blue, the interpolated cells are in yellow,
the required fields values at each interpolated cell. Dif-
and the hole cells are in red.
ferent interpolation schemes can be used, for instance, the
inverse distance is expressed as:
(i) Hole identification: Finding the cells which rep-
resent the solid bodies or the cells which are located X
ND
 
ξ I0 ¼ ðωDi Þ ξ Di (6)
outside the computational domain. These cells are i¼1
excluded from the calculation (hole cells2 in
OpenFOAM). where ξ I0 is the value of the field variable ξ at the in-
(ii) Fringe construction: Deciding the cells which shall terpolated cell I0, ND is the total number of donor cells for
receive the information from the other mesh(s). cell I0, ωDi is weight contribution by the donor cell Di, and
These cells are called the fringe, receptors, acceptors, finally, ξ Di is the value of ξ at the donor cell Di. The weight
or in OpenFOAM, they are called the interpolated is defined as:
cells2.  
1
(iii) Donors search: Identifying the cells which deliver dD j
ωDi ¼ N  (7)
the information to the interpolated cells (the donor P 1 
D

cells2 in OpenFOAM). Normally, those are the cells dDj

j¼1
which their centers are the closest to the interpolated
cells centers; however, other criteria can be used too. where dDi is the distance between the donor Di cell center
to the interpolated cell center, the denominator is the sum
The cells which are not hole cells and not interpolated of all the donors’ inverse distances.
are considered as calculated cells.2 For these cells, the
governing equations is solved. Thus, the donor cells are
under the calculated cells category.
The test case
The output information from this process is usually Eglin1 is by far the most dominating test case for code
called the Domain Connectivity Information (DCI).15,18 In validations for the store separation problem. That can be
the OpenFOAM code and literature, the names cellTypes attributed to the availability of the geometry, the experi-
and cell classifications19 have been used as well. mental test reports, and many published numerical results.
The DCI for the pitching airfoil case is presented in Basically, it is a wind tunnel test which was conducted in
Figure 2. For the overset mesh (Figure 2(a)), all the cells Arnold Engineering Development Center (AEDC) in 1990.
are calculated (in blue) except the cells at the outer The test object includes wing, pylon, and finned store
boundary (overset patch). These cells are interpolated (in (Figure 4). The test results are available at two free-stream
yellow) from the background mesh. For this case, there are Mach number (Ma): 0.95 and 1.2, while the angle of attack
no hole cells in the overset mesh. was fixed to 0 o for the two tests. In this work, the transonic
For the background mesh (Figure 2(b)), the hole cells (Ma = 0.95) test case has been used where the store tra-
(in red), which are the projection of the airfoil, are ex- jectory data is available from time = 0.0 s to time = 0.33 s.
cluded from the computation. The neighbor cells of the
hole cells are considered as interpolated cells. The re-
maining cells are the calculated cells.
Geometry
Once the DCI are computed, the donor cells (e.g., The full-scale geometry as per the test report1 has been
D0  D4 in Figure 3) for each interpolated cell (e.g., I0 in created using OpenVSP3 (for the airfoils) and Salome.4
Figure 3) are known. This is often called the interpolation Both codes are free and open-source. The combined wing,
stencil. pylon, and store geometry are shown in Figure 4.
4 Proc IMechE Part G: J Aerospace Engineering 0(0)

Figure 4. Eglin test case full-scale geometry.

The proposed strategy Figure 5. Block diagram showing the adapted foam-extend
overset library as a new “cellCellStencil” inside OpenFOAM in
The preliminary results of this work, which have been
blue color.
presented in Ref. 20, indicate clearly that the standard
OpenFOAM implementation has two main limitations:
The implementation in this way allows the selection of
(i) The flow solver: More accurate and stable com- the adapted OGA at run time, without changing any line of
pressible solver is required. code in the standard OpenFOAM. The developed interface
(ii) The OGA algorithm: A more efficient, faster and class contains “oversetMesh” object from the ported
robust algorithm is required that can classify the cells overset library as shown in Listing 1. After reading all
accurately for complex cases. Having an accurate necessary data from the case folder (e.g., fringes selecting
classification/DCI would improve the accuracy of the settings, hole patches, and the cells in each overlapped
results as well. mesh), the object is initialized (Listing 2). Any Open-
FOAM solver, which supports the dynamic meshes, calls
More details on the observed limitations may be seen in the “update()” function at each time step; thus, the OGA-
Appendix II. related part is written inside the inherited “update()”
function. The “oversetMesh” object provides the required
The flow solver. HiSA5 is an external (non-standard DCI information for each mesh (Listing 3). Finally, the
OpenFOAM) open-source and free solver which utilizes interpolation stencil depicted in Figure 3 and the inter-
the OpenFOAM libraries. It is a density-based and cou- polation weights from each donor (according to the se-
pled solver which can solve both transient and steady-state lected overset interpolation scheme) are returned to the
cases. It has been developed at the Council for Scientific flow solver (Listing 4).
and Industrial Research, South Africa (CSIR).6 HiSA Listing 1: Declaration of the adapted overset object
solver has been verified and used here together with the (foamExtendStencil.H).
overset method. The code is also parallelizable through
OpenFOAM parallel library.
Similar to most density-based coupled solvers, it solves
the following coupled vector form of conservation of
mass, momentum, and energy equations: Listing 2: Initialization of the adapted overset object
∂W (foamExtendStencil.C).
þ =  FðW Þ ¼ QðW Þ (8)
∂t
where W is the conservative variables vector, F(W) is the
flux vector, and Q(W) is the source terms vector. More
Listing 3: Sample code for obtaining the DCI infor-
information about the theory, code, and many validation
mation from the adapted overset object
cases can be found in Refs. 21 and 22.
(foamExtendStencil.C).
The adapted OGA algorithm. Under the OpenFOAM
platform, there is an alternative open-source overset li-
brary in the foam-extend fork.7 Its algorithm is more
robust and faster than the standard OpenFOAM overset
library. The first step was to test the capability of this
library. After that, the library was ported to OpenFOAM,
since many data structures were not compatible between
the two forks of OpenFOAM. The third step was to im-
plement this library as a new “cellCellStencil” inside the
openFOAM overset library as shown in the block diagram
in Figure 5. For this purpose, the interface class “foa-
mExtendStencil” has been developed.
Abuhanieh et al. 5

Listing 4: Sample code for transferring the interpolation function is called for each acceptor. Thus, the foam-
weights to the flow solver (foamExtendStencil.C). extend interpolation methodology, where one func-
tion call calculates the weights for all acceptors, was
not used.
(iii) In some cases, the overset is used only to assemble
the computational domain from different meshes,
and no motion is involved. In this case, running the
OGA algorithm every time step is unnecessary.
Listing 5: Sample code for using the “globalCellCells” There is an option in the foam-extend library to
function (foamExtendStencil.C). prevent the assembly again after the first time
(“cacheFringe”); however, it works only with
a specific fringe type (“overlap”). This option was
not used in this work; instead, this functionality was
implemented at a higher level inside the “foa-
mExtendStencil” class. In this way, the overall
OGA algorithm works only at the first time step,
while the interpolation (computationally much
cheaper than the OGA) is done at each time step as
usual. Furthermore, this functionality has been
extended by introducing a user-defined variable
A shortcoming has been discovered during the testing (“oversetFrequency”). This extension can be used
is that foam-extend library was not able to collect all the to reduce the overhead for the transient cases too
donors from all the processors in case of parallel run (step by executing the OGA every N time steps instead
(iii) in The overset/chimera method section). As a conse- of each time step. Lijewski23 compared the
quence, the results were changing according to the number overhead and the accuracy of assembling the OGA
of the used cores. every 2 and 10 time steps. According to his
This problem has been resolved by getting only the conclusion, the reduction in accuracy was small.
main donors (D0 in Figure 3) from the “oversetMesh”
object. Then, for collecting the remaining/extended do-
nors (D1  D4 in Figure 3) for each main donor, the Comparison with the standard OpenFOAM
“globalCellCells” function from the main class “cell- After checking all the standard solvers and the dynamic
CellStencil” have been utilized properly (Listing 5). This mesh libraries in OpenFOAM, the only available option
function uses the faces of the cell to find its neighbors to solve the store separation problem has been by using
(face-walk). Furthermore, it implements all the necessary the overRhoPimpleDyMFoam solver. This solver is
synchronization tasks between the processors in case of a segregated, pressure-based compressible solver which
parallel executions. supports the overset/chimera dynamic mesh library.
The final step was to integrate the new library with the
HiSA solver. During the early stages of the HiSA code
Mesh. To create the mesh, different open-source tools
investigation, it was observed by the authors that the solver
have been used, including Tetgen,24 SUMO,25 cfMesh,8
in principle shall be able to run any dynamic mesh library
and snappyHexMesh.9 After several trials, snappyHexMesh
from the standard OpenFOAM platform. Since the standard
(unstructured, octree-based, and hexa-dominant) has
OpenFOAM overset class (“dynamicOversetFvMesh”) is
been selected for mesh generation. Compared to the other
derived from the dynamic mesh class, and it uses the
meshing tools mentioned above, it retrieves the surface
“cellCellStencil” objects for performing the OGA op-
mesh with good quality while keeping the minimum cell
erations only, a proper implementation for the new OGA
volume relatively high and the number of cells relatively
was sufficient. Thus, the HiSA solver works with new
low. The minimum cell size has a direct effect on the
library without any modification on its source code. That
computation time of the overRhoPimpleDyMFoam solver
was an advantage of using the standard OpenFOAM as
since it is a segregated solver. Although it is an implicit
the common development platform.
solver, the maximum Courant number, which can guarantee
The differences between the adapted OGA library
a stable run is still limited because of the complexity of the
which was developed in this work and the original foam-
overset scheme. Increasing the minimum cell size reduces
extend library can be summarized as follows:
the maximum Courant number, which guarantees a more
stable run. The Courant number defined in OpenFOAM is
(i) The functionality for finding the extended donors
expressed as:
cells is different as explained previously.
(ii) In order to maintain the compatibility with the
1 X 
standard OpenFOAM, the overset interpolation Co ¼ Δtλ; λ¼ χ  (9)
functionality was implemented in a way that the 2V faces facei
6 Proc IMechE Part G: J Aerospace Engineering 0(0)

Figure 7. The residuals convergence for the developed

Figure 6. A slice from the 0.58 M cells mesh. The background
strategy and the standard OpenFOAM.
mesh is in blue and the overset mesh is in red.

where Co is the Courant number, Δt is the time step size, λ

is the characteristic time scale, V is the cell volume, and χ
is the volumetric flux at each face. The maximum Courant
number is the highest value of Co among all cells.
For comparison purposes with the developed strategy,
and considering the high computational time of the
overRhoPimpleDyMFoam solver, including the overset
overheads (see Appendix II), a very coarse mesh with
0.58 M cells has been generated and tested. A slice from
this mesh is shown in Figure 6. Figure 8. Pressure coefficient versus X/L for store at time =
0.0 s and roll angle = 5o for the standard OpenFOAM and the
Case setup and the standard OGA algorithm. For the sake of developed strategy using the 0.58 M cells mesh.
conciseness, the case setup and the details of the used
OGA algorithm are included in Appendix II. The case
parameters such as mass, center of mass, and moment of
inertia have been taken from the test report.1

Results
Steady-state studies
First, the steady-state Euler solution has been obtained.
Since the overset in OpenFOAM can be used only with
transient cases, a transient simulation (without motion) is
driven to steady-state by running the case until the re- Figure 9. Mach number plot at time = 0.17 s (plane: y normal
siduals reach convergence. The residuals are plotted in at the initial center of mass) for the standard OpenFOAM.
Figure 7, where the density-based HiSA solver used with
the developed strategy yields the density residuals, while
the pressure-based solver overRhoPimpleDyMFoam used for the (Cp) at this region in Figure 8. This example il-
with the standard OpenFOAM yields the pressure lustrates the direct effect of the OGA accuracy on the
residuals. results’ accuracy.
The pressure coefficient (Cp) at roll angle 5° is plotted
in Figure 8 for the standard OpenFOAM (over-
RhoPimpleDyMFoam plus the cellVolumeWeight OGA)
Transient simulations
and for the developed strategy (HiSA plus the adapted After initializing the field variables by the obtained steady-
OGA) compared with the experiment. Although the mesh state solution, the actual transient run started by applying
is coarse, the developed strategy solution is able to match 6-DOF motion solver. The Mach number at time = 0.17 s
the Cp trend of the experiment. However, for the standard over a plane is plotted in Figure 9 for the standard
OpenFOAM solution, there is a significant deviation OpenFOAM. A diffused solution can be observed in the
specially in the middle section. That can be attributed to store mesh near the outer patch (the outer boundaries of
the inaccuracy of the OGA algorithm, since the 5° roll the red box in Figure 6). This can be attributed to the first-
angle resides in the small region between the store and the order overset interpolation scheme which has been used in
pylon. Examining Figure 33(a) shows that this region was this case.
not solved by the flow solver; the cells there are either Unlike the standard OpenFOAM results, with the
interpolated or holes. Thus, the value of the fields there developed strategy shown in Figure 10, the outer patch of
take the free-stream values. That explains the value of zero the store cannot be distinguished, primarily due to the
Abuhanieh et al. 7

Figure 13. Linear displacement in the z-direction for the

Figure 10. Mach number plot at time = 0.17 s (plane: y store with comparison to the experiment for the standard
normal at the initial center of mass) for the developed strategy. OpenFOAM and the developed strategy.

Figure 11. Linear displacement in the x-direction for the Figure 14. Roll angle of the store with comparison to the
store with comparison to the experiment for the standard experiment for the standard OpenFOAM and the developed
OpenFOAM and the developed strategy. strategy.

Figure 12. Linear displacement in the y-direction for the Figure 15. Pitch angle of the store with comparison to the
store with comparison to the experiment for the standard experiment for the standard OpenFOAM and the developed
OpenFOAM and the developed strategy. strategy.

higher order interpolation scheme used (inverse distance).

Furthermore, the shocks are more observable, thanks to
the used second-order upwind scheme (AUSM+-up).26 A
brief description of this scheme is presented in Appendix III.
The linear displacements of the store versus time are
plotted in Figures 11–13. Compared with the experiment,
both results are in good agreement. However, the results of
the developed version of the code are better in all three
directions.
The angular displacements of the store versus time are
Figure 16. Yaw angle of the store with comparison to the
plotted in Figures 14–16, where a general agreement with
experiment for the standard OpenFOAM and the developed
the experiment can be observed for the two results. The
strategy.
developed strategy results are better for pitch and yaw
angles. However, for roll angle, the standard OpenFOAM
shows a better agreement with the experiment. That may That makes the roll angle too sensitive to the errors in
be attributed to the smaller value of the moment of inertia aerodynamic force calculations. Furthermore, with the
for rolling, compared to the pitching and yawing (27 kg.m2 absence of results with finer meshes for the standard
versus 488 kg.m2 for the Eglin case store at full-scale1). OpenFOAM, it is hard to evaluate the reliability of this
8 Proc IMechE Part G: J Aerospace Engineering 0(0)

Table 1. Timing comparison between the standard Results

OpenFOAM cellVolumeWeight OGA algorithm and the adapted
foam-extend one. Steady-state simulations
The standard The adapted First, the steady-state (Euler) solution has been obtained.
OpenFOAM foam-extend Standard/adapted The pressure coefficient (Cp) at roll angle 5° is plotted in
time ratio Figure 18, and compared with the experiment, the
# Of cores Time [s) # Of cores Time [s) [-] starCCM code results in Ref. 10, and the Fluent code
1 48.0 1 5.75 8.35 results in Ref. 11.
4 33.0 4 4.54 7.27 It can be seen that the developed strategy results for Cp
8 21.0 8 3.02 6.95 agree quite well with the experimental results and with the
two referenced numerical solutions. Apparently, the two
16 19.0 16 2.19 8.68
shocks at X/L; 0.27 and at X/L; 0.67 have been captured
well. The three Euler numerical results (and many others,
for example, Refs. 6, 8, and 9) show deviation in the small
region (X/L; 0.87  X/L; 0.94) or the last two measurement
particular behavior. Unlike the linear displacements, the points in the experiment report. This deviation can be
accuracy of the angular displacements depends on the attributed to the viscosity effect (see The viscosity effect),
accuracy of the calculated moments from CFD, which in since in literature, for example, Refs. 6 and 27, the Navier–
turn depends significantly on finding the accurate shock Stokes equations have been solved and there was an
location. That is why obtaining accurate results for angular improvement in that particular region for the same Mach
displacements is more challenging. number and roll angle. This roll angle (f = 5o) is of special
interest, since it is located in the tiny gap between the store
Code performance and pylon.28 Unlike the viscosity effect, adding the sting
(shown in Figure 19), which was gaged to the store to
The developed strategy improves the speed of solving the measure the pitching and yawing moments1 during the
store separation problems in two ways: testing, apparently did not improve the results
appreciably.8,9
(i) The adapted OGA algorithm is faster. Table 1 shows There is a third shock which appears at the tail of the
a comparison between the two algorithms (the store (at X/L; 0.98), which can be observed clearly in
standard OpenFOAM cellVolumeWeight and the Figure 20; however, there is no measurement point at this
adapted one) using the same node on TRUBA HPC10 location in the experimental data.
(mesh size 0.58 M) for one time step. It is clear that
the improved algorithm is much faster for this mesh
configuration (at least 7X). However, the effective
speedup with larger size meshes may need further
investigation.

(ii) The used flow solver is coupled and is able to solve the
problem with a higher Courant numbers, which allows
obtaining solutions with larger time steps.

Accuracy studies for the developed strategy. In this section,

the accuracy of the developed version of OpenFOAM and
its mesh independency have been studied by considering Figure 17. Slice from the fine mesh (4.48 M cells), the
different levels of mesh refinements and comparing the background mesh in blue and the overset mesh in red.
results with the Eglin test case experiments. The results are
also compared with the solutions of two popular com-
mercial CFD codes, starCCM,10 and Fluent,11 that have
also used an overset algorithm for this test case.
Four different meshes (coarse with 1.04 M cells, me-
dium with 1.62 M cells, fine with 4.48 M cells, and very-
fine with 6.20 M cells) have been generated for this study
using snappyHexMesh. A slice from the fine mesh (at
center of mass and y normal) is shown in Figure 17.

Case setup and the OGA algorithm results

Figure 18. Pressure coefficient versus X/L for store at time =
For the sake of conciseness, the case setup details and the 0.0 s and roll angle = 5o for the developed strategy (fine mesh),
used OGA algorithm results are included in Appendix II. starCCM, and Fluent.
Abuhanieh et al. 9

Figure 19. The used sting to measure the pitching and

yawing moments attached to the tail of the tested/metric store.1
Dimensions are in inches (1/20 scale).

Figure 21. Roll angle of the store. Comparison between the

Euler solutions for the four different meshes.

Figure 22. Pitch angle of the store. Comparison between the

Figure 20. Mach number contour at the lower side of the
Euler solutions for the four different meshes.
wing and store at time = 0.0 s for the fine mesh.

For the small oscillations observed in the present re-

sults, the authors have noticed similar oscillations in other
cases too while using snappyHexMesh for mesh gener-
ation. The octree concept used by snappyHexMesh does
not preserve the input surface mesh; instead, it snaps the
surface of the cells after cutting them. That can be the
reason for such oscillations which may be considered as
a minor post-processing issue.
In Figure 20, the Mach number contour for steady case Figure 23. Yaw angle of the store. Comparison between the
is plotted at the lower side of the wing showing the store as Euler solutions for the four different meshes.
well. It may be observed that multiple shocks are ap-
pearing there and they are not symmetrical around the
store. From this view (normal to z-direction), it can be
observed that the un-symmetry of the shocks around the
store is the main reason of yawing the store during
separation.

Transient simulations
The next step, is to start the transient simulation by
applying the 6-DOF motion solver. The effect of the
mesh refinement level on the solution has been studied
here. The angular displacements obtained by solving the Figure 24. Store location with time. The Mach number
Euler equations for the four meshes (coarse, medium, contours represent the fine mesh results and the experimental
fine, and very-fine) are compared in Figures 21–23, results are in blue.
respectively. As may be observed, in general, the solution
is converging monotonically toward the experimental
results as the mesh is refined. That indicates that the
discretization error is decreasing. Additionally, minor The change of store location with time is shown in
differences between the solutions of the fine and very- Figure 24. The blue color represents the experimental
fine meshes indicate that the mesh convergence of the results and the Mach number contours represent the ob-
solution is also reached. tained CFD results for the fine mesh. Minor deviations
10 Proc IMechE Part G: J Aerospace Engineering 0(0)

Figure 26. Roll angle of the store with comparison to the

Figure 25. The calculated center of pressure (the red sphere) experiment and two commercial codes for the fine mesh.
and the center of mass (the black sphere) at time = 0.05 s. The
contour plot shows the static pressure on the store.

from the experiment can be noticed, which reflects that

accurate displacements are obtained.
To highlight the importance of the applied external
forces to the store by the ejectors during separation, the
Center of Pressure (CoP) has been calculated at time =
0.05 s using the following relation:
Z
Figure 27. Pitch angle of the store with comparison to the
xp experiment and two commercial codes for the fine mesh.
CoP ¼ Z (10)
p

where CoP is the center of pressure coordinate vector, x is

the point coordinate vector, and p is the static pressure.
Figure 25 shows the location of the calculated CoP (the
red sphere). Since the center of mass for the store (the
black sphere) is located after the CoP and the lift force
direction is upward at this time step, the store is expected
to pitch nose-down. However, as per the results (e.g.,
Figure 22), the store is pitching nose-up. That is simply
due to the external forces, and the same can be verified Figure 28. Yaw angle of the store with comparison to the
easily by calculating the total moments on the store with experiment and two commercial codes for the fine mesh.
and without the ejector forces. To incorporate the external
forces by the ejectors, a new force restrained model
(“ejectorExternalForce”) has been developed in this work.
This model applies a fixed force on a specific point which
moves with the body up to a certain length (the stroke
length specified in the experiment).
Since the angular displacements are normally con-
sidered the critical results for this kind of analysis, the
developed strategy results for the fine mesh are compared
to the experiment, the starCCM code results in Ref. 10 and
the Fluent code results in Ref. 11 as shown in Figures 26–28,
for roll, pitch and yaw angles, respectively. Figure 29. Pressure coefficient versus X/L for store at time
It may be noticed that for the three angles, the de- = 0.0 s and roll angle = 5o for the fine mesh. Comparison
veloped strategy results are comparable to the results of between the Euler and the Navier–Stokes solutions.
starCCM and Fluent. Both codes used the overset as
a dynamic mesh technique to obtain the results presented
here.
The average y+ (the non-dimensional wall distance) for
The viscosity effect. To evaluate the viscosity effect on this this mesh is around 200; thus, wall functions were used.
particular test case, the Euler solution of the fine mesh Furthermore, the kω  SST (shear stress transport) tur-
(4.48 M) has been compared to the Navier–Stokes solution. bulence model29 has been used.
Abuhanieh et al. 11

expected that solving the Navier–Stokes equations with an

appropriate viscous mesh can provide better results for the
steady-state solution.
For the transient run, the angular displacements results
are shown in Figures 30–32 for roll, pitch, and yaw angles,
respectively. The work of Huang et al.30 can be considered
as a comprehensive study which analyzed the effect of the
mesh (size and topology), the model scale, the wing
symmetry, and the turbulence models. Their best Navier–
Stokes result for angular displacements (using Spalart–
Figure 30. Roll angle of the store. Comparison between the Allmaras turbulence model31) has been included along
Euler and Navier–Stokes solutions for the fine mesh. the work of Sunay et al.6 in these figures and compared
with the present results for both Euler and Navier–Stokes
simulations. Unlike the steady-state simulation, for the
trajectory prediction, it can be observed that solving the
viscous flow even with very well prepared meshes may not
necessarily provide better results than the Euler solution.
A similar observation can be made by reviewing the re-
sults of Figure 19 in Huang et al.30

Conclusions
A new strategy which involves modifying and adapting
Figure 31. Pitch angle of the store. Comparison between the the foam-extend overset library, and integrating it with the
Euler and Navier–Stokes solutions for the fine mesh. density-based coupled solver HiSA for solving the store
separation problem within the OpenFOAM platform has
been proposed and implemented. The new strategy is more
accurate and efficient than the standard OpenFOAM so-
lution for this problem; additionally, it is capable of
solving industry-scale cases/geometries. Furthermore, it is
faster at least by 7 times. A mesh refinement study has
been conducted, and an improvement in the results have
been observed with increasing the number of cells, which
is a good sign about the accuracy of the developed
strategy. The obtained results agree well with the exper-
imental results and are comparable with the results of the
Figure 32. Yaw angle of the store. Comparison between the two well-known commercial codes. A Navier–Stokes
Euler and Navier–Stokes solutions for the fine mesh. solution has been compared to the Euler solution. It is
shown that for the transonic regime, at least for this case,
the Euler solution for predicting the trajectory can provide
For the steady-state solution, Figure 29 shows the satisfactory results, which in practice can save time and
pressure coefficient (Cp) comparison at roll angle 5°. effort. Similar observations have been noticed in the
Compared to the experiment, for the Navier–Stokes so- literature.
lution, the second shock at X/L; 0.67 is slightly shifted to
the right. That can be attributed to the high y+ mesh used. Acknowledgments
A finer mesh in the viscous region can provide a better The numerical calculations reported in this paper were partially
result. With comparison to the Euler solution, the shock performed on the resources of the ULAKBIM High Performance
near the tail of the store at X/L; 0.98 is detected more and Grid Computing Center of The Scientific and Technological
sharply. However, since no measurement points are avail- Research Council of Turkey. We thank the Turkish Aerospace for
able at the test report in this region, improved compar- the support provided for this study.
isons of the Navier–Stokes solution with experiments
cannot be confirmed at this time. Declaration of Conflicting Interests
Unlike the developed strategy results for both Euler and The author(s) declared no potential conflicts of interest with
Navier–Stokes solutions, the results obtained by Sunay respect to the research, authorship, and/or publication of this
et al.6 and Pandya et al.27 matched better with the ex- article.
periment, for the last two measurement points. Further-
more, in their solutions, the second shock at X/L x 0.67 Funding
and the third shock at X/L x 0.98 are weaker, but still The author(s) received no financial support for the research,
respect the experimental results. Consequently, it is authorship, and/or publication of this article.
12 Proc IMechE Part G: J Aerospace Engineering 0(0)

ORCID iD 14. Benek JA, Steger JL and Dougherty FC. A flexible grid
Saleh Abuhanieh  https://orcid.org/0000-0002-3620-8546 embedding technique with applications to the Euler equa-
tions. In: 6th Computational Fluid Dynamics Conference.
MA, USA: Danvers; 1983.
Notes
15. Martin JE, Noack RW and Carrica PM. Overset grid as-
1. www.openfoam.com sembly approach for scalable computational fluid dynamics
2. www.openfoam.com/documentation/guides/latest/doc/ with body motions. J Comput Phys 2019; 390: 297–305.
guide-overset.html 16. Roget B and Sitaraman J. Robust and efficient overset grid
3. www.openvsp.org assembly for partitioned unstructured meshes. J Comput
4. www.salome-platform.org Phys 2014; 260: 1–24.
5. hisa.gitlab.io 17. Jingjing F and Chao Y. Enhancement and application of
6. www.csir.co.za overset grid assembly. Chin J Aeronautics 2010; 23(6):
7. www.sourceforge.net/p/foam-extend 631–638.
8. www.cfmesh.com 18. Noack RW, Boger DA, Kunz RF, et al. Suggar++: An improved
9. www.openfoamwiki.net/index.php/SnappyHexMesh general overset grid assembly capability. In: AIAA Computa-
10. https://www.truba.gov.tr/index.php/en/main-page/ tional Fluid Dynamics. San Antonio, TX, USA; 2009.
19. Chandar D. Development of a parallel overset grid frame-
References work for moving body simulations in OpenFOAM. J Appl
1. Heim ER. CFD Wing/pylon/finned Store Mutual Interference Computer Sci Mathematics 2015; 9(2): 22–30.
Wind Tunnel Experiment. ADB152669 Technical Report 1991. 20. Abuhanieh S, Akay HU and Bicer B. A new strategy for
2. Arnold RJ and Epstein CS. AGARD flight test techniques solving store separation problems using OpenFOAM
series on store separation flight testing. AGAR Dograph (abstract), In: 32nd International Conference on Parallel
1986; 5(300). Computational Fluid Dynamics. Nice: France, 2021.
3. Jamison KA. Grid-mode transonic store separation analysis (virtual).
using modern design of experiments. In: 31st Congress of 21. Heyns JA, Oxtoby OF and Steenkamp A. Modelling high-
the International Council of the Aeronautical Sciences. speed flow using a matrix-free coupled solver. In: 9th
Brazil: Belo Horizonte, 2018. OpenFOAM Workshop. Zagreb, Croatia, 2014.
4. Davids S and Cenko A. Grid based approach to store 22. Abuhanieh S, Bicer B and Sahin M. A validation for the
separation. In: 19th AIAA Applied Aerodynamics Confer- OpenFOAM-Hisa solver for drag prediction (abstract), In:
ence. Anaheim, CA, USA, 2001. 32nd International Conference on Parallel Computational
5. Panagiotopoulos EE and Spyridon DK. CFD transonic store Fluid Dynamics. Nice: France, 2021. (virtual).
separation trajectory predictions with comparison to wind 23. Lijewski L. Comparison of transonic store separation tra-
tunnel investigations. Int J Eng 2010; 3(6): 538–553. jectory predictions using the Pegasus/DXEAGLE and
6. Sunay YE, Gulay E and Akgul A. Numerical simulations of Beggar codes. In: 15th Applied Aerodynamics Conference.
store separation trajectories using the eglin test. Scientific Atlanta, GA, USA, 1997.
Tech Rev 2013; 63(1): 10–16. 24. Si H. Tetgen, a delaunay-based quality tetrahedral mesh
7. Dehghan M, Davari AR and Manshadi MD. Numerical generator. ACM Trans Math Softw 2015; 41(2): 1–36.
investigation on the weight, speed, and installation location 25. Tomac M and Eller D. Towards automated hybrid-prismatic
effects on fuel tank separation trajectory. Proceedings of the mesh generation. Proced Eng 2014; 82: 377–389.
Institution of Mechanical Engineers, J Aerospace Eng 2017; 26. Liou M. A sequel to AUSM, part ii: AUSM+-up for all
231(13): 2331–2344. speeds. J Comput Phys 2006; 214(1): 137+170.
8. Lijewski L and Suhs N. Chimera-eagle store separation. In: 27. Pandya MJ, Frink NT and Noack RW. Progress toward
AIAA Atmospheric Flight Mechanics Conference. SC, overset-grid moving body capability for USM3D un-
USA: Hilton Head Island, 1992. structured flow solver. In: 17th AIAA Computational Fluid
9. Demir H, Alemdaroglu N and Ozveren V. External store Dynamics Conference. Toronto, ON: Canada, 2005.
separation from fighter aircraft. In: RTO AVT Symposium on 28. Snyder DO, Koutsavdis EK and Anttonen JSR. Transonic
Functional and Mechanical Integration of Weapons and store separation using unstructured CFD with dynamic
Land and Air Vehicles. Williamsburg, VA, USA: RTO-MP- meshing. In: 33rd AIAA Fluid Dynamics Conference and
AVT-108, 2004. Exhibit. Orlando, FL, USA, 2003.
10. MacLucasa D and Gledhillb I. Time-accurate transonic CFD 29. Menter F, Kuntz M and Langtry R. Ten years of industrial
simulation of a generic store release case. R D J South Afr experience with the sst turbulence model. In: Proceedings of
Inst Mech Eng 2018; 34: 9–16. the fourth international symposium on turbulence, heat and
11. Khaware A, Sivanandham A and Gupta VK. Numerical mass transfer. Turkey: Antalya, 2003.
simulation of store separation trajectory for Eglin test case 30. Huang H, Blyth RH, Prior MA, et al. A comparison of
using overset mesh. In: AIAA SciTech Forum. FL, USA: approaches to multi-body relative motion using the Kestrel
AIAA Aerospace Sciences MeetingKissimmee; 2018. solver. In: AIAA Scitech 2019 Forum. San Diego, CA, USA,
12. Wadibhasme R. Exploration and Implementation of Various 2019.
Dynamicmesh in OpenFOAM (Master thesis). Centre for 31. Spalart P and Allmaras S. 30th Aerospace Sciences Meeting
Modeling and Simulation, Savitribai Phule Pune University; and Exhibit. RenoUSA: NV, 1992.A one-equation turbu-
2016. lence model for aerodynamic flows
13. Menon S. A numerical study of droplet formation and be- 32. Liou M and Steffen C. A new flux splitting scheme.
havior using interface tracking methods interface tracking J Comput Phys 1993; 107: 23–39.
method (PhD thesis). Mechanical and Industrial Engi- 33. Liou M. A sequel to AUSM: AUSM+. J Comput Phys 1996;
neering. University of Massachusetts Amhers; 2011. 129(2): 364–382.
Abuhanieh et al. 13

Appendix I Table 2. Solver settings.

Solver overRhoPimpleDyMFoam
Notations OpenFOAM version v2006
Co. Courant number [-] Flow Transient inviscid
CoP center of pressure coordinate [m] Discretization schemes First-order upwinding in space
d distance vector magnitude [m] and first-order implicit in time
Di donor cell i [-] (backward euler)
e specific internal energy [J/kg] Time step 1 × 104 s (maximum)
F applied forces vector [N] Maximum courant number 5.0
F(W) flux vector (convective flux and viscous flux) Pseudo time (inner loop) Maximum iterations = 30,
H angular momentum vector [kg  m2/s] settings tolerance = 1 × 104, there
Ii interpolated cell i [-] is no dual time scheme option
M applied moments vector [N  m] Free-stream Mach number 0.95
mV linear momentum vector [kg  m/s] Overset interpolation cellVolumeWeight (first-order)
ND number of donors [-] scheme
p static pressure [Pa]
q heat flux vector [W/m2]
Q(W) source terms vector
t time [s] Table 3. Boundary conditions.
v fluid velocity vector [m/s]
V cell volume [m3] Patch/field Pressure Velocity Temperature
W conservative fluid flow variables vector Farfield Free-stream Free-stream inletOutlet
x coordinate vector [m] pressure velocity
y+ non-dimensional wall distance [-] wing_and_pylon zeroGradient Slip zeroGradient
θ pitch angle (around y) [o] Store zeroGradient movingWall zeroGradient
λ the characteristic time scale [1/s] velocity
ρ density [kg/m3] oversetPatch Overset Overset Overset
σ mechanical stress tensor [N/m2] Symm Symmetry Symmetry Symmetry
τ viscous stress tensor [N/m2]
f roll angle (around x) [o]
χ face volumetric flux [m3/s]
ψ yaw angle (around z) [o] Table 4. dynamicMeshDict settings.
ξi field value at cell i
ωi interpolation weight of cell i [-] dynamicFvMesh dynamicOversetFvMesh
Δt time step size [s]. Solver sixDoFRigidBodyMotion

Appendix II
(A) The solution obtained by the overRhoPimpleDyMFoam
overRhoPimpleDyMFoam case setup is reasonable. However, the accuracy is not good
The details of the setup for the case solved with over- enough.
RhoPimpleDyMFoam are summarized in Tables 2–4. (B) A converged solution with a second-order upwind-
ing scheme was not obtained. This can be attributed
to the limited numerical stability of the used seg-
The standard OGA regated flow solver.
For this test case (Eglin 0.58 M cells) only the first-order (C) The computational time is high (842 core-hours for
cellVolumeWeight scheme worked with the over- the Eglin case with 0.58M cells), and more than the
RhoPimpleDyMFoam solver. In the OpenFOAM im- half of each time step execution time is consumed
plementation, the cellVolumeWeight is not only an by the overset processes. Thus, running a fine mesh
interpolation scheme, to be used by the interpolated cell to may easily turn the case to be non-practical to be
get the value from the donors cells, but it includes the solved.
OGA process, which has been described in The overset/ (D) Only the first-order overset interpolation scheme was
chimera method section. A slice showing the overset DCI used (cellVolumeWeight), since the other schemes
for the overset/store mesh (Figure 33(a)) and the back- (inverseDistance and leastSquares) OGA algorithms
ground mesh (Figure 33(b)) is shown in Figure 33. failed to classify the cells (either calculated, hole or
In this work, studying the current/available capabilities interpolated) correctly. A first-order overset in-
of the standard OpenFOAM to solve the store separation terpolation scheme is not enough to obtain accurate
problem reveals the following limitations: results.
14 Proc IMechE Part G: J Aerospace Engineering 0(0)

Figure 33. The overset DCI for the overset/store mesh (a) and the background mesh (b). The calculated cells are in blue, the
interpolated cells are in yellow and the hole cells are in red (plane: y normal at the initial center of mass). The observed limitations for the
standard OpenFOAM.

The developed strategy case setup. The adapted OGA results

The details of the setup for the case solved with the
HiSA solver and the developed strategy are summarized in For the adapted overset library, only one OGA algorithm is
Tables 5 and 6. available. However, for selecting the fringe (interpolated
cells) from each mesh, different methods can be selected
(manual, overlap, faceCells, etc.).
Unlike the standard OpenFOAM overset library, in this
implementation, the OGA algorithm is independent from the
interpolation scheme. In the present study, the interpolation
Table 5. Solver settings. schemes which have been implemented in the overset library
Solver HiSA (version 1.4.2) are: inverseDistance and the leastSquares as shown in Figure 5
in blue under the interface class. The prefix “foamExtend” has
OpenFOAM version v2006 been added to distinguish them from the ones used by the
Flow Transient inviscid standard library to be able to select them during run time. A
Discretization schemes Second-order upwinding in slice showing the DCI for the overset/store mesh (Figure
space (AUSM+-up) and 34(a)) and the background mesh (Figure 34(b)) is shown in
first-order implicit in time Figure 34. The overlap fringe has been used here.
(backward euler)
Time step 1 × 104 s (maximum)
Maximum courant number 60.0
Appendix III
Pseudo time (inner loop) Maximum iterations = 100,
settings tolerance = 5 × 103, AUSM+-up Scheme
using dual time scheme
AUSM+-up is a second-order upwinding scheme for shock
Free-stream Mach number 0.95
capturing, which is an extension to the AUSM (Advection
Overset interpolation scheme foamExtend_inverseDistance
Upstream Splitting Method) family of schemes.32,33 It has
(first-second-order)
been developed to improve the prediction of flows at low
dynamicMeshDict settings Same as the
Mach speeds, making the original AUSM scheme work ef-
overRhoPimpleDyMFoam
settings
ficiently at all speeds. The algorithm can be described as
follows26:

Table 6. Boundary conditions.

Patch/field Pressure Velocity Temperature

Farfield Characteristic farfield pressure Characteristic farfield velocity Characteristic farfield temperature
wing_and_pylon Characteristic WallPressure Slip Characteristic WallTemperature
Store Characteristic WallPressure movingWall velocity Characteristic WallTemperature
oversetPatch Overset Overset Overset
Symm Symmetry Symmetry Symmetry
 Abuhanieh et al. 15

Figure 34. The overset DCI for the overset/store mesh (a) and the background mesh (b) for the fine mesh. The calculated cells are in
blue, the interpolated cells are in yellow, and the hole cells are in red (plane: y normal at the initial center of mass).

(i) For each interface, calculate the Mach number’s left M1=2 ¼ Mþ 
ð4Þ ðML Þ þ Mð4Þ ðMR Þ
and right states
 (15)
uL=R Kp 2 pR  pL
ML=R ¼ (11)  max 1  σM ,0
a1=2 fa ρ1=2 a21=2

where ML/R are the Mach number left and right states, uL/ where Kp and σ are constants and pR and pL are the
R are the convective velocities (v n), and a1/2 is the speed pressure right and left states, respectively. ρ1/2 is the av-
of sound at the interface (can be obtained by averaging the erage of the density left and right states.
aL and aR).
(vii) Calculate the mass flux at the interface
(ii) Calculate the mean Mach number
ρL if M1=2 > 0
u2 þ u2 m_ 1=2 ¼ a1=2 M1=2 (16)
2 ρR otherwise
M ¼ L 2 R (12)
2a1=2
(viii) Calculate the pressure flux
(iii) Calculate the reference Mach number
p1=2 ¼ P þ 
ð5Þ ðML ÞpL þ P ð5Þ ðMR ÞpR
 
2  
Mo2 ¼ min 1,max M ,M∞2 2 ½0; 1 (13) ku P þ 
ð5Þ ðML ÞP ð5Þ ðMR ÞðρL þ ρR Þ fa a1=2 ðuL  uR Þ

(17)
where M∞ is the free-stream Mach number.
where Ku is constant. P ±ð5Þ is a fifth-order polynomial.
(iv) Find the scaling function value (ix) Finally, evaluate the total flux at the interface

fo ¼ Mo ð2  Mo Þ 2 ½0; 1 (14) ψL if m_ 1=2 > 0,

f 1=2 ¼ m_ 1=2 þ p1=2 (18)
ψR otherwise
±
(v) Evaluate the split Mach numbers (M (4)) for the left
Mþ 
ð4Þ ðML Þ and right Mð4Þ ðMR Þ Mach number states.
where ψ is a vector quantity convected by m._ In 1-D, ψ =
The used split Mach number is a fourth-order (1,u,H)T, where H is the total enthalpy. p is the pressure
polynomial function. flux vector, which contains only one pressure term, p = (0,
(vi) Calculate the interface Mach number p,0)T in 1-D.

View publication stats