Solid Composite Propellants (SCPs) are extensively used in the ﬁeld of propulsion for their chemical and mechanical stability in long-term storage, and for having simple production and operation processes. Computational modeling enables cost reduction, increased eﬃciency, and greater coverage of the conﬁguration space in the SCP design process. However, accurate and eﬃcient SCP modeling presents a number of numerical challenges. A primary obstacle in modeling these systems is capturing the complex evolving interface. Recently, it has been shown that the phase-ﬁeld method has a strong ability to model the combustion behavior of SCPs, implicitly capturing the topological evolution at a relatively low computational cost. Initial phase-ﬁeld methods show promise in their ability to do predictive regression modeling but require a number of approximations and heuristic modeling methods. This work presents a new formulation for the phase ﬁeld regression model that combines a thermal solver with an Arrhenius rate law to model interface regression using a comprehensive and physics-based approach. This improves the capability of previous methods by increasing the number of kinematic forces that are accounted for, and allowing the study of thermal diﬀusivity in the system.Italsoenablestheintegrationofafullycoupledsolid-ﬂuidinterface.Todemonstratethe eﬃcacyofthemodel,itisappliedandcalibratedtoahomogeneousmonopropellant(ammonium perchlorate). The model is validated for a range of temperatures, with a reasonable quantitative match to experimental data. It was also demonstrated that the model recovered the relationship between ignition time and heat ﬂux, with no ﬁtting required. Overall, the method showed great capability of reproducing experimental data by matching temperature, burn rate, and thermal diﬀusivity proﬁles. NOVELTY AND SIGNIFICANCE A new diﬀuse interface method is developed for calculating ignition and regression rates in SCPs. While this builds on previous work in diﬀuse interface modeling of SCPs, it is novel in its coupling to thermal transport, as well as its development of (and coupling to) a surrogate model for gas phase heat ﬂux. The model is shown to reproduce experimentalregressionratedatawithreasonableaccuracy.Moreover,itquantitativelycapturesignitionbehaviorwith noadditionalmodiﬁcationstothemodel.Themodelisimplementedinahighperformnacecomputationalframework that is able to resolve large, three-dimensional mesostuctures. AUTHOR CONTRIBUTIONS MM:Performing research, writing. ES:Performing research. PM:Performing research. JMQ:Supervision and writing.BR:Supervision and writing. 1. Introduction Solid Composite Propellants (SCPs) have been used as rocket and missile propellants for over a century. The earliest SCPs were composed of a mixture of charcoal, potassium nitrate, and sulfur [1, 2]. Over time, they have been further developed to enhance their propulsion properties. Despite increased eﬃciency and controllability of liquid propellants in the early 20th century, SCPs are still used because of their simplicity, reliability, longevity, and <Corresponding author brunnels@iastate.edu (B. Runnels) ORCID(s):0000-0001-9914-8811 (M. Meier); 0000-0003-0856-6343 (E. Schmidt); 0000-0003-4650-2930 (J. Matt Quinlan); 0000-0003-3043-5227 (B. Runnels) Meier, Schmidt, Martinez, Quinlan, Runnels Page 1 of 16 Diﬀuse interface method for solid composite propellant ignition and regression APHTPB Main diffusion flame AP decomposition Lean premixed flame Rich premixed flame 830 K1400 K3000 K Figure 1: Conceptual ﬂame structure of a burning AP/HTPB SCP. The AP (orange) particles undergo phase change and closely coupled self-deﬂagration reactions (light pink), which produce oxygen as a product. Heat transfer to the HTPB (navy) causes polymer degradation and decomposition into gaseous, short-chain hydrocarbons. These two regions of gas diﬀuse together creating a triple-point ﬂame (lean and rich premixed ﬂames shown in light red and light blue and primary diﬀusion ﬂame shown in red). The diﬀusion ﬂame typically encapsulates the AP regions due to the overall rich mixtures. Owing to the high temperature of the diﬀusion ﬂame and the low liftoﬀ height near the triple-ﬂame, the heat transfer to the solid is locally high, causing high local regression rates. cost [3]. Modern SCPs mainly consist of a mixture of an oxidizer and a fuel, with additional chemical additives to improve stability, cross-linking, and thrust [3]. A common SCP is the combination of ammonium perchlorate (AP) foranoxidizerandhydroxyl-terminatedpolybutadiene(HTPB),whichactsbothasabinderandafuel.ModernSCPs often have additive components such as aluminum particles to increase performance. Propellants such as AP/HTPB have a heterogeneous mesostructure, leading to a complex and time-varying burn behavior at the mesoscale (Fig. 1). Speciﬁcally,attheinterfacesbetweentheAPandHTPBalifteddiﬀusionﬂamesitsclosetothesolidsurfaceresulting in high heat transfer back to the solid. This high heat transfer locally promotes faster surface regression rates. The couplingbetweenSCPcompositionandburnratebehaviorisdiﬃculttoquantify,andhashistoricallybeendoneinan Edisonian approach, through costly and hazardous full-scale tests [4]. Computational simulations, on the other hand, are essential to eﬃciently and eﬀectively design propellants with reduced environmental impact and lower cost. Since the ﬂame structure and chemical mechanisms are complex and the nonlinear interface topology must be accountedfor,itisachallengetocomputationallysimulatetheregressionofSCPs.Diﬀerentmethodshavepreviously beenproposed[5,4,6,7,8,9],buttheygenerallyassumeasimpliﬁedinterfacialgeometryinordertomaketheproblem tractable. Again, one of the main issues is that SCPs have a heterogeneous structure at the mesoscale, formed mostly by agglomerations of AP particles in a HTPB matrix. Since these materials decompose diﬀerently and the resulting ﬂame structure depends on the concentration of the fuel and oxidizer, which in turn depends on the location of the materials’ interface. This interface, where vigorous reaction and high heat transfer occurs, is geometrically complex andever-changing.Furthermore,thepreferentialregressionneartheinterface,duetotheheattransfer,tendstoforman undulating surface geometry. This can lead to another diﬃculty in regression simulations, the solid surface topology can even change. While AP is a solid crystalline compound that is stable at room temperature, at elevated temperatures it can act as an oxidizer or even react as a monopropellant. The speciﬁc decomposition mechanism of AP is complex due to its pressure and temperature dependence. At low pressure, the crystal structure of AP undergoes phase transformation before decomposing, but at pressures higher than 2MPa the material directly sublimates. In this case, the AP can supportself-deﬂagration.Theﬂametemperatureandsurfacetemperaturehavebeenfoundtobeapproximately1400K and800Krespectively[6,10].HTPB,ontheotherhand,isapolymercompoundformedbyverylongcarbonchains. It is commonly used in SCPs, in part due to its mechanical properties and that it has a curing point lower than the decomposition temperature of AP [11]. The HTPB is gasiﬁed in a pyrolysis process where the long carbon chains breakintosmallerhydrocarbons.SincethisdecompositionprocessissigniﬁcantlyslowerthantheburningofpureAP orthecombustionofAP/HTPBattheirinterface,chunksofbindercanbeejectedfromtheSCPduringthecombustion Meier, Schmidt, Martinez, Quinlan, Runnels Page 2 of 16 Diﬀuse interface method for solid composite propellant ignition and regression [12,13].Additionally,APandHTPBaredenseandhaveverylowthermalconductivity,actingasinsulatorsduringthe burningofSCPs.Asaresult,ithasbeenobservedthatheatpropagationduringcombustiononlyaﬀectsthetemperature a short distance ahead of the surface while a large portion of the materials’ temperature is unchanged [14]. AlthoughtheburnrateofSCPmixturesofAP/HTPBhavebeenexperimentallymeasuredinthepastusingoptical and break-wire methods [15, 16], these methods can be costly and hazardous, and do not provide surface topology information. Therefore to understand the interfacial burning of a SCP more fundamentally, a sandwich structure is often used in which layers of fuel and oxidizer are burned. Speciﬁcally, the test samples have a layer of HTPB sandwichedbetweenlayersofAPoneitherside.Comparedwithcomplexmixturesinaparticle/matrixstructure,these experimentsareeasiertoperformandtheresultsaremorerepeatable.IncontrasttotheﬂamestructureforanAP/HTPB SCPthatisthree-dimensionalandinherentlyunsteady,thesandwichpropellantexperiments’ﬂamestructureisnearly two-dimensional —which simpliﬁes observations of the ﬂame—and is nearly steady. As a result, the detailed surface geometry in relation to the interface can be quantiﬁed. In previous studies, the description of this combustion process has been divided into low-pressure premixed ﬂame, where gasiﬁed binder and oxidizer are mixed at the molecular level,andhigh-pressurecombustionwheretheﬂamebehaviorssplitintooxidizingproducts+ammoniaandintothree diﬀerent ﬂames: (i) premixed binder rich ﬂame, occurring at low pressure (ii) Oxidizing Products + Ammonia, and (iii) Binder products + AP oxidizing products. Where (ii) and (iii) occur at high pressure [17, 18]. One reason that simulations have been used to predict the surface regression of burning SCPs is to reduce the environmental impact of the combustion and improve their design at a relatively low cost. Modeling this regression requires two major steps: (i) evaluate the regression rates of the propellants and (ii) evolve the surface morphology. Traditionally, the regression rates have been evaluated using a pressure power-law in the form of r= .P_P0/nas proposed by Summerﬁeld [19]. Although this method has performed well in previous studies, it is obtained by ﬁtting the experimental data and fails to capture the physical unsteady ﬂame behavior and its inﬂuence on the regression. Furthermore, predictions cannot be made on SCPs where the power-law behavior has not been experimentally measured. Another issue is that tracking the surface morphology can be particularly challenging for geometries such as packed AP particles in a HTPB matrix, where the interface between the two materials is complex. The traditional approach for this problem has been to use an explicit approach, parameterizing and tracking the interface as it moves through the simulation’s mesh grid. In such cases, the solid/gas boundary, or unburnt/burnt boundary, is treated as a sharp interface [20, 4]. The level-set method is the most well-developed approach that uses the Eulerian view of thesurfaceregression[5,21,22].Thismethodapproximatesanisocountourmappingbasedonlocalregressionrates, thus it eﬀectively also uses a sharp interface model. The most severe limitations of the level-set methods are that it is very diﬃcult to capture topological interface changes and the underlying assumption in the level-set model is that the interface geometry is near-planar. In order to avoid these issues, the phase-ﬁeld method for simulating regression was recently introduced by Kanagarajan et al.[12], in which a diﬀuse interface model is used to track the solid/gas interface. In the phase-ﬁeld method, an order parameter ﬁeld is used (usually called .x/) to represent the domain occupied by the unburned materialasthesupportof =1.Iftherangeof is[0;1],thenthesupportof =0representsthedomainoccupiedby theburned/gasphasematerials.Thboundaryisconsequentlyrepresentedimplicitlyastheregioninwhich 0<<1, with a (usually controllable) diﬀuse width of . The advantage of this approach, similar to the level-set method, is that the implicit boundary eﬃciently captures complex morphology and topological transitions. Unlike the level-set method’sﬁeldfunction,however,thephase-ﬁeldmethod’sorderparametercanbeinterpretedphysicallyasapointwise reactioncoordinatethatevolvesaccordingtoathermodynamicpotential.Suchapotential,generallycalledafreeenergy functional, may include contributions from multiple physical systems, such as chemical energy, stored mechanical energy, and thermal energy. The diﬀuse boundary formulation is also readily amenable to solving PDEs driven by complexboundaries,includingsolidmechanics[23],ﬂuidmechanics[24],andheattransport(discussedbelow).Since the diﬀuse interface must be adequately resolved, the diﬀuse boundary methods do incur additional computational cost.Thiscostmaybeeﬀectivelylimited,though,throughtheuseofadaptivemeshreﬁnement.Limitationsfordiﬀuse interface methods are normally associated with length scale parameters eﬀects over the solution and computational complexity and cost. As it has been demonstrated in reference [12], the phase-ﬁeld equations can be formulated to assure convergence of the simulation and that advanced meshing techniques can be used with the model to maintain computational eﬃciency. Thisworkpresentsaphase-ﬁeldmodelfortheregressionoftheburninterfaceduringignitionanddeﬂagrationof an SCP with arbitrary morphology. This work builds on the previously introduced phase-ﬁeld model in [12] is used, but the heuristic pressure power-law kinetics are replaced with a “full-feedback” coupling that accounts for the eﬀect Meier, Schmidt, Martinez, Quinlan, Runnels Page 3 of 16 Diﬀuse interface method for solid composite propellant ignition and regression gas solid ap htpb Figure 2: Problem schematic illustrating the role of the two ﬁeld variables and .=1represents the solid/unburned region and is the initial condition setting for the problem, while = 0implicitly deﬁnes the gaseous region. The diﬀuse interfaceAlocated in the region 0< <1. c). =1indicates the oxidizer AP, while =0indicates the binder HTPB. The transition region between oxidizer and binder, indicated by 0< <1, is also diﬀuse in order to smoothly transition material properties with no grid dependence. [12] of mass and energy exchange across the boundaries between the ﬂuid and solids phases. This reduces the number of ﬁtting parameters, increases the physicality and predictiveness of the model. It also enables the ignition process to be simulated, which was previously not possible, and provides the complete thermal response thoughout the domain. Thoughnotexploredhere,italsomakesitpossibletocoupletoareactingﬂuidsolversothatthecombustioneﬀectson theregressioncanbedirectlycalculated.Whiletheimplementationofaﬂuidsolverisbeyondthescopeofthiswork, a ﬂuid stand-in has been implemented as a placeholder to both simulate the eﬀects of that solver and also to assure convergence of the solution. The remainder of this paper is structured as follows: in Section 2 the phase-ﬁeld model is constructed based on the chemical potential behavior of the solid propellant combustion model, and the treatment of multiple species is discussed. Furthermore, the derivation of the diﬀuse heat equation is shown and the construction of a stand-in ﬂuid phase is shown. In Section 3 the adaptive mesh reﬁnement and computational solution framework is presented. In Section 4 a qualitative analysis of the model performance is discussed based on four application cases: (i) homogeneoussubstancesregression,(ii)oxidizerandfuelbinderheterogeneoussandwich,(iii)packedoxidizerspheres in a fuel binder matrix, and (iv) Muller packing of oxidizer in a fuel binder matrix. Finally, Section 5 synthesizes the contributions of this work. 2. Diﬀuse interface thermally driven regression model In this section, a full-feedback phase model is developed and applied to the problem of burn interface regression. The model builds on the work originally proposed by Kanagarajan et al.[12], which is a general reference for this section. The theoretical underpinnings of the phase-ﬁeld theory are ﬁrst introduced, and then the diﬀuse boundary method for thermal evolution is presented, along with a surrogate model for estimating the heat ﬂux at the burn interface. 2.1. Free energy functional and regression rate kinetics Followingpriorwork,acontinuumregion ÏRn,n=2;3,asimplyconnectedmeasurablesubset,isestablished containingSCPinboththesolidandﬂuidphases.Theorderparameter, ËH1. ;[0;1]/,isaﬁeldimplicitlytracking the physical state of the material over as a function of space and time. In the present work, material is considered to be either solid (when = 1) or gaseous (when = 0). No intermediate liquid melt phase is considered, though the present model may be extended to include it without adjustment to the method. The boundary between the solid phase and the gaseous phase is consequently tracked implicitly, and the representative sharp interface is taken as the setA= ^xË :.x/ = 0:5`. The diﬀuse boundary is deﬁned as the open set A= ^xË ;0< .x/<1`. The volume of the diﬀuse boundary is given by ðAðùðAð, whereis a regularizing length scale parameter, Meier, Schmidt, Martinez, Quinlan, Runnels Page 4 of 16 Diﬀuse interface method for solid composite propellant ignition and regression and the approximation is exact as long as small compared to the smallest radius of curvature of A. The diﬀuse characterofAcanbeinterpretedintwoways:thediﬀuseregionmayberegardedasaphysicaldomainoverwhichthe solid-to-gas reaction occurs; i.e. interpreting as a reaction coordinate along the lines of the thickened ﬂame model [25, 26]; alternatively, Ais considered strictly as a numerical regularization. The present work adopts this second interpretation, and thus the “error” of the solution, O./, vanishes as 0. This is referred to as the sharp interface limit. The solid phase is a heterogenous composite consisting of an AP oxidizer embedded in an HTPB matrix. The modeling framework is suﬃciently robust to account for other species, such as metal fuels or burn rate enhancers. Here, the simplest SCP composition in order to streamline the formulation of the model is considered. As with the solid/gasinterface(deﬁnedimplicitlyby ),itisadvantageoustodeﬁnethegeometryoftheSCPimplicitlyaswell.A “diﬀusespeciesﬁeld” ËH1. ;[0;1]/isused,inwhichtheregion ^xË :.x/=1; .x/=1`correspondstothat occupied by solid AP; ^xË :.x/=1; .x/=0`corresponds to the region occupied by solid HTPB. (Note that the ﬁeldisnotassignedmeaninginthegaseousregion.)Like ,thediﬀusespeciesﬁeldcontainsdiﬀuseboundaries thatareregularizedtohaveadiﬀusewidth,designated .Thediﬀusivenessof isconsideredtobeapurelynumerical regularizationsincethephysicalboundaryisnearlyperfectlysharp.Unlike ,thediﬀusespeciesﬁeldissimplyassigned rather than evolved, and so the diﬀusiveness constraint is constant and determined at initialization. As with most such models, this phase ﬁeld model is built upon the deﬁnition of a free energy functional f: H1. ;[0;1]/Rthatdrivestheevolutionoftheorderparameter.Here(following[12]),thefollowingorderparameter is used: f[]=Ê w./+1 22ð(ð2 dx; (1) wherewis the chemical potential, can be interpreted as a Lagrange multiplier or non-dimensionalization factor,  is the energy at the solid/gas interface. The governing equation for follows by the adoption of a kinetic relation. In somecases,anonlinearkineticlawmayberequired[27,28].Inmostcasesincludingthis,however,linearkineticsare suﬃcient[29].Therefore,theorderparameterevolutionissetproportionallytothegradientofthefreeenergy,leading to the following kinetic evolution equation d dt=*L df d=*L4  dw d*5 ; (2) whereLis a rate constant referred to as “mobility.” The mobility directly controls the interface velocity, and it is essential that it be computed correctly in order to produce an accurate regression rate. Historically, the regression rate has been connected to pressure via a power-law relation, which correctly captures the pressure dependence of regressioninanumberofdiﬀerentSCPconﬁgurations.Previouswork[12]showedthataspecies-dependentpressure power-law can be used to estimate regression in a phase-ﬁeld model by assigning diﬀerent power-law constants to the fuel binder, the oxidizer, and the interface region. While this heuristic approximation produces reliable results, it ultimatelyfailstocapturethefulleﬀectofthecouplingbetweensolidandgasphases;thisrendersthemodelunreliable for scenarios in which more complex solid/gas interactions are at play. In this updated model, the pressure power law is replaced with an Arrhenius law, L=Aexp0*Ea RuT1 ; (3) whereAis the pre-exponential factor, Eais the activation energy, Ruis the universal gas constant, and Tis the temperature. This form for the mobility is more physically sound, and relies on physical constants that have physical interpretability. Of particular note is the absence of an explicit pressure term from the mobility calculation. Here, explicitpressuredependenceisreplacedbyimplicitdependencethroughthetemperatureevolution,whichisinbetter keeping with the physics of the problem. As it has been demonstrated in [12], the model accuracy has low dependence on the choice of chemical potential and can eﬀectively be ﬁtted as a forth order polynomial function, with prescribed values and zero derivative at =0, ½, and1. To track multiple species, the properties at the interface between materials are computed using a mixture rule on the diﬀuse species parameter (Fig. 2). 2.2. Thermal evolution driven by diﬀuse boundary conditions As interface regression relies on temperature (Eq. (3)), it is essential to accurately evolve temperature within the solid phase as driven by heat ﬂux from the ﬂuid phase across the diﬀuse boundary. Special care must be taken to Meier, Schmidt, Martinez, Quinlan, Runnels Page 5 of 16 Diﬀuse interface method for solid composite propellant ignition and regression Parameter name Symbol Value Base pressure [MPa] pp 1.0 Burned energy density [MPa] w0 0:0 Unburned energy density [MPa] w1 1:0 Activation energy density [MPa] w1_2 2:0 Diﬀuse width [mm] " 0.0005 Chemical potential factor [-]  0.001 Interface energy [MPa]  1.0 Shape constant [-] 0.02726 Table 1 AP and HTPB parameters for the solid phase model imposeaﬂuxboundaryconditionusingadiﬀuseboundary.Onewayistoparameterizesharpinterfaceastheisosurface A=^xË :=½`andthendrivethethermalevolutionaccordingly;thisiscumbersome,however,andultimately counteractstheadvantageofthediﬀuseboundaryformulation.Instead,boundaryconditionsondiﬀuseboundariescan be imposed through source terms obeying the constraint that the sharp interface solution is recovered exactly in the limitas0.Thiswasproposedin[30],anddevelopedmorerecentlyfortheEulerequations[24]andtheequations of linear elasticity [23], and the formulation for the heat equation is presented here. Fourier’s law and the thermal diﬀusion equations for a discrete boundary are given by k(T*q=0 xË unburned (4) c)T )t*( q=0 xË unburned (5) T*T0=0 xË)D unburned (6) n q+q0=0 xË)N unburned (7) where unburnedis the limit of the support of as0,)D unburnedand)N unburnedare the boundaries over which Dirichlet and Neumann conditions are applied respectively, q0is the heat ﬂux, and nis the outward-facing normal. In the diﬀuse boundary context, the explicit boundary conditions are replaced with source terms: k(.T/*q=(kT0 (8) c)T )t*( q=ð(ðq0; (9) whereT0andq0arefunctionsoverthediﬀuseboundary,havingconstantvaluealongthenormaldirectionequaltothe prescribed values in the discrete boundary case. Applying the product rule, and observing that (=*nð(ð, yields k(T*q «­¯­¬ f=ð(ðnkT*T0«¯¬ g(10) 0 c)T )t*( q «­­­­­­¯­­­­­­¬ f1 =ð(ðn q+q0«­­¯­­¬ g: (11) Itwasdemonstratedin[24]that,assumingthequantitiesindicatedas fandgarebounded,andthat (gisboundedin the normal direction, then f=0for all xË unburnedandg=0for all xË) unburnedin the limit as 0. Since fandgin Eqs. (10) and (11) correspond to the governing equation and boundary conditions for Fourier’s law and the heatequation, respectively, Eqs.(8) and(9) are concludedto be equivalentto Eqs. (4)to (7)in the sharpinterface limit. Special care must be taken to ensure boundary condition consistency in the diﬀuse interface formulation as on a discreteboundarydomain.WhenNeumann(heatﬂux)boundaryconditionsaredesired,thenthegoverningequations to be solved are Eqs. (4) and (9). It is worth remarking that one should not attempt to solve Eq. (8) and Eq. (9) Meier, Schmidt, Martinez, Quinlan, Runnels Page 6 of 16 Diﬀuse interface method for solid composite propellant ignition and regression 0.02.55.07.510.012.515.017.520.0Heat Flux - [kW / cm^2] 200 150 100 50 0 50 100 Distance to interface [m] HTPB AP0.1MPa (data) 1MPa (data) 2MPa (data) 6MPa (data) 10MPa (data)0.1MPa (fit) 1MPa (fit) 2MPa (fit) 6MPa (fit) 10MPa (fit) (a)Heatﬂuxfromtheﬂametothematerialsurface compared to [17] η(x) q(x) ϕ(x)(b) BSAMR showing selective meshing for packed spheres in 2D based on temperature and gradients Figure 3: Heat ﬂux and computational domain simultaneously, because this is the same as imposing Neumann and Dirichlet conditions simultaneously, which is unphysical and would lead to numerical error. In the event that one desires to enforce a Dirichlet boundary condition across the diﬀuse interface, one should solve Eqs. (6) and (8). For other cases such as mixed or Robin boundary conditions, there are numerous strategies for consistently applying diﬀuse boundary conditions. In the present work, however, only heat ﬂux conditions are of interest, and so the ﬁnal form of the heat equation to be used here is )T )t=1 c( .k(T/+ð(ðq0 c(12) Equation(12)eﬀectivelyevolvesthetemperatureinthesolidphaseofthesystem,butthermalevolutionisundeﬁned where=0andcaneasilybecomeunstableasonlythetotaltemperatureﬁeldacrossbothphasesiswell-deﬁned.Itis thus necessary to implement an artiﬁcial “stand-in” ﬂuid state to evolve the temperature outside the support of . The stand-inisimplementedhereinthesamemannerasthatproposedin[24];thetotaltemperatureevolvedbyintroducing an auxiliary equation: Ttotal=Tsolid+.1*/Tfluid; (13) whereTfluidis an assigned value for the ﬂuid stand-in and Tsolidis the temperature in the solid phase. The quantity .Tsolid/is evolved according to a forward ﬁnite diﬀerence equation d.Tsolid/ dtóóóót=ðt*1)Tsolid )tóóóót+Tsolidóót*1) )tóóóót; (14) which ensures the stability of Ttotalwithout aﬀecting the solution inside the unburned domain. It is worth noting that whileTfluidis an assigned value in this work, it can easily be replaced by a full ﬂuid-phase solver. 2.3. Surrogate model for heat ﬂux In order to eﬀectively evolve Eq. (12) the surface heat from the ﬂame into the solid ( q0) must be known. There aremultiplepossiblewaystodetermine q0.Onewayistosimultaneouslyrunagasphasecalculationtodeterminethe local heat ﬂux at the boundary. Such a calculation is certainly possible, but is numerically expensive and requires the incorporation of an entirely self-contained gas phase model. This work seeks to use an approach that approximates the ﬁdelity of a full gas-phase model, without the computational expense. By assuming a constant ﬂuid pressure over the course of the regression, the surface heat ﬂux can be approximated as a function of pressure and concentration of species. Here, heat ﬂux data is taken from Gross and Beckstead [31], who reported the surface heat ﬂux during AP/ HTPB regression. By modeling a sandwich structure, it is possible to estimate the heat ﬂux for pure HTPB (Fig. 3a, dashed lines—left), pure AP (Fig. 3a, dashed lines—right), and the spike in heat near the AP/HTPB interface due to theclose-proximitydiﬀusionﬂame(Fig.3a,dashedlines—center).HeatﬂuxesforpureAP,pureHTPB,andthepeak heat ﬂux at the interface are approximated using the dimensionless polynomial ﬁts  qAP.p/=CAP 1p+CAP 0*0 ; (15a) Meier, Schmidt, Martinez, Quinlan, Runnels Page 7 of 16 Diﬀuse interface method for solid composite propellant ignition and regression  qHTPB.p/=CHTPB 1p+CHTPB 0*0 *1; (15b)  qD.p/=CD 2p2+CD 1p+CD 0*1 2 qAP*1 2 qHTPB; (15c) whereCare ﬁtting parameters and 0is the reference length scale for the diﬀuse interface between species. The parameterdetermines the diﬀusiveness of the solid-solid boundary between AP and HTPB. The form of the ﬁt is chosensothatthetotalintegratedheatﬂuxisinvariantwithrespecttothevalueof .Theseareassembledintoageneral approximation function for the heat ﬂux as a function of species composition ( Ë[0;1]), pressurep, and mass ﬂux  m: q0. ;p;  m/=Q qAP + qHTPB.1* /+.0_/. qD/4 .1* /tanhL  m. ;/  mAP 0+.1* / mHTPB 0M ;(16) where mAP 0, mHTPB 0areﬁttingparameterscorrespondingtothereferencemassﬂuxacrossthesolid/gasboundaryinthe computationalsimulation.Themassﬂuxtermapproximatestheheatﬂuxintheunknowndomainbetweentheno-burn case (where heat ﬂux is zero) and the steady case (as reported). The form provides an interpretable yet close ﬁt to the computational data (Fig. 3a) when  m= m0, implying that the burning is occurring evenly and at maximum rate at all regions. The mass ﬂux is determined by,  m. ;/=. /) )t; (17) where. /returnsthedensityofAPorHTPBbasedonalinearmixturerule,andthederivativeof isdeterminedby Eq.(2).Thisclosesthefeedbackloop:theheatﬂuxisdeterminedbytheregressionrate,determinedbytheArrhenius rate law, determined by thermal diﬀusion, determined by heat ﬂux. 3. Computational methods The accuracy of phase ﬁeld models depends on suﬃcient grid resolution across the diﬀuse interface. It has been previously shown that 8 to 16 cells can be used to provide optimal results for SCPs regression modeling [24, 12]. Implementing this level of reﬁnement meshing onto large domains with relatively small interface length scales can become computationally expensive. Adaptive mesh reﬁnement (AMR) is a general class of discretization methods that allows for the construction of a meshing grid with multiple levels of reﬁnement that can be dynamically adjusted based on reﬁnement and CFL criteria. Block-structured AMR (BSAMR) is a particular type of AMR that enhances performance by splitting the reﬁnement levels into independent overlapping grids, using averaging and interpolation methodstoestablishcommunicationbetweendiﬀerentlevelsofmeshing.Thisstreamlinescomputationandpromotes excellent scalability on distributed memory platforms. BSAMR also enables the use of temporal sub-stepping, which in turn allows for the number of reﬁnement levels to be changed arbitrarily without violating the CFL condition. Theworkpresentedhereisbasedonmethodsimplementedintheauthors’in-housesoftware(Alamo[32]),which isbasedontheopen-sourcelibraryAMReX[33].TheuseofBSAMRtotheproblemofSCPregressionisparticularly eﬀective due to the insulating properties of AP and HTPB (as described in [14]) that result in localization of the heat transportaswellastheorderparameterevolution.ItshouldbenotethatBSAMRisusedtoresolveonlytheburnfront, not necessarily the entire geometry of the SCP simultaneously, which would lead to prohibitive cost and minimal improvement in accuracy. To eﬀect reﬁnement, every 10 timesteps, cells are tagged or de-tagged if the reﬁnement criterion is met (or not). Here, the reﬁnement criterion is: ð(ððxðg0:001 or (18) ð(Tðð(xðg5:0 and>0; whereðxðis the local magnitude of the grid spacing. This triggers reﬁnement at the burned/unburned interface (wherethegradientof isnonzero)andwherevertherearelargethermalgradientswithinthesolidphase.Thus,high resolution is maintained at the interfaces and the eﬀects of thermal diﬀusivity are captured at suﬃcient resolution to ensure accuracy (Fig. 3b). Theformulationisfullyexplicit,andsothesolverusesanexplicitforwardEulertimeintegratortoevolveboththe kinetic equation (Eq. (2)) and the heat equation (Eq. (12)). Spatial derivatives of ﬁrst and second order are computed Meier, Schmidt, Martinez, Quinlan, Runnels Page 8 of 16 Diﬀuse interface method for solid composite propellant ignition and regression Parameter name Symbol Combined AP HTPB Thermal Conductivity J msK k 0.4186 0.1463 Density kg m3  1957 920 Speciﬁc Heat J kgK cp 1297.9 2418.29 Thermal Diﬀusivity m2 s 1:610*76:510*8 Pre-exponetial Factor [-] A 1:451051:4103 Activation Energy [K] EA 11000 7500 Base Mass Flux kg m2  m0 1000 5000  qconstant [-] C0 CD 0=0:3225 CAP 0=0.42 CHTPB 0=0.323  qpressure factor [ MPa*1] C1 CD 1=2.797 CAP 1=0.46 CHTPB 1=1.114  qpressure square factor [ MPa*2]C2 CD 2=-0.09906 Ref. interface width [ m] 0 3.0 3.0 3.0 Heat ﬂux scaling factor [ W_m2]Q 107 Table 2 Thermal properties of AP and HTPB and combined diﬀuse interface region. usingcentralﬁnitediﬀerencemethods,andonelayerofghostcellsisusedtoenforcedomainboundaryconditionsand tocommunicatebetweencourse/ﬁnelevels.ForwardEulerintegrationissubjecttoinstabilityunlesstheCFLcondition is met. Here, the CFL condition for a propagation speed, u, is expressed as uóóóót xóóóóf1:0: (19) The computational domain for all cases 4mm1mm.1mm/ for 2D (3D) simulations. VisIt[34]isusedforpostprocessingandvisualization.Allregressionratesarecalculatedbydeterminingthefarthest progression of the = 0:5isosurface. This mimics the experimental “trip wire” method, that also bases regression rates on the farthest progression of the interface (rather than, for instance, the average location of the interface). 4. Results ThissectionpresentsresultsfortheapplicationofthemodeltoAPmonopropellantandAP/HTPBsolidcomposite propellants.ItshouldbenotedthattheformofthemodelisgeneralandnotuniquetotheseparticularSCPspecies.In thiswork,however,theexamplesarelimitedtoAP/HTPB.Allparametersarekeptthesameacrossalltests;thenumber ofcalibratedparametersareminimalandlimitedmainlytotheexponentialpre-exponentialfactorandactivationenergy. Physical properties are taken from known values. To determine the values of the ﬁt parameters for AP, a one- dimensional calibration is performed. The pressure range for which the model is considered to be valid is 2MPa to 6MPa, as this is the range in which AP can sustain self-deﬂagration for those pressure levels (and for which experimental data is available). Calibration was accomplished through a grid search in parameter space and the judicious use of the bisection method to identify correct parameter values. In total, approximately 700 simulations were performed in order to evaluate best-ﬁt values for the z0and m0over the valid pressure range. General model parameters are presented in Table 1; parameters speciﬁc to particular species are presented in Table 2. 4.1. Thermal transport and deﬂagration of pure AP A one-dimensional test of heat transfer in AP is considered to validate the diﬀuse boundary method for heat-ﬂux driventemperatureevolution.ThedomainisahomogeneousAPsample,andtheinterfaceisdiﬀusedbutnotpermitted to regress (i.e., the pre-exponential factor is set to zero). The system is exposed to a heat ﬂux of 1106W_m2, and evolved by Eq. (12) over 0.5seconds with AMR active. The thermal evolution is measured over time, and compared to the sharp-interface solution for transient conduction described in [35] (Fig. 4a). Very close agreement is observed between the diﬀuse boundary solution and the sharp interface solution. The main diﬀerence is a slight delay in the Meier, Schmidt, Martinez, Quinlan, Runnels Page 9 of 16 Diﬀuse interface method for solid composite propellant ignition and regression 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 Distance [m]300350400450500550600T emperature [K] 0.5s (Exact) 0.4s (Exact) 0.3s (Exact) 0.2s (Exact) 0.1s (Exact) 0.01s (Exact)0.5s (Alamo) 0.4s (Alamo) 0.3s (Alamo) 0.2s (Alamo) 0.1s (Alamo) 0.01s (Alamo) (a) Comparison between Alamo results (markers) and the exactsolution(solidlines)forthetemperatureevolutionfor a 1D pure AP domain under constant heat ﬂux. 0.0000 0.0001 0.0002 0.0003 0.0004 Distance [m]300400500600700800T emperature [K]5 ms 10 ms 15 ms 20 ms 25 ms 30 ms 35 ms 40 ms 45 ms 50 ms(b)Evolvingtemperatureproﬁleover50msforanexternalpressure of 4MPa and stand-in ﬂuid temperature of 600K. The dashed line indicatesthetemperatureproﬁleintheﬁctitiousgas-phasestand-in. Figure 4: Veriﬁcation and demonstration of the diﬀuse boundary heat transport model diﬀuse boundary solution, causing the temperature to be slightly lower. This can be attributed to a time shift that results from the diﬀusiveness of the interface: the initial heat ﬂux is slightly limited and this causes the evolution to slightlylagtheexactproﬁle.Althoughtheabsoluteerrorgrowswiththeincreasingtime/solutionvalue,the%errorin the solution decreases. Furthermore, apparently increased diﬀusiveness of the boundary also decreases the % error. Next,theevolutionoftemperaturedrivenbyheatﬂuxwithamovingboundaryisconsideredinordertodetermine thesteady-statebehaviorofthemodelforthemassﬂux heatﬂux thermalevolution regressionrate massﬂux cycle.Unlikepreviousphaseﬁeldmodelsforregression,thepresentmodelisstableintheunburnedconﬁguration,i.e., becausethereis nomassﬂux,there isnoheatﬂux.Inthecurrent simulations,itisnecessary toinitiatetheregression of the solid through an imposed heat ﬂux, mimicking the physical process of ignition. Starting from a equilibrium, room-temperature solid, the sample is subjected to a heat ﬂux of 1106W_m2, which initiates a rise in temperature and subsequent regression of the interface (Fig. 4b) Thermal behavior of the solid propellant is as expected, with the temperatureattheinterfacesimilartothatdescribedintheliterature,stayingwithintherangeof500Kto900K[3].The thermalproﬁleinthesolidregion(Fig.4b—solidlines)isconsistentwiththethermalconductivityofAP,andreaches asteadystateafterabout10ms.Thethermalproﬁleisincludedforthestand-ingaseousphaseaswell(Fig.4b—dashed lines) to illustrate the eﬀect of the stand-in: without the stand-in, the thermal evolution in the burned phase generally becomesunstable.Ifthismodeliscoupledtoafullyreactingﬂuidsolver,theﬂuidstand-inwouldbereplacedwiththe simulated ﬂuid temperature in the =0domain. Regression rates were measured for pure AP in the 1–6MPa pressure range and compared to experimental data from Price et al.[36] (Fig. 5a). The correlation between the model and experimental data is good. Interestingly, the regression rate was relatively insensitive to the ﬁtted model parameters. For instance, small changes in the Arrhenius ratelawtendedtohavelittleeﬀectontheregression,whichappearstobeduetothefactthatchangesintheregression rateinturnchangesthetemperatureaswell.Asubstantialimprovementinthepresentmodelisitsabilitytoapproximate thepressurecutoﬀforAPdeﬂagration.althoughthemodelpredictsasmallburnrateat1MPa,itisclosetotheobserved behavior. This can be attributed to the physically realistic feedback between the mass ﬂux and heat ﬂux in the model. 4.2. Burning of AP/HTPB sandwich composite The “sandwich” conﬁguration is useful for experimentally or numerically investigating deﬂagration behavior in a simple yet nontrivial conﬁguration. In the present work, the sandwich structure provides a useful test case on which to further calibrate the model. Here, a laminate of HTPB with thickness equal to 110*4meters is placed in an AP matrix. All of the parameters for AP in the previous section are left unchanged; only the HTPB parameters are varied. Therefore, the only variables calibrated using this case are the length scale of the diﬀuse interface between AP and HTPB 0, and base mass ﬂux for mixed substance  m0. The size of the interface plays an important role in Meier, Schmidt, Martinez, Quinlan, Runnels Page 10 of 16 Diﬀuse interface method for solid composite propellant ignition and regression 1 2 3 4 5 6 Pressure [MPa]012345678Burn Rate [mm/s] Pure AP - Experiment Pure AP - Model (a) Pure AP 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Pressure [MPa]2.02.53.03.54.04.55.05.5Burn Rate [mm/s] Sandwich - Experiment Sandwich - Model (b) AP/HTPB Sandwich Figure 5: Alamo regression rate prediction. (a) Comparison of experimental data for pure AP with Alamo’s calculated average burn rate for the pressure range of 1 to 6MPa. (b) Comparison of Alamo results and experimental data for a sandwich set up for pressure ranges of 0.8MPa through 4.0MPa. Laser IgnitionPure AP (a) (b) (c) (d) Figure 6: The ignition energy (laser-induced heat ﬂux) is applied to a middle square in a pure AP domain for the ﬁrst 0.025seconds and then shut oﬀ (a,b, and c). Once the ignition temperature is achieved the burn can self-propagate through the surface and regress throughout the domain (d). assuring convergence of the results, since it controls the model’s ability to capture the eﬀects of higher heat ﬂuxes at the AP/HTPB interface. A minimum of eight interface cells is required for proper evaluation of those eﬀects. Regression rates were calculated for the AP/HTPB sandwich in the 0.8–4MPa pressure range, and compared to correspondingexperimentaldatafrom[36](Fig.5b).SincephysicalpropertiesattheAP/HTPBinterfacearecomputed usingthelinearmixturerule,theburningrateatthematerialinterfacevariesbasedontheburnrateofpureHTPB.Itis observedthattheregressionratesforthesandwichcaseshowagoodmatchbetweentheexperimentsandsimulations. 4.3. Ignition in pure AP The model’s capability is further demonstrated with two sets of example simulations: (i) The model’s ability to sustain deﬂagration, and (ii) the eﬀects of laser power on the ignition process. In the ﬁrst one, the model is applied in a 3-dimensional domain with the laser turned on for the initial 0.01seconds of integration over a smaller region of the body’s left face (0.1 by 0.1mm square face). The burn front propagation is simulated over a 0.1-second interval (Fig. 6). It is observed that the burning starts locally while the laser is on but the system can self-sustain after it is shut oﬀ att=0:025s. It can also be observed that areas outside of the initial laser spot eventually receive lateral heat through the mass/heat ﬂux feedback cycle and conduction, causing it the regression interface to propagate along the body’s left face. An analysis of the eﬀect of laser power on the system is also performed. The laser is set as a pulse that starts addingenergytothesystemat0.0secondsandstopsat0.025seconds,withatotalintegrationtimeof0.05seconds.It is observed that the AP would sustain self-deﬂagration if the total energy added to the system reaches the activation energy. For the pure AP case, the total energy added is enough to raise the surface temperature to 400K, and the minimum laser energy to reach that temperature under 0.025seconds was found to be 0:9106W_m2. If the laser powerisbelowthatvalue,thetemperaturepeaksat0.025secondsandthencools.Thisresultisparticularlyinteresting Meier, Schmidt, Martinez, Quinlan, Runnels Page 11 of 16 Diﬀuse interface method for solid composite propellant ignition and regression 0.00 0.01 0.02 0.03 0.04 0.05 Time [s] 0 1 2 3 4 5 6 7Burn Rate [mm/s] mmmmm 5.0e7 /2 1.0e7 /2 5.0e6 /2 1.0e6 /2 0.9e6 /2 0.8e6 /m2W W W W W W (a) Burn rate variation with laser power 1022×1023×1024×1026×102 Heat Flux [cal/s/cm2]103 6×104 2×103 3×103 4×103 Square Root of Ignition Time [s1/2]Baer Eq. (2 MPa) Baer Eq. (3 MPa) Baer Eq. (4 MPa) Baer Eq. (6 MPa) Model (2 MPa) Model (3 MPa) Model (4 MPa) Model (6 MPa) (b) Ignition time and dependence on laser heat ﬂux Figure 7: (a) The laser is turned on at time t=0seconds and turned oﬀ at time t=0:025seconds. If the laser power is high, the burn rate increases. If the laser energy is too low, no ignition occurs. If ignition occurs, the burn proceeds at a uniform rate after laser is deactivated even if over driven by a high input power. (b) Comparison between the simulations and the values from literature [37] indicates realistic ignition times and scaling. sincesimilarhavepreviouslybeenreported[10],andnotuningofthemodelwasmadetomatchthisobservedignition temperature.Moreover,itisobservedthattheburnratedirectlydependsontheamountofenergyaddedtotheinterface as it regresses (Fig. 7a). It is also apparent that after the laser is shut oﬀ, the burn rate of all simulations that reached the ignition energy approach the self-sustained burn rate. Oscillations in the eﬀective burn rate can be observed as a consequenceoftheslowpropagationofenergythroughconductionandquicklocaltemperatureincreaseastheevolving interface regresses. Evaluating the solid propellant response to diﬀerent levels of ignition energy is crucial for allowing the proper design of igniters. Baer and Ryan have developed correlations between ignition time and ignition energy based on variables such as surface heat ﬂux, pressure, and material composition [37], which relates the ignition time and the ignition energy ﬂux to ignition time by assuming a homogeneous, semi-inﬁnite body with uniform initial temperature and constant heat ﬂux (see Equation 4 in [37]). Since the proposed assumptions are natural conditions for a computational simulation, simulations are performed for pure AP cases in the pressure range of 2MPa through 6MPa for six diﬀerent laser power levels, for a total of 24 samples. Lower laser power was tested but did not initiate regression. Baer’s equation is not directly dependent on pressure but its eﬀects are captured by the surface reaction frequency factor and the transmitted energy per unit reaction, which is equivalent to the pre-exponential factor in the Arrheniuslaw.SimulationresultsarecomparedtothesolutionofBaer’sequationrelatingignitiontimetoheatﬂux,and areasonableagreementisobserved(Fig.7b).Again,itshouldbenotedthatthemodelusedthecalibrationcoeﬃcients found in the previous section to produce these results. 4.4. Burning of packed AP spheres in HTPB matrix Akeyadvantageofthediﬀuseinterfaceapproachisitsabilitytoaccountforvariouspropellantcompositionswith arbitrary morphology. To demonstrate this, a number of 3-dimensional simulations were run to determine regression ratesinanSCPconsistingofAPspherespackedintoaHTPBmatrix.Thepackedsphereconﬁgurationsweregenerated usingaMillerpackingscheme,asdescribedin[38,40],whichconsistsofinsertingspheresinthedomainwithaninitial low volume fraction and then resizing the spheres to achieve the target mass fraction. To ensure accurate comparison to experimental data, 4 diﬀerent 3-dimensional packings have been created with mass fractions between 60% and 80% and total AP spheres of 500 and 1000, for a total of 40 diﬀerent packings (Fig. 8). To generate 2-dimensional packings,the103-dimensionalsamplesareslicedinto10cross-sections.Sliceswitha2-dimensionalpackingfraction less than 40% were removed, leaving a total of 80 2-dimensional samples. Simulations of ignition and burn for 2- dimensional and 3-dimensional domains are run for a simulation time of 0.05seconds. Ignition was eﬀected with a laserpowerof 1:0106W_m2appliedtotheboundaryduringtheinitial0.025seconds.Thetimesteponthecoarsest levelofthesimulationissetto 510*6s.Thesamplesweresimulatedfor4diﬀerentpressures,resultinginatotalof Meier, Schmidt, Martinez, Quinlan, Runnels Page 12 of 16 Diﬀuse interface method for solid composite propellant ignition and regression (a) (b) 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 Sphere Radius [m]050100150200250300350Number of Spheres1000 Spheres 500 Spheres (c) 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 Sphere Radius [m]050100150200250Number of Spheres1000 Spheres 500 Spheres (d) Figure 8: 3D samples of Muller Packed AP spheres [38] and particle distribution. (a) 60% AP mass fraction. (b) 83% AP mass fraction. (c) AP particle size histogram for 60% AP mass fraction. (d) AP particle size histogram for 83% AP mass fraction. 320 2-dimensional simulations and 160 3-dimensional simulations. The burn rate represents the average value of all instantaneous burn rates between every time step and the error bar shows its standard deviation. The regression rate results, along with error bars indicating standard deviation in instantaneous burn rate, are comparedwithexperimentaldata[39](Fig.9).Asigniﬁcantvariationintheaverageburnrateacrossdiﬀerentsamples withthesamepressureandcompositionisobserved.Thisisexpectedsincetheburnrateishighlyinﬂuencedbyparticle distribution and particle size. Although all samples had relatively similar total particle size distribution, the particle size at the region at the section regressed during 0.05seconds have a higher variance, leading to the higher impact on regressingbehavior.Thetemperaturewasdeterminedtobeintherangeof700Kthrough1200K,whichiswithinthe expected range. The burn rate dependency on physical characteristics of SCPs, such as mass fraction and particle size, are well reported, and the simulations in this work have been able to mimic those behaviors. For instance, as the AP mass fractiondecreasesbelow55%theHTPBvolumefractionishighenoughthatdeﬂagrationceasesandonlyaveryslow regression is observed. This is often observed in low-pressure simulations where the thermal and chemical forces are not strong enough to initiate the regression if no AP particles are present in the front interface. If this pure HTPB interface is located in the middle of the domain, the result is a high increase in the local temperature with minimal or no interface movement. In such cases, the local temperature at the interface propagates very quickly, which can cause CFL condition issues. Robustness requires careful setup of integration time or an adaptive time step based on local properties.Particlesizeimpactonburningratehasalsobeendemonstratedinthesesimulations,decreasingtheaverage particle size reﬂected in an increase in burn rate. This phenomenon is widely reported in the literature [13, 41, 42] as the increase in the surface area between AP and HTPB accelerates the regression rate due to the close proximity of Meier, Schmidt, Martinez, Quinlan, Runnels Page 13 of 16 Diﬀuse interface method for solid composite propellant ignition and regression 1 2 3 4 5 Pressure [MPa]1 0123456Burn Rate [mm/s] Kohga (MF 60%) Kohga (MF 70%) Kohga (MF 80%) 60±2.5 % 65±2.5 % 70±2.5 % 75±2.5 % 80±2.5 % (a) 1 2 3 4 5 Pressure [MPa]0123456Burn Rate [mm/s]Kohga (MF 60%) Kohga (MF 70%) Kohga (MF 80%) Alamo (MF 60%) Alamo (MF 65%) Alamo (MF 70%) Alamo (MF 75%) Alamo (MF 80%) (b) Figure 9: Comparison of Alamo burn rate with experimental data from [39] for a packed spheres set up. (a) Average and standard deviation for 320 samples of 2D simulations with mass fractions in the range of 55% through 85%; (b) Average Burn rate for 40 Samples of 3D simulations with AP mass fractions of 60% through 80%; Spheres represent pickings with 1000 AP particles and squares are packings of 500 AP particles. the premixed ﬂame. This is observed to be true for most cases, but large AP particles in low-pressure conﬁgurations increase the burn rate since the self-deﬂagration of AP has a faster burn in those cases. Overall, the results show reasonable accuracy when compared to experimental data, especially considering the minimal number of ﬁtted parameters. The model prediction is more accurate for larger AP mass fractions, and less accurate for more sparsely packed compositions. This is a result of the reduced order heat ﬂux estimation that is based only on local composition, and does not take into account the nonlocal characteristics of the burn interface. For the higher AP loadings, the interface is more planar, and therefore more geometrically similar to the sandwich conﬁgurationfromwhichtheheatﬂuxdataismeasured.Forlowerloadings,ontheotherhand,theinterfacecomplexity is substantially increased, which alters the amount of heat transferred into the solid. 5. Conclusion In this work, a computational model is presented that predicts regression rate and interface morphology using a diﬀuse interface model, This work builds on the previously developed phase ﬁeld model for regression, but replaces theheuristicpressure-dependencewithaphysically-motivatedArrheniuslawinformedbylocaltemperature.Theheat equation is formulated, using the diﬀuse boundary method, to impose heat ﬂux boundary conditions using an order parameter-basedsource-term.Thisiscoupledwithastand-inmodeloftheﬂuidphaseinordertoregularizethesolution and prevent numeric overﬂow. The local temperature determines the local mobility of the regression surface, and is now directly linked to physical properties and can account for both thermal and chemical kinematic eﬀects. The thermal behavior of the burning of solid composite propellants is successfully modeled and the temperature proﬁle hasgreatagreementwithexperimentaltestingofthosematerials.RegressionrateswerecalibratedusingpureAPand AP/HTPBsandwichcases.Basedonthosecalibrations2-and3-dimensionalsimulationsofpackedspherescaseswere performed and the regression behavior has shown reasonable qualitative results that followed experimental data for diﬀerent pressure and material compositions. It was also demonstrated that the model accurately captures ignition behavior of AP when subjected to thermal loading. The primary limitation of the present model is its reliance upon a reduced-order estimate of heat ﬂux from computational sandwich tests. This simpliﬁcation is responsible for the error in the predicted regression rate, which is consistent with the expected correlation between regression surface complexity and error. To obtain heat ﬂux information that is accurate for complex regression surfaces, the gas phase must be resolved simultaneously. Such an implementation is outside the bounds of the present work; however (as previously stated) the two phases can be readilycoupledacrossadiﬀuseboundary[24,23].Suchacoupledsolverconstitutesfuturework;however,thepresent results justify the use of the model for analysis of SCP regression and morphological dependence. Meier, Schmidt, Martinez, Quinlan, Runnels Page 14 of 16 Diﬀuse interface method for solid composite propellant ignition and regression (a) 0.01 seconds (b) 0.02 seconds (c) 0.1 seconds (d) 0.2 seconds Figure 10: Evolution of a 65% AP mass fraction samples over the course of 0.2seconds. Color corresponds to surface temperature. 6. Acknowledgements AuthorsMM,ES,PM,andBRacknowledgesupportfromtheOﬃceofNavalResearch,grant#N00014-21-1-2113. PMacknowledgessupportfromtheCaltechSummerUndergraduateFellowship(SURF)Program.Thisworkusedthe INCLINE cluster at the University of Colorado Colorado Springs. INCLINE is supported by the National Science Foundation, grant #2017917.