Interior ballistics two-phase reactive flow model applied to large caliber guided projectile-gun system

11 Pages·2.24 MB·English

by  Mahmoud Rashad, XiaoBing Zhang, Hazem El Sadek.

Hindawi Publishing Corporation
Advances in Mechanical Engineering
Volume 2014, Article ID 698032, 10 pages
http://dx.doi.org/10.1155/2014/698032
Research Article
Interior Ballistics Two-Phase Reactive Flow Model Applied to
Large Caliber Guided Projectile-Gun System
MahmoudRashad,XiaoBingZhang,andHazemElSadek
SchoolofEnergyandPowerEngineering,NanjingUniversityofScienceandTechnology,Nanjing210094,China
CorrespondenceshouldbeaddressedtoMahmoudRashad;[email protected]
Received11December2013;Accepted7February2014;Published29April2014
AcademicEditor:JiinY.Jang
Copyright©2014MahmoudRashadetal. This is an open access article distributed under the Creative Commons Attribution
License,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperly
cited.
Transientcomplexphenomenatakeplaceinagunduringinteriorballisticcycle.Understandingthesephenomenaclearlyand
describingthemathematical modelsaccurately arecrucialtopredictthebehavior ofgunsystemconsideringfiringsafetyand
performance.AmathematicalmodelbasedonEulerian-Eulerianapproachforreactivegas-solidflowarisingduringinteriorballistic
cycleinsidelargecalibernavalgunguidedprojectilesystemwasdeveloped.Themodelincludedthegoverningequationsofmass,
momentum,andenergyforbothphasesaswellastheconstitutivelaws.Thesystemofequationswassolvedusingsecond-order
accurateMacCormacktechnique.One-dimensionalshocktubemodelwasutilizedtotesttheabilityofthenumericalalgorithm
insolvingtheinitialboundaryvalueproblemforthesystemofequationswithshockwavebehavior.Thenumericalmethodwas
verifiedusingexactsolutionofatestproblem.Themovingcontrolvolumeconservationmethodwasusedtohandlethemoving
boundaryandaself-adaptingmethodwasusedtoexpandthecomputationaldomaininordertofollowtheprojectilemotion.
Numericalresultswerevalidatedwithexperimentaldata.Theinteriorballisticsperformanceofa130mmnavalgunwasclosely
predictedusingthepresentedmodelandthenumericalcode.
1.Introduction with minimal computer resources, they cannot address the
physical hydrodynamics of the problem as manifested in
Increasing demands for superior ballistic performance of ignition-induced pressure waves [6]. Continually incidence
large caliber naval guns and up-to-date guided munitions of catastrophic failures in conventionalguns was ultimately
developmentsaretheleadingstimulustogothroughrealistic traced to ignition and combustion instability problems [7–
investigations for the complex phenomenon that occurs 9]. Great demands to conquer and analyze these problems
during interior ballistic cycle. Events inside gun chamber motivatedrevolutionarydevelopmentoftwo-phaseflowinte-
possessatransientbehaviorcoupledwithcomplexreactions rior ballistic models which have the advantage of assessing
underhighpressureandtemperatureconditions.Atthesame theignitionstimulusonflamespreadingandtheformation
time guided munitions must survive against the rigorous of pressure waves. Two-phase flow problems are solved
loads during gun firing and deliver their intended utility at throughtwoapproaches:theEulerian-Eulerianapproachand
thetarget. Eulerian-Lagrangian approach. Models based on Eulerian-
Earlyinteriorballisticsmodels,whichappearedinabun- Eulerian approach [10–14], also called two-fluid models,
dance as computer codes, are of the zero-dimensional consider the gas and solid phase as two continuum flow
(lumped parameters) variety [1]. These models exemplify phases, with each having its mass, momentum, and energy
such assumptions as uniform and instantaneousignition of equations, respectively. These models have the advantages
the entire propellant charge, with combustion taking place in such cases where the density of solid particles possesses
inasmoothlyvarying,well-stirredmixture[2–5].Although high values and the solid phase volume fraction could be
lumpedparametersmodelshavetheadvantagesofsimplicity aneffectiveparameter.So,suchmodelsareappropriateand
aswellasitguaranteeparametricanalysisandoptimization popular in most interior ballistic numerical simulations.
2 AdvancesinMechanicalEngineering
phases and the energy conservation law for the gas phase.
Since the solid propellant ignition is depending on surface
temperature and the model assumes that propellant grains
are incompressible, the energy equation for the solid phase
wasreplacedbytheinterphaseheattransferequationsforthe
granularpropellants.Thestatementofaconstantdensityfor
the solid propellant particles serves as the equation of state
forthesolidparticles;thatis,𝜌𝑝 = constant.Thegoverning
1 2 3 4 5
equationsforvariousphasesarepresentedasfollows.
(1) Breech block (4) Guided projectile Massconservationequationofthegasphaseis
(2) MP propellant grains (5) Barrel
(3) Igniter 𝜕(𝜑𝜌 𝐴)
𝑔 +∇⋅(𝜑𝜌 𝑢 𝐴)=𝑚̇ 𝐴+∑𝑚̇ 𝐴. (1)
Figure1:Schematicillustrationofpropellantchargeinlargecaliber 𝜕𝑡 𝑔 𝑔 𝑐 ign
NGGPS.
Massconservationequationofthesolidphaseis
On the other hand, models based on Eulerian-Lagrangian 𝜕[(1−𝜑)𝜌 𝐴]
𝑝 +∇⋅[(1−𝜑)𝜌 𝑢 𝐴]=−𝑚̇ 𝐴. (2)
approachtreatthegasphaseascontinuum,whichismodeled 𝜕𝑡 𝑝 𝑝 𝑐
by the local averaged equations at the macroscopic scale,
while the solid phase is treated as an individual particle Momentumconservationequationofgasphaseis
whichcanbetracedinspaceandtime[15–20].Thisapproach
hasadvantagesthatitcanhandlepolydispersedparticlesize 𝜕(𝜑𝜌 𝑢 𝐴)
𝑔 𝑔
+∇⋅(𝜑𝜌 𝑢 𝑢 𝐴)
distributionsaswellaspredictingtwo-phaseflowsincluding 𝜕𝑡 𝑔 𝑔 𝑔
(3)
particles that have large accelerations. But such models are
complex and need huge computer resources as well as long =−𝑓𝐴+𝑚̇ 𝑢 𝐴+∑𝑚̇ 𝑢 𝐴−(𝜑𝐴)∇𝑝.
𝑠 𝑐 𝑝 ign ign
timeduringruncomparedtoEulerian-Eulerianmodels.
In this work, a mathematical model based on Eulerian- Momentumconservationequationofsolidphaseis
Eulerian approach for a gas-solid flow of solid granular
propellant and its products of combustion inside the gun 𝜕[(1−𝜑)𝜌𝑝𝑢𝑝𝐴]
barrel during interior ballistics cycle in large caliber naval 𝜕𝑡
gunguidedprojectilesystem(NGGPS)wasestablished.The
(4)
model includes the balance equationsof mass, momentum, +∇⋅[(1−𝜑)𝜌 𝑢 𝑢 𝐴]+∇[(1−𝜑)𝑅 𝐴]
𝑝 𝑝 𝑝 𝑝
andenergyforbothphasesaswellasnecessaryconstitutive
=𝑓𝐴−𝑚̇ 𝑢 𝐴−(1−𝜑)𝐴∇𝑝.
laws. A second-order accurate MacCormack technique was 𝑠 𝑐 𝑝
utilizedtosolvethehyperbolicsystem ofequations[21].G.
A. Sod shock tube test was utilized to verify the ability of Energyconservationequationofthegasphaseis
MacCormack algorithm to solve the initial boundary value
𝜕[𝜑𝜌 𝐴(𝑒 +𝑢 ⋅𝑢 /2)]
problem (IBVP) for the system of equations with shock 𝑔 𝑔 𝑔 𝑔
wavebehavior[22].Themovingcontrolvolumeconservation 𝜕𝑡
(MCVC) method was established to handle the moving
𝑝 𝑢 ⋅𝑢
boundaryandaself-adaptingmethodwasusedtoexpandthe +∇⋅[𝜑𝜌 𝑢 𝐴(𝑒 + + 𝑔 𝑔)]
computational domain in order to follow the movement of 𝑔 𝑔 𝑔 𝜌𝑔 2
theprojectiledownthegunbore[23–25].Figure1illustrates
𝑝 𝑢 ⋅𝑢
schematicallytheconfigurationofatypicalNGGPScharge. =−𝑄 𝐴−𝑓𝑢 𝐴+𝑚̇ 𝐴(𝑒 + + 𝑝 𝑝)
Thistwo-phaseflowcodewasdevotedtosimulatetheevents 𝑝 𝑠 𝑝 𝑐 𝑝 𝜌 2
𝑝
during interior ballistic cycle inside a real large caliber
NGGPS.Thepredictednumericalresultswerevalidatedwith 𝜕(𝐴𝜑)
+∑𝑚̇ 𝐻 𝐴−𝑝 .
experimentaldata. ign ign 𝜕𝑡
(5)
2.TheoreticalModel
Inthepreviousequations,𝜑isthevolumefractionofthe
The physical events, criterion of ignition, flame spreading, gasphase,𝜌𝑔,𝜌𝑝aregasandsoliddensity,𝑢𝑔,𝑢𝑝aregasand
andcombustionofpropellantgrainsinsidethegunchamber solidvelocity,𝑃,𝑒𝑔arepressureandinternalenergyofthegas
can be referred to in [26]. Based on charge configuration phase,𝑚𝑐istherateofgasmassgenerationduetopropellant
showninFigure1,atwo-phaseflowmathematicalmodelwas combustion,𝑚 isthemassflowrateofgasfromthevent-
ign
establishedbasedonEulerian-Eulerianapproach.Interaction holes on the igniter, 𝐻 is stagnation enthalpy of the gas
ign
processesbetweenphasesareaddedtosingle-phaseconser- flowfromthevent-holesontheigniter,and𝑓𝑠,𝑄𝑝,and𝑅𝑝are
vation laws via source terms. The governing equations are interphase drag, interphase heat transfer, and intergranular
madeofthemassandmomentumconservationlawsforboth stress,respectively.
AdvancesinMechanicalEngineering 3
The interphase heat transfer equations (6)–(8) for the And the rate of gas mass generation due to propellant
granularpropellantsreplacethegoverningequationforsolid combustion was calculated depending on the propellant
phaseenergyandweredevelopedasfollows: burningrateequationasfollows:
𝜌 (1−𝜙)𝑆 𝑞 𝑑𝑟
𝑄 = 𝑃 𝑝 𝑝 =−𝑑̇
𝑃 𝑀 𝑑𝑡
𝑃
𝑆 =𝜎 𝑆 (6) 𝑑̇=𝑏𝑝𝑛
𝑝 𝑖 (10)
𝑀 =(1−𝜓)𝑀, 𝜌2(1−𝜙)𝑆 𝑑̇
𝑝 𝑖 𝑚̇ = 𝑝 𝑝 ,
𝑐 𝑀
𝑝
where𝜎,𝜓arethepropellant’srelativeburnedsurfaceratio
and the propellant’s relative burned quantity, respectively.
where 𝑏 is theburningrate coefficient and 𝑛 isthe burning
They were determined from the shape function relations
of seven perforated propellant grains. 𝑆𝑖, 𝑀𝑖 are the initial rateexponent.Theinteractionbetweenthetwophases,such
as interphase drag and intergranular stress, is necessary to
burning surface and the initial mass of propellant grain,
close the governing equations above. For more details, we
respectively.Consideringthemechanismsofconvectionand
radiation,theinterphaseheatexchange𝑞iscalculatedas referto[18].
𝑞=(ℎ +ℎ )(𝑇 −𝑇 ) 3.NumericalApproach
𝑝 re 𝑔 ps
0.4Re2𝑝/3Pr1/3𝜆𝑓 4𝛾 3.1. MacCormack Technique. Equations (1)–(5) for one-
ℎ𝑝 = 𝑑 , Pr= 9𝛾−5 dimensional two-phase reactive flow can be written in the
𝑝 form of conservation laws for simplifying the numerical
𝑑𝑝𝜌𝑔󵄨󵄨󵄨󵄨󵄨𝑢𝑝−𝑢𝑔󵄨󵄨󵄨󵄨󵄨 solution:
Re𝑝 = 𝜇 𝜕𝑈 𝜕𝐸
+ =𝑆, (11)
(7) 𝜕𝑡 𝜕𝑥
𝐶 𝑇3/2
1 𝑔
𝜇=
𝐶2+𝑇𝑔 where𝑈,𝐸,and𝑆aretheconservedvariables,thefluxvector,
andthesourcevector,respectively.
𝐶3𝑇𝑔3/2 Thecomponentsof𝑈aretheconservedvariablesandcan
𝜆 =
𝑓 𝐶 +𝑇 bewrittenas
4 𝑔
ℎ =𝐶 (𝑇 +𝑇 )(𝑇2+𝑇2), 𝜑𝜌 𝐴
re 5 𝑔 ps 𝑔 ps 𝑔
(1−𝜑)𝜌 𝐴
wttrhahenetrshefeeℎrrr𝑝meidasltbthyrearnhasedfaeiatrtticrooanne,ffisfReceri𝑝erenidts,bt𝜇hyeicsoRtnhevyeenmcotiloodlnsec,nuℎulrameribsvetihsr,ceo𝜆hs𝑓ietaiyst, 𝑈=(𝑈𝑈𝑈123)=(((((( 𝜑𝜌𝑔𝑢𝑔𝐴𝑝 )))))). (12)
𝑇𝑔isthegastemperature,and𝛾isthespecificheatratio.𝑇ps 𝑈4 (( (1−𝜑)𝜌𝑝𝑢𝑝𝐴 ))
isthesolidparticlesurfacetemperatureanditwascalculated 𝑈 ( )
5
as 𝑢2
𝑔
𝜑𝜌 (𝑒 + )𝐴
𝑔 𝑔 2
𝑑𝑇 2𝑞√𝑎 (√𝑡+𝛿𝑡−√𝑡) ( )
ps = 𝑝 , (8)
𝑑𝑡 𝜆 √𝜋 𝛿𝑡
𝑝 Andthecomponentsof𝐸arethefluxfunctions:
where𝑎𝑝 istheaccelerationofpropellantparticlesand𝜆𝑝 is 𝜑𝜌 𝑢 𝐴
thecoefficientofconductivity. 𝑔 𝑔
The Noble-Abel equation of state for the gas phase
(1−𝜑)𝜌 𝑢 𝐴
mixture was adopted in this formulationto account for the 𝐸 ( 𝑝 𝑝 )
1 ( )
nonidealbehaviorofthegasphase: 𝐸=(𝐸𝐸23)=(((( 𝜑𝜌𝑔𝑢𝑔2𝐴 )))). (13)
𝑝= 𝑅𝑇𝑔 , (9) 𝐸𝐸4 ((( (1−𝜑)(𝜌𝑝𝑢𝑝2 +𝑅𝑝)𝐴 )))
(1/𝜌 −𝛼) 5
𝑔 𝑢2 𝑝
𝑔
𝜑𝜌 𝑢 (𝑒 + + )𝐴
where𝑅isthegasconstantand𝛼isthegascovolume. ( 𝑔 𝑔 𝑔 2 𝜌𝑔 )
4 AdvancesinMechanicalEngineering
Andthecomponentsof𝑆arethesourcetermfunctions:
𝑚̇ 𝐴+∑𝑚̇ 𝐴
𝑐 ign
−𝑚̇ 𝐴
𝑆 ( 𝑐 )
1 ( )
𝑆 ( )
2 ( −𝑓𝐴+𝑚̇ 𝑢 𝐴+∑𝑚̇ 𝑢 𝐴−(𝜑𝐴)∇𝑝 )
𝑆=(𝑆3)=(( 𝑠 𝑐 𝑝 ign ign )). (14)
𝑆 ( )
4 ( 𝑓𝐴−𝑚̇ 𝑢 𝐴−(1−𝜑)𝐴∇𝑝 )
𝑆 ( 𝑠 𝑐 𝑝 )
5
𝑝 𝑢 ⋅𝑢 𝜕(𝐴𝜑)
−𝑄 𝐴−𝑓𝑢 𝐴+𝑚̇ 𝐴(𝑒 + + 𝑝 𝑝)+∑𝑚̇ 𝐻 𝐴−𝑝
𝑝 𝑠 𝑝 𝑐 𝑝 𝜌 2 ign ign 𝜕𝑡
( 𝑝 )
Thenonlinearhyperbolicsystemofdifferentialequationsof wheretheweightingfactor,𝑘,haslimitedvalueswhichcannot
theestablished two-phasereactiveflowmodel(11)issolved betoosmalltopreventresultsfluctuations.
using MacCormack technique which is an explicit finite Thenecessaryconditionforconvergenceofthepresented
differencemethodwithsecond-orderaccuracyinbothspace numerical scheme in this study is the Courant-Friedrichs-
and time via two steps, predictor and corrector [21]. At Lewy (CFL) stability criterion [32]. For our model this
eachstepofthetimemarchingsolution,asmallamountof conditioniscalculatedas
aaprtpifirocaiaclhviisscoapsiptylieisdatdodesmdaonodthatShheuflmuacntufialttieornesnihnanthceemperne-t 𝑐ΔΔ𝑥𝑡 ≤1, with 𝑐=max(󵄨󵄨󵄨󵄨󵄨𝑢𝑔󵄨󵄨󵄨󵄨󵄨+󵄨󵄨󵄨󵄨󵄨𝑢𝑝󵄨󵄨󵄨󵄨󵄨+|𝑎|), (20)
dictedresultsateachtimestep[27–31].
The predictor step based on the right-hand side with where𝑐isthelargestvelocityinthesystemand𝑎isthesound
forwarddifferenceis speed.
𝑈𝑡+Δ𝑡 =𝑈𝑡− Δ𝑡 (𝐸𝑡 −𝐸𝑡)+Δ𝑡𝑆𝑡+V𝑡, (15) 3.2. Adaptive Grid and Computational Domain. The finite
𝑖 𝑖 Δ𝑥 𝑖+1 𝑖 𝑖 𝑖
computationaldomainofthepresenttwo-phaseflowmodel
wheretheartificialviscositytermV𝑡iscalculatedas isdividedinto102gridpointsinthechamberregionbefore
𝑖
the projectile starts to move. Then, after the base pressure
V𝑖𝑡 = 𝐶𝑥𝑝󵄨󵄨󵄨󵄨𝑡𝑝𝑖𝑡+−1−2𝑝2𝑡𝑝+𝑖𝑡+𝑝𝑡𝑝𝑖𝑡−1󵄨󵄨󵄨󵄨(𝑈𝑖𝑡+1−2𝑈𝑖𝑡+𝑈𝑖𝑡−1), (16) rmeaocvhee.sThtheerefsotarret,tphreesesquureati(o3n0MofPma)ottihoenporfotjhecetiplreosjetacrtitlsetios
𝑖+1 𝑖 𝑖−1
essentialforrepresentingthecellextensionattheprojectile
where 𝐶𝑥 is an arbitrary specified parameter with typical basewhichcanbewrittenas
valuewithinrange(0.01–0.3). 𝑑V
𝑝𝐴=𝜑 𝑚 pr, (21)
And the corrector step based on the right side with 1 pr 𝑑𝑡
rearwarddifferenceandsubstitutingthepredictedvaluesof
thetimederivative𝑈𝑡+Δ𝑡attime𝑡+Δ𝑡is where𝑝ispressureattheprojectilebase,𝐴isthebarrelcross-
𝑖 sectionalarea,𝜑1isthecoefficientofsecondarywork(related
𝑈𝑡+Δ𝑡 = 1[𝑈𝑡+𝑈𝑡+Δ𝑡− Δ𝑡 (𝐸𝑡+Δ𝑡−𝐸𝑡+Δ𝑡)+Δ𝑡𝑆𝑡+Δ𝑡] to each gun system), 𝑚pr is the projectile mass, and Vpr is
𝑖 2 𝑖 𝑖 Δ𝑥 𝑖 𝑖−1 𝑖 the projectile velocity. At the moment when the projectile
startstomove,thedisplacementoftheprojectileintheaxial
+V𝑡+Δ𝑡, directioncanbewrittenas
𝑖
(17) 𝑥𝑗𝑛+1 =𝑥𝑗𝑛+VprΔ𝑡. (22)
where V𝑡+Δ𝑡 is the calculated artificial viscosity using the
𝑖 At that time, the last mesh cell of the computational
predictedvaluesasfollows:
domain begins to expand. To overcome the problem of
V𝑡+Δ𝑡 = 𝐶𝑥󵄨󵄨󵄨󵄨󵄨󵄨𝑃𝑡𝑖++1Δ𝑡−2𝑃𝑡𝑖+Δ𝑡+𝑃𝑡𝑖−+1Δ𝑡󵄨󵄨󵄨󵄨󵄨󵄨 disyanpapmliiecdehxperaenisni.oAnsotfhtehipsrcoejellc,tailegrmidovaedsa,ptthaetiolenngatlghoorifththme
𝑖 𝑃𝑡+Δ𝑡−2𝑃𝑡+Δ𝑡+𝑃𝑡+Δ𝑡 last mesh cell of the computational domain Δ𝑥2 will also
𝑖+1 𝑖 𝑖−1 (18) increaseasshowninFigure2andthegridwillbeself-adapted
𝑡+Δ𝑡 𝑡+Δ𝑡 𝑡+Δ𝑡 accordingtothesetwoconditions.
×(𝑈 −2𝑈 +𝑈 ).
𝑖+1 𝑖 𝑖−1 If Δ𝑥2 ≤ 1.5Δ𝑥, the number of the mesh cells is still
constant but the value of Δ𝑥1 is changed to be equal to the
TheShumanfilterapproachforthecalculatedconserved valueofΔ𝑥2.
variablescanbewrittenas If Δ𝑥2 ≤ 1.5Δ𝑥, a new mesh cell will be added and
𝑈𝑖𝑡 = 𝑘+11(𝑈𝑖𝑡−1+𝑘𝑈𝑖𝑡+𝑈𝑖𝑡+1), (19) tohfeΔv𝑥a2lu−eΔo𝑥f.Δ𝑥1 will be calculated to be equal to the value
AdvancesinMechanicalEngineering 5
j−2 j−1 Δx j At the start pressure value, the projectile starts to move
2
tn+1 whiletheleftboundaryisstillstationary.Thetwo-phaseflow
confinedinthechamberbehindtheprojectilestartstomove
asthepressureandchambervolumeareincreasing.MCVC
Trajectory of
method is established to handle the computational domain
projectile motion
attherightboundary.Applyingthemassconservationonthe
tn controlvolumeofthelastmeshcellbehindtheprojectilebase
Δx asshowninFigure2,theporosityofpoint𝑗attime𝑡𝑛+1can
j−2 j−1 1 j
becalculatedviathefollowingprocedure.
Figure2:Schematicofmovingcontrolvolume. Thesolidphasemassesinsidethecontrolvolumeattimes
𝑡𝑛and𝑡𝑛+1arecalculatedas
Table1:Inputdatausedinnumericalcomputation. 𝑚𝑝𝑛 =[𝐴𝜌𝑝(1−𝜑)]𝑛𝑗−1/2(𝑥𝑗𝑛−𝑥𝑗𝑛−1), (24)
Parameter Value
Geometry 𝑚𝑝𝑛+1 =[𝐴𝜌𝑝(1−𝜑)]𝑛𝑗−+11/2(𝑥𝑗𝑛+1−𝑥𝑗𝑛−+11). (25)
Guncaliber(m) 0.130
Barrellength(m) 6.29 Themassofsolidphasethatflowsintothesectionarea𝐴is
Chamberlength(m) 0.7 givenby
Chambervolume(m3) 0.01366
𝑛+1/2
Projectilemass(Kg) 33.4 Δ𝑚𝑝 =[𝐴𝜌𝑝(1−𝜑)𝑢𝑝]𝑗−1/2Δ𝑡. (26)
Igniterlength(m) 0.5
Igniterdiameter(m) 0.16
The burnet mass of solid phase from the control volume
Numberofholesoftheigniter 60 duringtimeΔ𝑡is
Holediameter(m) 0.004
Physicalproperties Δ𝑚̇ =−𝐴𝑚̇𝑛(𝑥𝑛−𝑥𝑛 )Δ𝑡. (27)
br 𝑐 𝑗 𝑗−1
Blackpowdermass(Kg) 0.08
Propellantchargemass(Kg) 10.6
Then, the mass of the solid phase in the control volume at
Propellantcovolume(m3/kg) 0.001 time𝑡𝑛+Δ𝑡canbegivenas
Propellantdensity(kg/m3) 1600
PPrrooppeellllaannttfigonrciteio(Jn/ktegm) perature(K) 906510500 𝑚𝑝𝑛+1 =𝑚𝑝𝑛 +Δ𝑚𝑝+Δ𝑚̇br. (28)
Propellantgeometry Cylindrical
In the Eulerian coordinate system, the next relation can be
Initialporosity 0.56
deduced:
Startpressure(MPa) 30
Burningratecoefficient 1.257∗10−8 𝑥𝑗𝑛−+11 =𝑥𝑗𝑛−1. (29)
Burningrateexponent 0.83
Finally,from(25)and(29),theporositycanbeexpressedas
3.3.StationaryandMovingBoundaryConditions. Inthisone-
𝑚𝑛+1
dimensional two-phase flow model, a criterion of reflective 𝜑𝑛+1 =1− 𝑝 .
boundary condition is utilized at the left (breech position) 𝑗−1/2 𝐴𝜌 (𝑥𝑛+1−𝑥𝑛 ) (30)
𝑝 𝑗 𝑗−1
andright(projectilebase)boundariesbeforetheshotstarts
tomove.Figure3representstwofictitiouspointsatthegrid
points 𝑖 = 0 and 𝑖 = 𝑁 + 1, respectively. These points are But
requiredinthecalculationsatthefixedboundariesandthey 1
are located outside the left and right physical domains. The 𝜑𝑗𝑛−+11/2 = 2(𝜑𝑗𝑛−+11+𝜑𝑗𝑛+1). (31)
primitivevariablesofbothphasesarecalculatedatthefixed
boundaryas Thus,anexpressionfortheporosityofpoint𝑗atthetime𝑡𝑛+Δ𝑡
canbepresentedas
𝑢 =−𝑢 , 𝑢=(𝑢 ,𝑢 ),
0 1 𝑔 𝑝
𝑞 =𝑞 , 𝑞=(𝜌 ,𝑝,𝜑,𝑇), (23) 𝜑𝑗𝑛+1 =2𝜑𝑗𝑛−+11/2−𝜑𝑗𝑛−+11. (32)
0 1 𝑔
where𝑢isthevelocityvectorforbothphasesand𝑞isavector Applyingthesamepreviousstepsusingenergyandmomen-
representing the other primitive variables of the flow (gas tum conservation equations, the remaining primitive flow
density,pressure,porosity,andtemperature). variables(𝑢𝑔𝑛+1,𝑢𝑝𝑛+1,𝜌𝑔𝑛+1,and𝑒𝑔𝑛+1)canbederived.
𝑗 𝑗 𝑗 𝑗
6 AdvancesinMechanicalEngineering
Projectile base
Breech position
i=0 1 2 N N+1
Figure3:Schematicofcomputationaldomain.
4.ResultsandDiscussion downthebore.Astheprojectilemovesararefactionwaveis
formedandtravelsintheoppositedirectiontowardsthegun
4.1.NumericalVerifications. Acommontestfortheaccuracy breech(firstrarefaction)duetothechambervolumeincrease
of computational fluid dynamics codes, one-dimensional following the projectile motion. But the combustion rate at
model introduced by G. A. Sod (shock tube) [22], is used the breech position is higher than the traveling velocity of
to test the ability of MacCormack algorithm to solve the the rarefaction wave at that time; therefore a second shock
IBVPforthesystemoffirst-order,nonlinearcoupledpartial waveisformedandtravelsagaintotheshotbasefollowedby
differential equations with inhomogeneous term possessing the second rarefaction wave. Similarly, the third shock and
shockwavebehavior.Thistestconsistsofaone-dimensional rarefactionwavesareformedandtravelbetweenthebreech
Riemannproblemforatubeofunitlength−1 ≤ 𝑥 ≤ 1[33]. and the shot base. As the pressure increases, the projectile
TheCFLnumberistakentobe0.7forstabilityrequirements velocity increases resulting in attenuation to the pressure
of the solution and the final time is 𝑡 = 0.2. The initial wave strength. The chamber pressure reaches its maximum
conditionsfortheleftandrightstatesat𝑡=0are atabout8.4msandthenstartstodecreasegraduallyuntilthe
endoftheinteriorballisticscycle.
1.00, for 𝑥≤0.0,
𝜌(𝑥,0)={
0.125, for 𝑥>0.0,
4.2.2.Two-PhaseTemperatureProfile. Figure6representsthe
1.00, for 𝑥≤0.0, temperature distribution of the gas phase and propellant
𝑝(𝑥,0)={0.10, for 𝑥>0.0, (33) particles (solid phase). It is shown from Figure 6(a) that
thegastemperatureisgraduallyincreasedintheearlystage
0, for 𝑥≤0.0, ofcombustioncomparedtothepropellantparticlesthathave
𝑢(𝑥,0)={ amorerapidburningrate.Attimeabout8.6msthegastem-
0, for 𝑥>0.0.
peraturereachesitsmaximumvalue.Alsoitisclearthatthe
temperature distribution near the breech position is higher
Figure4showsthecomputedresultsofthetestcompared
thanthatbehindtheshotbase.Withtheprojectilemotionthe
totheexactsolution.Itisclearfromthisfigurethatthereis
gastemperatureslowlydecreases.Theignitionwavedirection
agoodagreementbetweenthenumericallycalculatedresults
ofpropagationthroughthepropellantchargebedaswellas
andtheexactsolutionoftheEulersystemofequations.
thewavefrontisclearlyrepresentedinFigure6(b).Theentire
propellantchargebedistotallyignitedattimeabout3.25ms
4.2.ApplicationoftheModeltoaRealLargeCaliberNGGPS.
asitpossessespropellantsurfacetemperatureofabout615K.
The presented gas-solid reactive flow model in this work
is applied to real 130mm naval gun launching a guided
projectileutilizingpropellantchargewithignitionsystemas 4.2.3. Two-Phase Velocity Profile. The predicted velocity
described schematically in Figure 1. The calculated results distributions of gas and solid phase during the interior
of some interior ballistic parameters are compared to the ballistic cycle are shown in Figure 7. It is shown that there
experimentaldata. Theinputdata forcomputationscarried arenegativevelocitiesinvent-holesregionintheearlystage
outbythecodeisgiveninTable1. ofignitionattimesabout1.7msto4ms.Thecentralreason
of these negative velocities is the tendency of gas phase
combustionproductstopropagatetowardsthebreechandthe
4.2.1.PressureHistory. Figure5showsthecalculatedpressure
projectilebaseduetopressuregradientinsidegunchamber.
history during the interior ballistics process. As the jet
Butregardingtheirhighinertia,thesolidphaseparticlesare
of burned gases flows out of the igniter holes into the
clearlyobservedtohavelowervelocitythanthatofgasphase.
mainpropellantcharge,thepropellantaroundthevent-holes
Thislagofvelocityisobviousinthenegativevelocityregion
region starts to burn first and the pressure increases in the
aswellasduringthewholeinteriorballisticscycle.
gun chamber at about 2.4ms as shown in Figure 5(a). Due
to the pressure gradient, the flame continues to propagate
viatheunburnedregions.Atthattimeastrongshockwave 4.3. Experimental Validation. Due to complication and
(first shock) occurs towards the direction of projectile base transient behavior of the interior ballistics process, many
asshowninFigures5(c)and5(d).Byincreasingtheignited parameters cannot be measured experimentally except the
area, the flame continues to propagate and pressure in the chamber pressure and the projectile velocity at the muzzle.
chambercontinuouslyincreases.Atabout2.9msthepressure Table 2 represents the comparison between experimental
reaches the start value and the projectile starts to move and computational results for maximum chamber pressure
AdvancesinMechanicalEngineering 7
1.0 1.0
0.8
0.8
0.6
Pressure 0.6 Velocity 0.4
0.4
0.2
0.2
0.0
0.0
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
X(m) X(m)
Exact Exact
Test1 Test1
(a) (b)
1.0
0.8
y 0.6
nsit
e
D
0.4
0.2
−1.0 −0.5 0.0 0.5 1.0
X(m)
Exact
Test1
(c)
Figure4:Comparisonbetweenthenumericalandexactsolutionin1DRiemannproblem.
Table2:Comparisonbetweennumericalresultsandexperimentalresults.
Ballisticparameter Experimentalresults Simulationresults Errorprobability
Maximumchamberpressure(MPa) 319.00 319.45 0.14%
Muzzlevelocity(m/s) 857 861.5 0.52%
andprojectilemuzzlevelocity.Thevalidationshowsagood guidedprojectiles.Therefore,atwo-phaseflowmathematical
agreementbetweenexperimentalandcomputationalresults. model for the solid granular propellant and its products
of combustion inside large caliber NGGPS during interior
5.Conclusions ballisticcyclewasdeveloped.Theignitionofballpropellant
intheignitertubeiscoupledwiththetwo-phaseflowmodel
Themainobjectiveofthisresearchistopredicttheinterior ofthegunchamber.Theusednumericalmethodisconfirmed
ballistics parameters such as chamber pressure and muzzle andvalidatedsuccessfullyusingSodshocktubetestinorder
velocitytosustainanaccurateandsafelaunchingprocessfor to check its accuracy and aptitude to solve the complex
8 AdvancesinMechanicalEngineering
350 18 Breech Muzzle
75 32 320
300 16 288
107
250 14 128 128 256
160
Pressure (MPa) 121500000 Time (ms)118620222432152986917 252603171prTorjaejcetcilteo rmyo otfion 211129622048 Pressure (MPa)
96
50 4 64 43 64
2 11 32
0
0
0 2 4 6 8 10 12 14 16 18 0 1 2 3 4 5 6
Time (ms) Distance (m)
Breech pressure
Projectile base pressure
(a) Pressurehistoryversustime (b) Pressureprofileon𝑥-𝑡diagram
320
30 First Second Third 8 Third Pressure Third 288
20 rarefactionrarefactionrarefaction shock decrease rarefaction 256
wave (MPa) 100 me (ms) 64 Ssehcoocnkd First raSreefcaocntidon 211296420 ure (MPa)
Pressure −−1200 sFhiorsctk Ssehcoocnkd sThhoircdk Ti 2 shock rareFfiarcsttion 1962648 Press
Pressure
increase 32
−30
0 0
0 2 4 6 8 10 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Time (ms) Distance (m)
(c) Pressurewavetimediagram (d) Initialpressurewavedistribution
Figure5:Pressurehistoryduringinteriorballisticcycleinthe130mmNGGPS.
Breech Muzzle
18 2450 4.0 615
1586 1514
16 1586 1370 2306 3.5 575
14 1658 1298
Time (ms) 112086 222130849002 65011877340 prTorjeacjetcilteo mryo otfion 2111074115348024 mperature (K) Time (ms) 3221....5050 456Wfroanvet 36Ip1grnoitpiaogna wti5oa9vn1e 544439516566 mperature (K)
4 1874 Igniter vent-holes 866 Te 1.0 376 Te
2 region 578 0.5 337
434
290 297
0 1 2 3 4 5 6 0.0 0.1 0.2 0.3 0.4 0.5 0.6
Distance (m) Distance (m)
(a) Gasphasetemperature (b) Solidphasetemperature
Figure6:Numericalresultsofthetwo-phasetemperaturedistributionon𝑥-𝑡diagram.
AdvancesinMechanicalEngineering 9
1186 Bree1c5h6 412 572 636764 796M82u8zzle 876302 1186 Bre5e4ch162 227 400421 576280 72M4uzzle 871002
60 316 700
ms) 111420 124253248444 540604Trajectory of 643074468 city (m/s) ms) 111420 124005270 3476844645T7r2ajectory of 549864 ocity (m/s)
Time ( 86 92 projectile motion 29220 hase velo Time ( 86 1176 projectile motion 327780 phase vel
4 −36 Igniterre vgeionnt-holes −−31664 Gas p 4 Igniter vent-holes 162 Solid
region
2 −292 2−32 54
−420 −54
0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
Distance (m) Distance (m)
(a) Gasphasevelocity (b) Solidphasevelocity
Figure7:Numericalresultsofthetwo-phasevelocitieson𝑥-𝑡diagram.
two-phase flow problem. A grid adaptation algorithm was [5] J. Z. Ming, Interior Ballistics in Guns, Beijing Institute of
developedandusedtodealwiththeprojectilemotionwith Technology,Beijing,China,2004.
MCVC boundary approach. The main advantages of this [6] A.W.Horst,E.K.George,andP.S.Gough,“Newdirectionsin
model are (1) predicting times and distributions of shock multiphaseflowinteriorballisticmodeling,”BRLReport3102,
waveoccurrencethatcandisturbtheguidedprojectileand(2) USAARDC,BallisticResearchLaboratories,AberdeenProving
estimatingaccurateandfastinteriorballisticparameters.The Ground,Md,USA,1990.
performanceofa130mmNGGPSwascloselypredictedand [7] A.J.BudkaandJ.D.Knapton,“Pressurewavegenerationingun
agoodagreementbetweenexperimentalandcomputational systems:asurvey,”Tech.Rep.BRL-MR-2567,BallisticResearch
resultswasobtained. Laboratories,AberdeenProvingGround,Md,USA,1975.
[8] A. W. Horst, I. W. May, and E. V. Clarke, “The missing link
ConflictofInterests betweenpressurewavesandbreechblows,”Tech.Rep.ARBRL-
MR-02849,USAARRADCOM,BallisticResearchLaboratory,
The authors declare that there is no conflict of interests AberdeenProvingGround,Md,USA,1978.
regardingthepublicationofthispaper. [9] I.W.MayandA.W.Horst,“Chargedesignconsiderationsand
their effect on pressure waves in guns,” Tech. Rep. ARBRL-
TR-02277, USA ARRADCOM, Ballistic Research Laboratory,
Acknowledgments
AberdeenProvingGround,Md,USA,1980.
The research was supported by the Research Fund for [10] P. S. Gough and F. J. Zwarts, “Modeling heterogeneous two-
phasereactingflow,”AIAAJournal,vol.17,no.1,pp.17–25,1979.
the Natural Science Foundation of Jiangsu province
(BK20131348) and Key Laboratory Fund (Grant no. [11] M.R.Baerand J.W.Nunziato, “Atwo-phasemixturetheory
9140C300103140C30001),China. forthedeflagration-to-detonationtransition(DDT)inreactive
granular materials,” International Journal of Multiphase Flow,
vol.12,no.6,pp.861–889,1986.
References
[12] C. R. Woodley, S. Billett, C. Lowe, W. Speares, and E. Toro,
“TheFHIBSinternalballisticscode,”inProceedingsofthe22nd
[1] M.M.Rashad,X.B.Zhang,andH.Elsadek,“Numericalsim-
InternationalSymposiumonBallistics,pp.322–329,Vancouver,
ulationofinteriorballisticsforlargecaliberguidedprojectile
Canada,November2005.
navalgun,”JournalofEngineeringandAppliedScience,vol.60,
no.2,2013. [13] R. Acharya and K. Kuo, “Implementation of approximate
[2] P.G.BaerandJ.M.Frankle,“Thesimulationofinteriorballistic riemannsolvertotwo-phaseflowsinmortarsystems,”Journal
performanceofgunsbydigitalcomputerprogram,”BRLReport ofAppliedMechanics,vol.77,no.5,ArticleID051401,9pages,
1183,ARDC,BallisticResearchLaboratories,AberdeenProving 2010.
Ground,Md,USA,1962. [14] C.ChengandX.Zhang,“Modelingofinteriorballisticgas-solid
[3] K.D.FickieandJ.Grosh,“Atechniqueforcodeaugmentation,” flow using a coupled computational fluid dynamics-discrete
Tech.Rep.BRL-MR-3622,USArmyBallisticResearchLabora- elementmethod,”JournalofAppliedMechanics,vol.80,no.3,
tory,AberdeenProvingGround,Md,USA,1987. ArticleID031403,6pages,2013.
[4] W.R.FredrickandS.R.Timothy,“Alumped-parameterinterior [15] P. S. Gough, “Initial development of core module of next
ballistic computer code using the TTCP model,” BRL Report generationinteriorballisticmodelNGEN,”Tech.Rep.ARL-CR-
3710,ARDC,BallisticResearchLaboratories,AberdeenProving 234,USArmyResearchLaboratory,AberdeenProvingGround,
Ground,Md,USA,1988. Md,USA,1995.
10 AdvancesinMechanicalEngineering
[16] P.S.Gough,“Formulationofanext-generationinteriorballistic [33] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid
code,”inProceedingsofthe28thJointArmyNavyNasaAirForce Dynamics: a Practical Introduction, Springer, New York, NY,
Combustion Subcommittee Meeting,pp.321–337,SanAntonio, USA,3rdedition,2009.
Tex,USA,November1991.
[17] M.J.NuscaandP.S.Gough,“Numericalmodelofmultiphase
flowsappliedtosolidpropellantcombustioningunsystems,”in
Proceedingsofthe34thAIAA/ASME/SAE/ASEEJointPropulsion
ConferenceandExhibit,Cleveland,Ohio,USA,July1998,AIAA
paperno.98-3695.
[18] Y.X.YuanandX.B.Zhang,MultiphaseHydrokineticFoundation
of High Temperature and High Pressure, Harbin Institute of
Technology,Harbin,China,2005.
[19] C.MaandX.Zhang,“Simulationofcontaminationprevention
for optical window in laser ignition systems of large-caliber
guns,”JournalofAppliedMechanics,vol.78,no.5,ArticleID
051014,2011.
[20] J.S.Jang,H.G.Sung,T.S.Roh,andD.W.Choi,“Numerical
study on properties of interior ballistics according to solid
propellant positions in chamber,” in Proceedings of the 26th
InternationalSymposiumonBallistics,pp.721–730,Miami,Fla,
USA,September2011.
[21] J.D.Anderson,ComputationalFluidDynamics:theBasicsand
Applications,McGraw-Hill,NewYork,NY,USA,1995.
[22] G. A. Sod, “A survey of several finite difference methods for
systemsofnonlinearhyperbolicconservationlaws,”Journalof
ComputationalPhysics,vol.27,no.1,pp.1–31,1978.
[23] S.G.Ahmed,“Anewalgorithmformovingboundaryproblems
subjecttoperiodicboundaryconditions,”InternationalJournal
ofNumericalMethodsforHeat&FluidFlow,vol.16,no.1,pp.
18–27,2006.
[24] I.Demirdzˇic´andM.Peric´,“Finitevolumemethodforpredic-
tion of fluid flow in arbitrarily shaped domains with moving
boundaries,” International Journal for Numerical Methods in
Fluids,vol.10,no.7,pp.771–790,1990.
[25] A. van Dam and P. A. Zegeling, “A robust moving mesh
finite volume method applied to 1D hyperbolic conservation
lawsfrommagnetohydrodynamics,”JournalofComputational
Physics,vol.216,no.2,pp.526–546,2006.
[26] H.Krier,S.Rajan,andW.F.vanTassell,“Flamespreadingand
combustioninpackedbedsofpropellantgrains,”AIAAJournal,
vol.14,no.3,pp.301–309,1976.
[27] K. Alexander and L. Yu, “New adaptive artificial viscosity
methodforhyperbolicsystemsofconservationlaws,”Journal
ofComputationalPhysics,vol.231,no.24,pp.8114–8132,2012.
[28] S.KawaiandS.K.Lele,“Localizedartificialviscosityanddiffu-
sivityschemeforcapturingdiscontinuitiesoncurvilinearand
anisotropic meshes,” Center for Turbulence Research Annual
ResearchBriefs,2007.
[29] E.J.Caramana,M.J.Shashkov,andP.P.Whalen,“Formulations
ofartificialviscosityformulti-dimensionalshockwavecompu-
tations,”JournalofComputationalPhysics,vol.144,no.1,pp.70–
97,1998.
[30] S. P. Nelson and M. L. Weible, “Three-dimensional Shuman
filter,”JournalofAppliedMeteorology,vol.19,no.4,pp.464–469,
1980.
[31] A.HartenandG.Zwas,“SwitchednumericalShumanfiltersfor
shockcalculations,”JournalofEngineeringMathematics,vol.6,
no.2,pp.207–216,1972.
[32] R. Courant, K. O. Friedrichs, and H. Lewy, “On the partial
differenceequationsofmathematicalphysics,”IBMJournalof
ResearchandDevelopment,vol.11,pp.215–234,1964.