Solid composite propellants (SCPs) are ubiquitous in the ﬁeld of propu lsion. In order to design and control solid rocket motors, it is critical to understand and accurately pred ict SCP regression. Regression of the burn surface is a complex process resulting from thermo-chemical -mechanical interactions, often exhibiting extreme morphological changes and topological transitions. Diﬀuse interfac e methods, such as phase ﬁeld (PF), are well-suited for modeling processes of this type, and oﬀe r some distinct numerical advantages over their sharp-interface counterparts. They also provide a convenie nt framework for incorporating multiple multiphysical dynamics. In this work, we present a phase-ﬁeld me thod for modeling the regression of SCPs with varying species and geometry. We construct the model from a the rmodynamic perspective, leaving the base formulation general. A diﬀuse-species-interface ﬁeld is e mployed as a mechanism for capturing complex burn chemistry in a reduced-order fashion, making it possi ble to model regression from the solid phase only. The computational implementation, which uses block-str uctured adaptive mesh reﬁnement and temporal substepping for increased performance, is brieﬂy discus sed. The model is then applied to four test cases: (i) pure AP monopropellant, (ii) AP/PBAN sandwich, (iii) AP/HTPB san dwich, and (iv) spherical AP particles packed in HTPB matrix in two and three dimensions. In all c ases, reasonable quantitative agreement is observed, even when the model is applied predictive ly (i.e., no parameter adjustment), as in the case of (iv). The validation of the proposed PF model demonstrates it s eﬃcacy as a numerical design tool for future SCP investigation. Keywords: Solid composite propellants, Numerical simulation, Phase ﬁeld modeli ng, Diﬀuse interface methods 1. Introduction Solid composite propellants (SCPs) are composed of a mixture of fuel and o xidizer species that are unmixed at the molecular level. An example of a modern SCP oxidizer/f uel combination is ammonium perchlorate (AP) particles in a hydroxyl-terminated polybutadien e (HTPB) matrix. SCPs such as this have become the propellant of choice in the ﬁeld of solid rocket propulsion f or tactical missiles and launch vehicle ∗Corresponding author Email address: brunnels@uccs.edu (Brandon Runnels) 1 boosters. Even though SCPs have slightly lower performance than com plex liquid rocket engines, they are simple, reliable and chemically/mechanically stable in long-term st orage. In solid rocket motors the burn rate and its dependence on the gas dynami c ballistics is a critical variable for a number of reasons. First, the burn rate and grain surface area deter mine the propellant ﬂow rate and, as a result, the motor thrust. Second, how the burn rate is changed by t he local pressure and ﬂow near the burning surface is an important factor in preventing combustion inst abilities. The burn rate is a function of the properties of the solid components and their mass fractions, enthal pies of phase change, chemical reaction mechanism, ﬂame structures, heat transfer mechanisms, and intern al ballistic gas dynamics. Speciﬁcally, the most important parameters are the AP mass fraction and the AP particle si zes [1]. An extensive set of experiments varying the particle size and concentration of AP on the bu rning rate was performed by Bastress [2]. While there have been numerous experimental measurements of SCP burn rates, many of the details of the combustion and detailed solid-ﬂuid interaction remain unknown. One of the confounding issues is the complex combustion mechanism. I n addition to being an oxidizer, pure AP can support a self deﬂagration wave, for which the burn rates hav e been measured [3–7]. On the other hand, the fuel/binders cannot act as a monopropellant, rather, h eat transfer from the gas-phase reactions causes the long polymer chains of the binder to pyrolyze and s ublimate into mainly shorter chain hydrocarbons. When AP particles are suspended in a matrix of binder s uch as HTPB or polybutadiene acrylonitrile (PBAN) the resulting ﬂame structure is complex and un steady, which makes the experimental burnratemeasurementchallenging. Despitethis, asubstantialamoun tofofhistoricaldatahasbeengathered over the past several decades [8, 9]. The two most common methods for meas uring the burn rate are with break-wires [10] or optically [9, 11–14]. Experimental investigation of sandwi ch structures has made it possible to distill the binder/oxidizer interactions in a simpl e yet nontrivial SCP geometry. For instance, the burn rates have been measured in a composite structure made of a s ingle layer of binder sandwiched between two layers of pure AP [4, 15, 16]. These data can be used to validate numerical burn rate models. There are several numerical models that have been developed to est imate the burn rate. A simple steady-state combustion of SCPs was developed as a two parameter, power law in pressure r=a(P/P0)n [17]. Even though this model fails to capture important features such as the relation between the unsteady ﬂame structure and the burning rate, it remains very useful. Anothe r model was developed that introduces a mechanism to drive the energy release at the burning surface and the heterogeneous reaction of the oxidizer [18]. The next increase in model complexity comes by considering th e ﬂame structure. For instance, models have been introduced where the assumed structure is three ﬂame s instead of one: (a) premixed binder rich ﬂame, (b) lean premixed oxidizer (AP) ﬂame, and (c) the diﬀusi on ﬂame anchored by the other two premixed ﬂames [19–21]. The relation between ﬂame structure and burn rat e calculations in modeling has been extensively studied [22–24]. In some gas-phase models, erosive burning near the surface of the propellant have been incorporated into a three-ﬂame model to inclu de the eﬀect of the high-velocity gas ﬂow [25, 26]. Since real SCPs are have a wide distribution of AP particle s izes, models have been created to 2 include a bimodal distribution of AP particles, for instance, the pe tite ensemble model, which also includes the temperature sensitivity of the SCP [27]. An analytical method base d on pressure and temperature correlation has been developed to study the dependence of SCP deﬂagrat ion rate on pressure and initial temperature [28] for SCPs. One of the key ﬁndings experimentally was t hat the burning rates of the AP composite propellants increase with increasing AP concentration and de creasing AP diameter [29]—a ﬁnding that is supported with current modeling eﬀorts [30]. Experimental observation indicates that the burn surface is far from p lanar [31]. In general, material heterogeneity produces a highly complex burn front, sometimes pro ducing particles of polymer binder or oxidizer that are advected away by the ﬂow [32]. To model this behav ior, it is necessary to adopt a numer- ical approach capable of resolving complex geometry and topological transitions . Explicit meshing of the regressing surface is cumbersome and can be prohibitively diﬃcult , especially when topological transitions occur. Eulerian approaches, in which the interface moves through th e grid, are more attractive. Two such methods, both treating the interface as a sharp discontinuity, hav e been predominantly used to study regres- sion. The ﬁrst, originally proposed by Hegeb et al. [33] and since applied by Massa et al., [34] and others [35], is to use a coordinate transformation to map the regressing surface to a planar surface. Although this works well for near-planar surfaces, it can fail when the surface exper iences a topological transition, which can make it impossible to construct a smooth map. The second is to use t he level-set (LS) method, which relies on the evolution of a surrogate ﬁeld, an isocontours of which may be regarded as a sharp surface. This was originally applied by Wang et al. [36] and others [37, 38] to modeling surfac e regression in SCPs. Here, surface regression is eﬀected by the calculated local regression rate , which is used to update the surrogate ﬁeld and then a new interface is approximated by the updated isocont our. While sharp interface methods are often desirable due to their increased accuracy, we suggest that th e problem of interface regression does not particularly beneﬁt from this enhanced precision due to the lack of sharpness in the physical interface itself. As such, some of the cumbersome numerical tasks needed for sh arp interfaces (front tracking, explicit boundary conditions, nondiﬀerentiability, etc.) needlessly comp licate the solution. Inthiswork,wepresentadiﬀuseinterfacemodelforsurfaceregre ssion. Herethesolid/gasphasetransition is treated as continuous at very small scales. The phase ﬁeld (PF) mode l provides an attractive framework for modeling combustion [39–41] using a Gibbs-type construction. It als o provides an elegant means for tracking complex surface evolution by modeling interfaces as diﬀus e. PF models have been employed in a wide range of settings including microstructure evolution [42–44], frac ture mechanics [45, 46], dendritic growth [47], and solid-ﬂuid interactions [48–50]. PF provides a robust mean s for thermodynamic modeling of combined multiphysics processes, while retaining the ability to capture complex morphological behavior. It also provides interpretability to model parameters (e.g., order p arameter as a reaction coordinate, chemical potential as a reaction energy, diﬀuse boundary as a reaction zone, etc.). Indeed, the great advantage of the phase ﬁeld method is the ease with which it can be extended to treat fully complex ﬂuid phase phenomena (mixing, diﬀusion, and reactions) as well as solid phase behavior (the rmal and elastic response). In this 3 work, the scope is restricted to burn rate prediction based on ambien t pressure only using a heuristic rate law, and does not include concurrent simulation of ﬂuid or solid phase m echanics. It is entirely possible, though, to couple this present model to both solid and ﬂuid phase mode ls, replacing the heuristic rate law with fully resolved solid and gas dynamics. (We refer the reader to con current work [51], which directly addresses the challenges of solid/ﬂuid coupling in diﬀuse boundary mo dels.) The PF can caues a substantial increase in the computational resolution n ecessary to resolve the solution acrossthediﬀuseboundary. Throughthestrategicapplicationofmeshing techniqueslikelocaladaptivemesh reﬁnement, however, it is possible to achieve good performance wi th suﬃciently high accuracy. The remainder of the paper is structured as follows: In Section 2 we d evelop the chemical and thermo- dynamic description of the phase ﬁeld solid propellant burn model, t he derivation of the burn speed from the model, and the incorporation of multi-species eﬀects though a diﬀ use species parameter. The presented model, like most diﬀuse interface methods, requires the applic ation of specialized computational techniques to solve in a tractable manner. Therefore, in Section 3, we brieﬂy dis cuss the computational framework that is used to implement the model and generate results. In Section 4, we examine four distinct cases: pure AP, AP with a thin “sandwich” laminate of PBAN, AP with a thin “sandwich” lami nate of HTPB, and packed AP spheres in an HTPB matrix. The ﬁrst three cases are used to cal ibrate our model, which is then applied to the practical example of packed AP spheres in an HTPB matrix in two and three dimensions. In all of these cases, we demonstrate the eﬃcacy of the model via comparison to legacy experimental data. We conclude this work in Section 5 with a discussion synthesizing re sults and the limitations of the model. This section also outlines the path forward for integrating other physical dy namics into the proposed framework. 2. Phase ﬁeld solid burn model In this section we present the phase ﬁeld model as applied to the pr oblem of SCPs. Consider a continuum regionΩﬁlledwithburnedorunburnedSCP.Ωisgenerallyunintere sting, andisinpracticeoftenarectilinear computational domain, as the geometry of the solid phase is resolved implic itly (Fig. 1, left and center). Over this region we deﬁne a ﬁeld variable, commonly referred to as an order par ameterη: Ω→[0,1]. Physically, ηmay be regarded as a pointwise reaction coordinate that tracks the progres s of the burn reaction as a function of space and time. In this work we regard the problem as being tw o-phase, considering solid and ﬂuid only; that is, no melting is considered. It is important to note that in a case with more species or phases,ηcan no longer be interpreted as a reaction coordinate and should be replac ed with an alternative construction such as a set of volume or mole fractions. The order parameter functions simultaneously as a thermodynamic reac tion coordinate and an indicator function describing the geometry of the unburned region. We now begi n our construction of a phase ﬁeld model through the establishment of a free energy functional f:C2(Ω,[0,1])→R, given by, f[η] =/integraldisplay Ω/bracketleftBig λw(η)+1 2ε2κ|∇η|2/bracketrightBig dx, (1) 4 1 0.5 Ea ΔHr burned regionunburned regioninterface region burnedunburned Figure 1: The order parameter ηrepresentations—(a) Order parameter η= 1 shows unburned region, η= 0 shows the burned region and the interface region is 0 < η <1, (b) The representation of order parameter ηin terms of burned, unburned and interface regions along the length of the coupon, (c) The rep resentation of free energy w(η) in terms of burned, unburned and interface regions wherethesquarebracketson fshouldbeunderstoodtoindicatedependenciesontheenclosedargu mentsand their spatial/temporal derivatives. Normal invariant notation ( ∇indicates the gradient, ∆ the Laplacian, etc.) is employed. The ﬁrst term in the integrand indicates the c hemical potential, w, which is generally a diﬀerentiable multiwell function with minima at 0 and 1 (Fig. 1, righ t). The coeﬃcient λis, for practical purposes, a scaling factor, but also admits interpretation as a Lagrange mul tiplier that may be scaled arbitrarilytobringthesolutionclosetothesharpinterfacelimit . Thesecondtermintheintegrandintroduces nonlocality by penalizing gradients in η. The coeﬃcient κcan be interpreted as the energy of the interface (i.e., surface tension) between the burned and unburned region; h ere, we simply view it as a numerical regularization. The variable ε, then, is a parameter that controls the length scale of the solution, and w ill be discussed in more detail subsequently. The free energy functional fprovides the basis for the construction of a kinetic evolution equati on. We adopt the classical ansatz, following non-equilibrium thermodynamic theory as well as the cl assical Ginzburg- Landau theory [52] that the order parameter evolves proportionally to the gradi ent of the free energy. A kinetic variable, called the “mobility” ( L) is used scale the ﬂow of η. As will be shown subsequently, dividing by the length scale εis necessary to prevent dependence of the ﬂame speed on the numeri cal interface region width. The equation governing the evolution is given by the followin g: ∂η ∂t=−L εδf δη=−L/bracketleftbiggλ ε∂w ∂η−εκ∆η/bracketrightbigg , (2) whereδ/δηis the variational derivative. The relatively low order of Eq. (2) makes it readily solvable using optimized explicit methods, as discussed in Section 3. Coupling t o additional physics occurs through the modiﬁcation of Eq. (1) to include other thermodynamic contributions such as thermal or elastic free energies, which are then integrated into the governing equation through the vari ational derivative. As we will show, however, the reactive and diﬀusive components alone suﬃce to captur e a broad range of SCP behavior. 2.1. Determining burn rate in the limit Of particular interest in this present work is the connection of the governing equation to the burn rate in the sharp-interface limit. To begin, we derive the relationship b etween one-dimensional burn rate and model 5 1 0 1 2 3 4 5 x 0.0 0.2 0.4 0.6 0.8 1.0g(x) = 1.0 = 0.5 = 0.1 = 0.05 = 0.01 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 = (x ct )/ 0.0 0.2 0.4 0.6 0.8 1.0g( )Figure 2: Illustration of 1D burn rate calculation. Diﬀuse i nterface converges to the sharp interface limit as ε→0, and the interface moves from left to right at a constant velocity c(left). All functions can be collapsed onto a single functio ng(ξ), where the derivative is nonzero only in the shaded ( −ǫ/2,ǫ/2) region (right). parameters in a single-species material, in order to interpret the model parameters ( L,ε, etc) in terms of measurable properties. (We comment that the proof presented in thi s section, though important, is fairly technical. If desired, the reader may skip to the salient result , Eq. (9), without a loss of continuity.) For the case of steady burning, we approximate ηε(x,t) =g(ξ) ξ=x−ct ε, (3) wheregis a monotonically increasing C2function satisfying: (1) g(ξ) = 0,g′(ξ) = 0 when ξ≤ −ε/2; (2) g(ξ) = 1,g′(ξ) = 0 when ξ≥ε/2; (3) antisymmetry of g(ξ)−1 2aboutξ= 0 (Fig. 2). The constant cis the rate of regression, which will be solved for here. We do not currently specify any limits on ε, but will show that the solution satisﬁes the governing equations and scales with ε. Now, substituting into the governing equation produces −c εdg dξ=−L/bracketleftBigλ εdw dg−κ εd2g dx2/bracketrightBig . (4) The nonlinearity induced by the chemical potential renders Eq. (4 ) unsolvable in the general case, which is undesirable as we seek a general relationship between the chemic al potential, mobility, and burn width. Therefore we seek a solution in the weak form by integrating both sides over the region [ −ε/2,ε/2]: c/integraldisplayε/2 −ε/2dg dξdξ+Lκ/integraldisplayε/2 −ε/2d2g dx2dξ=L/integraldisplayε/2 −ε/2dw dgdξ. (5) 6 By inspection, and application of the properties of g, we ﬁnd that the ﬁrst integral reduces to unity, and the second to zero. Next, applying a change of basis, the integral on the right- hand side becomes /integraldisplay1 0dw dgdξ dgdg. (6) From this we derive the ﬁrst important result, that ﬂame speed cdoes not depend on the barrier height (activation energy) of w. The proof of this follows by noting that wcan be additively decomposed into two functions that are symmetric and antisymmetric about g= 0.5;w(g) =ws(g)+wa(g). Barrier height will generally contribute to wsbut not to wa. Because wsis symmetric, it then follows that d ws/dgis antisymmetric about g= 0.5. Recalling that g−0.5 is antisymmetric about g= 0.5, we see that d ξ/dgis symmetric, and consequently /integraldisplay1 0dws dgdξ dgdg= 0, (7) concluding the proof of the above result. Since the speed depends on the antisymmetric part of the chemical potential only, we can re-write it in terms of a template function ˆ wawhere ˆwa(0) = 0,ˆwa(1) = 1 so that w=w0+(w1−w0)ˆwa+w1/2ws (8) wherew0,w1, andw1/2are parameters that can be calibrated for the material. This leads ﬁnally t o the relation c=λL(w1−w0)γ[ˆwa], (9) whereγis a constant determined by evaluating Eq. (6) for the choice of templat e function. From this we see that, in the one-dimensional case, the mobility is proportional to t he regression rate divided by the energy diﬀerence associated with the reaction. As ε→0, we obtain the sharp interface limit, for which the above analysis shows the one-dimensional case to be unchanged. As the di ﬀuse thickness goes to zero, the interface at any point can be approximated as planar (due to the C2restriction on η) within a suﬃciently small neighborhood, allowing the interface motion at all points to be re garded locally as one-dimensional. Therefore, we conclude that Eq. (9) holds as a general statement of the re lation between speed and mobility in the sharp-interface limit. 2.2. Chemical potential specialization We showed in the previous section that the behavior of the model is i nsensitive to the precise choice of chemical potential w(up to scaling), and is in fact impervious to the magnitude of the energy barrier. 7 Therefore we adopt, for convenience, a polynomial approximation to w: w(x,η) =N/summationdisplay i=1ai(x)ηi(x). (10) Proper segregation of the solution to η= 0 and η= 1 requires that w(0) andw(1) be stable local equilibria. To capture the energy of the reaction, the chemical energy at those two p oints (w0,w1, respectively) must be speciﬁed. Finally, an energy barrier must be prescribed as well ( w1/2). Consequently, a suitable polynomial must satisfy the following conditions: w(0) =w0∂w ∂η/vextendsingle/vextendsingle/vextendsingle η=0= 0 w1/2=w(1/2) w(1) =w1∂w ∂η/vextendsingle/vextendsingle/vextendsingle η=1= 0.(11) In order to satisfy the above conditions, wmust be of degree four (having ﬁve coeﬃcients) or greater. A fourth order polynomial is used here for simplicity, and because it has been shown to work well for other models of this type (see, for instance, [53, 54]). Solving the above eq uations for a0, ...,a4produces  a0 a1 a2 a3 a4 = 1 0 0 0 0 0 −11 16 −5 18−32 14 −8 16 −8 /bracketleftBiggw0 w1/2 w1/bracketrightBigg . (12) It is possible to tailor wto a form that is informed by the chemistry, in which the form of wcorresponds to the reaction energy as a function of reaction coordinate η. It is also possible to extend this potential to capture multi-phase (solid/liquid/gas) behavior through the additi on of another minima at a value of η designated to be a liquid phase. In this work, the aim is to show the e ﬃcacy of a simple model in capturing regression rates; these extensions will be considered in future w ork. 2.3. Mobility calculation In composite propellants, multiple constituent species interact chemically to sustain deﬂagration in the thermodynamic regime of interest. Ammonium perchlorate will burn as a monopropellant, whereas fuel binders (HTPB/PBAN) decompose but do not burn on their own [55]. Therefor e a regression rate model must capture the eﬀect of the species concentration and spatial distr ibution on regression rate. We introduce a spatial variable, called the diﬀuse species ﬁeld ( φ(x,t)) to determine pointwise concentration of oxidizer and binder. We adopt the convention that φ= 1 corresponds to pure oxidizer, φ= 0 to pure binder. In principle, φcan be evolved in time to represent morphological changes (such as agglomerat ion of particles) during the burn. In this work, however, the diﬀuse species ﬁel d does not evolve in time. It is deﬁned at all times and for both the burned and unburned regions. The purpose of φis to determine material-speciﬁc properties (like mobility) that are used to evolve η, and it is generally prescribed as an initial condition. Material characteristics can then be represented simply using a m ixture rule. (Note that the use of a mixture rule is not to eﬀect homogenization , but rather to allow a smooth transition between material types 8 binder oxidizercombusion degradationmixing and reaction ﬂame diﬀuse species ﬁeldeﬀective burn speedPhysical Simpli ﬁed Figure 3: Schematic illustrating the concept of the diﬀuse s pecies concentration simpliﬁcation. (Left: Physical) At t he boundary, binder and oxidizer diﬀuse and mix at the species boundary, w hich enables combustion to occur. Far from the boundary, the binder degrades and the oxidizer burns, and the result is a ﬂa me that is close to the solid interface near the species bound ary, but further away from the oxidizer interface and even furthe r from the binder interface. Physically this produces diﬀer ing regression rates in the three regimes of interest. (Right: S impliﬁed) The diﬀuse species ﬁeld φis used to adjust the eﬀective regression rate depending on the species—generally, slow i n the binder, faster for the oxidizer, and fastest in the inte rfacial region. in an explicitly resolved SCP geometry. For a discussion of homogenizati on, and some of the associated challenges, the reader is directed to [56].) We aim to use continuous m aterial variability to capture the eﬀect of material heterogeneity on regression rate. Physically, this is a cons equence of the mass ﬂuxes and diﬀusion causing mixing and facilitating reactions. The diﬀering local mole f ractions that result from the ﬂux from diﬀerent propellant species cause locally diﬀerences in the reac tion rate and ﬂame height above the solid. These diﬀerences greatly inﬂuence the heat ﬂuxes into the solid an d hence the regression rate. Far from the binder/AP boundary, the AP can support a self deﬂagration ﬂame, w hich ensures that there is some heat transfer from the ﬂame back to the solid—driving the surfac e regression. Also far from the interface, but over the binder, any heat fed back to the solid goes to transforming the species and the phase to gaseous fuel species. This is an endothermic process and results in a generally slow binder surface regression. Near the interface, the gaseous fuel and oxidizers have strong gradient s and as a consequence diﬀuse rapidly towards each other. Speciﬁcally, a diﬀusion ﬂame results from the di ﬀusion of fuel and oxidizer together (in stoichiometric proportions) away from the solid surface, because the heat transfer is too fast back to the solid but temperatures are too low for combustion. In other words, the diﬀusion ﬂame is slightly lifted from the solid surface. This lifted diﬀusion ﬂame, and the lean/rich prem ixed ﬂames that sit at the base of the diﬀusion ﬂame, induce very steep temperature gradients normal to th e surface. Since this causes very high heat transfer from the ﬂame to the solid, near the interface, the local regression rate is considerably higher near the solid/solid interface. To capture this phenomenon, it is necessary to modify the property ﬁeld in such a way as to capture local interface-speciﬁc properties; namely, the pointwise burn rate is Leﬀ(x,t) =Lbin(1−φ(x,t))+Loxφ(x,t)+4φ(x,t)(1−φ(x,t))Lint (13) whereLbinandLoxare the mobility of the pure binder and pure oxidizer, respective ly, andLintis an interface mobility correction term. In this way, we capture the eﬀect of inte rface mixing and reaction chemistry 9 (Fig. 3 left) using a reduced order surrogate model for pressure-in formed burn speed (Fig. 3 right), without incurring the cost of a full gas-phase simulation. The diﬀusiveness of the species boundary is determined by a parameter, ζ, which controls the degree of mixing of the two species. The indivi dual mobilities are calculated by combining a pressure power law for cwith Eq. (9), yielding the expression Li=ni(p/ppi)ri λ(w1i−w0i)γ(14) where variables with an isubscript take on diﬀerent values for diﬀerent regimes (binder, o xidizer, and interface). The constants niandriare calibrated power law parameters, λandγ0are model constants, p is the pressure, and w1i,w0iare burned and unburned energies, respectively. Values for all param eters and model constants are enumerated in Table 1. 2.4. Summary We conclude this section by reiterating the essential components of the model, distilling the governing equations, initial conditions (ICs), boundary conditions (BCs), and ot her equations together in a streamlined presentation. We always assume that simulation time starts at t= 0. The solution domain is denoted Ω, and contains two types of boundaries: Dirichlet boundaries ( ∂1Ω) and periodic boundaries. The following is a complete description of the model for the evolution of η: Evolution of η:∂η ∂t=−Leﬀ(x)/bracketleftBigλ ε∂w ∂η−εκ∆η/bracketrightBig ∀x∈Ω, t≥0 (15a) Unburned IC: η(x,0) = 1.0 ∀x∈Ω (15b) Dirichlet BC: η(x,t) =η0(x,t) ∀x∈∂1Ω, t≥0 (15c) Chemical potential: w(η) =/bracketleftbig 1η η2η3η4/bracketrightbig 1 0 0 0 0 0 −11 16 −5 18−32 14 −8 16 −8 /bracketleftBiggw0 w1/2 w1/bracketrightBigg η∈[0,1] (15d) Diﬀuse species ﬁeld: φ(x,t) =φ0(x) ∀x∈Ω, t≥0 (15e) Mobility calculation: Leﬀ(x,t) =Lbin(1−φ(x,t))+Loxφ(x,t) +4φ(x,t)(1−φ(x,t))Lint ∀x∈Ω, t≥0 (15f) Pressure power law: Li=ni(p/p0i)ri λ(w1i−w0i)γi= bin,ox,int.(15g) The primary equation to be solved is the reaction-diﬀusion equation for the evolution of η, Eq. (15a), a nonlinear parabolic partial diﬀerential equation. The entire domain is initialized in an unburned state (Eq. (15b)), and burn is initiated through boundary conditions (Eq. (15c) ) by setting η0= 0 on the xmin face and η0= 1 on the xmaxface. The other faces are periodic. The nonlinearity is introduced by the choice of nonlinear potential w(Eq. (15d)), which is deﬁned through three model parameters w0,w1/2,w1, corresponding to the burned energy, energy barrier, and unburned e nergy, respectively. The diﬀuse species 10 Figure 4: Sample result (from packed spheres case) illustra ting the use of block-structured adaptive mesh reﬁnement at the boundary. The base level grid is very coarse (left). Multipl e BSAMR levels allows for reﬁnement at the ηboundary (left popout). Reﬁnement is based only on η(not onφ), andφis updated only on an as-needed basis to enhance performance (right popout). ﬁeld is based on a prescribed input φ0(x) (Eq. (15e)) and is used to calculate the mobility Leﬀusing a mixture law (Eq. (15f)). The mixture law is determined based on in dividual single-species mobilities as well as an “interface mobility”, all of which follow a pressure-dependent p ower law (Eq. (15g)), with parameters ni,p0i,λ,γgiven as model inputs. The chemical potential weight λ, interface diﬀuse width ε, and interface energy parameter κare all given as model inputs. Because the main governing equation is par abolic, the solution is generally well-behaved for a suﬃciently small timeste p, and can therefore be solved using basic explicit numerical integration techniques. We pause here to note the physical phenomena that are not resolved by th is model. The model does not explicitly simulate the ﬂuid phase, nor does it include any of the s olid-phase mechanics. Instead, the eﬀects of the solid and ﬂuid phase are captured through a reduced order model for the interface mobility, which controls the regression. As will be shown subsequently, it is possi ble to come close to experimental results even with these simpliﬁcations. It is important, however, to emp hasize that a key feature of this model is that since it is an energy-based framework and uses a diﬀuse interface structure, it is easy to incorporate multiphysics eﬀects, including solid and ﬂuid phase phenomena. 3. Computation One of the primary limitations of diﬀuse interface methods is the com putational cost. By deﬁnition, there is a distinct scale separation between the application scale and the diﬀuse interface scale; indeed, if the scales are not properly separated, the model is unlikely to be acc urate. Computationally this presents a challenge, as full resolution of the entire domain is both prohibitive ly expensive and at the same time unnecessary. Adaptive mesh reﬁnement (AMR) is a computational strate gy particularly well suited for diﬀuse interface methods. Speciﬁcally, AMR selectively reﬁne s the diﬀuse interface and coarsens the mesh elsewhereasappropriate. Block-structuredAMR(BSAMR)isapartic ularlyattractivetypeofAMR(Fig.4). The BSAMR strategy is to divide the domain into distinct levels, eac h corresponding to a level of reﬁnement, 11 and then facilitates communication between levels through averaging an d interpolation between timesteps. This data structure is particularly advantageous for performance, as the individual patches can easily be divided across a distributed memory computational platform. BSAMR also allows for temporal subcycling, automatically reducing the timestep and increasing the number of it erations on ﬁner levels. This eliminates common CFL diﬃculties in explicit methods that can arise during me sh reﬁnement. In this work, we use an in-house code “Alamo” to implement and test the diﬀuse interface model [57]. Alamo is built on the AMReX library [58], and implements a variety of phys ical models ranging from ﬂuid mechanics to microstructure evolution. A fully explicit, forward Euler temporal scheme is used to evolve the kinetic equation. The ﬁnite diﬀerence method is used to comp ute all spatial derivatives; speciﬁcally, a second order central diﬀerence stencil is used for the Laplacian of η. . The coarse level timestep is set to 2.5×10−5seconds, and a temporal substep factor of 2 is used between reﬁnement levels, for a ﬁnest-level timestep of approximately 2 ×10−8seconds. Because the governing equation is parabolic, it behaves well when integrated by an explicit solver as long as the timestep is s uﬃciently small enough to prevent instability. Here, the timestep is chosen to be at least a factor of two smaller than that required by the CFL condition. For more discussion on timestep selection in phase ﬁel d type models, the reader is referred to [44, 59]. Regridding is set to occur every 100 timesteps (on the coarse level) and occurs whenever the following condition holds: |∇η||∆x| ≥0.01, (16) where ∆xis the vector of grid spacings (Fig. 4 left popout). These values are based on previous phase-ﬁeld type models implemented in Alamo by some of the authors [43, 45, 57, 60], and are c hecked to ensure that there is no mesh dependency observed in the solution. A full disc ussion of the regridding methods used by Alamo and AMReX is outside the scope of this work. Instead, the reader is d irected to the original AMReX [58] paper, and the associated documentation, and the previously menti oned BSAMR phase ﬁeld papers for a broader discussion of regridding methods in phase ﬁeld applications . Reﬁnement does not currently occur at solid/solid boundaries, unle ss they intersect the burn interface, in ordertomaximizeperformance. Thespeciesﬁeld φisnotevolvedintime, butitcanstillbecomputationally expensive to calculate. The packed spheres case, for instance, can r equire signiﬁcant time to calculate on a highly reﬁned domain, depending on the number of spheres in the sys tem (Eq. (17)). In the BSAMR scheme, φmust be either stored at full resolution, or recalculated on-the-ﬂy as domains are reﬁned and coarsened. In the present implementation, φis calculated on the ﬁnest level only, and is simply averaged/interpol ated elsewhere as the grid evolves. This is acceptable since the evoluti on of the system always takes place on the ﬁnest level only, and is constant elsewhere regardless of the value of φ. When plotting the stored value of φ, however, it may appear that φis diﬀerent in the pre-burn and post-burn regions, due to diﬀeren ces in the way it was interpolated. (Fig. 4 right popout). 12 Table 1: AP, PBAN, HTPB parameters for the solid phase model Parameter name Symbol AP PBAN HTPB AP+ AP+ PBAN HTPB Flame speed exponent (-) n 1.042 0.0 0.0 0.0 0.0 Flame speed coeﬃcient (mm/s) r 1.222 0.1 0.5 10.01 12.5 Base pressure (MPa) pp 1.0 Burned energy (MPa) w0 0.0 Unburned energy (MPa) w1 1.0 Activation energy (MPa) w1/22.0 Diﬀuse width (mm) ε 0.0005 Chemical potential factor (-) λ 0.001 Interface energy (MPa) κ 1.0 Shape constant (-) γ 0.02726 The two-dimensional simulations are generally tractable for a single nod e or high-performance worksta- tion. All of the 2D simulations presented here were run simultaneousl y on the STAMPEDE2 supercomputer, with each simulation using a single node with 48 MPI tasks. The time t o complete varied between runs, but generally averaged approximately 2 hours of wall-clock computation time for each second of simulation time. The three-dimensional simulations, naturally, required substant ially more computational time to complete. The domain size was reduced to 30% of the height of the 2D simulation size to reduce the computational time. The 3D calculations were performed on the INCLINE cluster at the Uni versity of Colorado Colorado Springs. Each calculation was run on 4 nodes with 128 processors per node for a total of 512 MPI tasks. The time to complete also varied between runs, but generally average d approximately 33 hours of wall-clock computation time for each second of simulation time. 4. Results In this section, we apply the computational phase ﬁeld model to study a variety of SCP conﬁgurations, with the aim of establishing its validity and predictive ability. I n this study, we use AP for the oxidizer and consider binders of both PBAN and HTPB. Four diﬀerent cases are consid ered: the coupon of pure AP (pure AP); a thin layer of PBAN sandwiched between AP (AP/PBAN sandwich ); a thin layer of HTPB sandwiched in between AP (AP/HTPB sandwich); and AP sphere particles w ith HTPB binder (AP/HTPB packed spheres) in two and three dimensions. 4.1. Pure ammonium perchlorate burning We begin with a calibration of the model for the pure AP case to establish ed, experimentally measured pressure-dependent AP regression rates. The pressure-depende nt mobility is calculated using Eq. (9), and is reported in Table 1 along with parameters used in this work. Asimulati on domain of 4mm ×0.5mm was used with a base mesh of 16 ×2 cells and 7 levels of BSAMR reﬁnement. Burn rates are calculated and c ompared to legacy experimental data reported in [4]. To achieve consistency w ith experimental measurements, such as those reported in [8, 15, 61], which use a break-wire method to measure b urn speed, amounts to measuring the instantaneous eﬀective position of the interface as the location of t he leading point of the burn front. 13 4mm0.5mm 0.1mm(a) Dimensions of sandwich structure consisting of AP ( gold) and HTPB or PBAN binder layer ( blue). 0.6 0.8 1.0 2.0 4.0 6.0 Pressure (MPa) 2.0 3.0 4.0 6.0 Flame speed (mm/s) Pure AP (Price 1983) Pure AP (Model) AP/PBAN Sandwich (Price 1983) AP/PBAN Sandwich (Model)(b) Comparison of model results for ﬂame speed in pure AP and AP+PBAN sandwich structure Figure 5: AP/PBAN sandwich burn results Reasonable agreement to experimental values are recovered (Reported in Figs. 5b and 7a). This is unsurprising, given the analytic relation between mobility, chemi cal potential, and burn rate; nevertheless, it indicates that the model will produce realistic results for th is type of monopropellant burn. 4.2. AP and PBAN sandwich In realistic SCPs, the behavior induced by the SCP geometry (compl ex interface surface of binder with the packed AP particles) results in complex burn behavior that elud es simplistic analysis. Therefore, to better understand the interactions between oxidizer and binder , sandwich structures provide a simpliﬁed yet nontrivial case study. Sandwich structures have been examined in a variety of experimental reports with varied compositions [62–68], over the span of the past several decades. The ge ometry of sandwich structures is simple: a laminate of binder between two blocks of AP (Fig. 5a). Here, a 4mm×0.5mm geometry is used with a laminate thickness of 0.1mm. As with the pure AP case, a 16 ×2 base grid is used with 7 levels of mesh reﬁnement, and a diﬀuse interface thickness of ζ= 0.015 is used. Pressures ranging from 500kPa to 4MPa are considered, and parameters are adjusted by comparison to the val ues reported by Price et al. [4]. In this simulation the parameters for AP are held constant with the value s that were found for the pure AP combustion and we calibrate the binder parameters only. The resulti ng burn rate calculations appear to be consistent with experimentally measured values, with an average error of approximately 10% for PBAN, compared to less than 5% for pure AP. (Fig. 5b). Observed variation in the v elocity is minimal. The 10% error bars are added to reﬂect the anticipated uncertainty in the comp utational results. An interesting point to note is that, although a pressure power law form had been assume d for both the PBAN and mixture regions, the best exponential coeﬃcient was found to be zero, i.e., a constant value was determined to best match the experiments. A key feature of this model is its ability to capture the complex regr ession morphology in addition to the regression rate. A comparison was made between experimentally obse rved burn front shapes (Fig. 6a) and numerical simulation (Fig. 6b). For the numerical results, traces are plotted at every 10ms; blue traces 14 (a) Experimental observations of burn proﬁles for low press ure (a,d), medium pressure (b,e), and high pressure (e,f); repr oduced from Price et al. [4, Fig 1] (b) Traces of the burn front evolving over time. Blue traces a re plotted at time intervals of ∆ t= 100ms. From left to right the pressures are 0.5, 0.8, 2, and 4MPa. Figure 6: AP/PBAN Sandwich: comparison between experiment ally observed and numerically predicted interface morphol ogy every 100ms. The simulated burn behavior, qualitatively, follows the morphologic al trends observed in experiments. At low pressures, AP burns very slowly because the pressure is t oo low to maintain monopropellant burning. This manifests as a deep grooved burn proﬁle, caused by the diﬀeren tial between the burn speeds in the interface region and the pure AP region. Remarkably, we even observe the slight dimple present in the binder layer, which is the result of an excess of binder beyond that which i s needed to react locally with AP near the interface. At high pressure, where AP is able to eﬀectively se lf-react, an opposite trend is observed: the AP burns faster than the binder, leaving a distinctive protrusion of the binder. 4.3. AP and HTPB sandwich The setup for a sandwich structure composed of AP with an HTPB binder is identical to the AP/PBAN case except for the selection of parameters for the binder and binder i nteraction. Again, the AP parameters are unchanged. Experimental data from Knott et al. [15] is used for calibrat ion of the binder parameters. It should be noted that we include the experimental data for all binde r thickness for comparison from this work, observing that the diﬀerence between them is relatively sm all. Parameters for HTPB are determined by grid search. As with the PBAN case, it has been determined that the pressure law exponent was negligibly small, meaning that a con stant value for the burn rate suﬃced to produce the observed results. Comparison to experimental data i ndicates an approximate error of 5– 20% between the simulated regression rate and the observed rate (Fig. 7a). The 10% error bars have been added here as well. It is of interest that an inﬂection in the burn rate is observed at approximately 8kPa, matching a similar trend in the experiment. Despite the similar mathematical form to the AP/PBAN case, this demonstrates that nonlinear phenomena can be captured with this model. In this case, we attribute this change to the more extreme morphology at the interface. 15 0.2 0.4 0.6 1.0 2.0 4.0 6.0 Pressure (MPa) 1.0 2.0 3.0 4.0 6.0 Flame speed (mm/s) Pure AP (Price 1983) Pure AP (Model) AP/HTPB Sandwich (Knott, 2002) AP/HTPB Sandwich (Model)(a) Comparison of model results for ﬂame speed in pure AP and A P+HTPB sandwich structure. (b) Traces of the burn front evolving over time. Blue traces are plotted at time intervals of ∆ t= 100ms. From left to right the pressures are 0.5, 0.8, 2, and 4MPa. Figure 7: Regression rates in AP/HTPB sandwich structures w ith varying pressure The burn proﬁle is qualitatively similar to that calculated for AP and PB AN (Fig. 7b). In both cases, the interface regresses most quickly at the solid/solid interface, d ue to the enhanced burn chemistry enabled by the mixing of the binder and AP mass ﬂuxes. As the pressure is incr eased, the AP matrix is able to sustain a reaction independently. At high enough pressure, the pur e AP regression rate exceeds that of the solid/solid interface, essentially leaving a ﬂat AP surface. In all cas es, a protrusion formed by the HTPB binder is left behind. At low pressure the binder is quickly c onsumed, but at high pressure it is persistent. We note that such a protrusion would, in reality, be unstable and unlik ely to remain attached as indicated by the model. We expect, however, that the quantity of HTPB not eﬀecti vely burned would likely correspond to the volume left behind in this model. 4.4. Packed 2D AP spheres in HTPB matrix In this section we test the predictive capability of the model by c onsidering the practical case of packed AP spheres (with two diﬀerent mass concentrations) embedded in an HT PB matrix. The initial geometry of the oxidizer particles are modeled as spheres similar to previou s studies (Fig. 8a). This test case is representative of practical SCP, and we emphasize that the model param eters are not adjusted to match the results to experimental values; the same parameters determined i n the previous sections are used, and the only change in this work is the geometry of the binder and matrix. The 2D simulation domain is speciﬁed to be 4mm ×1.6mm and is packed with spheres having a uniform radial distribution of r= 55µm. Mass fractions of 78% and 64% AP are considered, which corresponds to volume fractions of 62% and 45%, respectively. In two dimensions, the corresponding area fractions are 74% and 53%, respectively. A 3D mass fraction of 75% was used, instead of 80% as in th e corresponding experiment, because of the inherent diﬃculty in calculating pack ings with packing factors greater than 62%. Although AP particles are generally irregular, we follow precedent [69] in the use of non-overlapping spheres 16 4mm1.6mm .25mm .25mm(a) 1 2 3 4 5 6 7 Pressure (MPa) 0 2 4 6 8 10Flame speed (mm/s) Kohga (MF 80%) Kohga (MF 75%) Kohga (MF 70%) Kohga (MF 64%) Kohga (MF 63%) Kohga (MF 62%) 2D Model (MF 78%) 2D Model (MF 64%) 3D Model (MF 78%) 3D Model (MF 64%) (b) Figure 8: Results for packed AP spheres in an HTPB matrix. (a) Visualization showing the generated 3D packing and the corresponding 2D slice. The inset 0.25mm section is the regi on used for the 3D simulation. (b) Results for the model with mass fractions of 78% and 64%, compared to experimental data reported in [9]. 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Time (s) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5Burn position (mm) 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Time (s) 0 2 4 6 8 10Burn rate (mm/s) 64%, 1.0MPa 78%, 1.0MPa64%, 2.0MPa 78%, 2.0MPa64%, 3.0MPa 78%, 3.0MPa64%, 4.0MPa 78%, 4.0MPa64%, 5.0MPa 78%, 5.0MPa64%, 6.0MPa 78%, 6.0MPa64%, 7.0MPa 78%, 7.0MPa Figure 9: Burn positions (left) and instantaneous interfac e velocities (right) calculated using central diﬀerence for this calculation. (We note that the eﬀect of non-spherical particle s, as well as multi-modal and lognormal radial distributions, is non-trivial; this will constitute futur e work.) The domain is ﬁlled with a unimodal distribution of spheres corresponding to the aforementioned volume fractions. In order to faithfully represent the 3D sample with a 2D model, a 3D domain is packed with spheres using t he Jodrey and Tory algorithm [70] implemented in OpenMC [71]. A slice of the 3D spheres in the domain is then used to generate the eﬀective unimodal distribution in 2D. To generate the diﬀuse spec ies ﬁeld corresponding for the packed spheres, the formula is used φ(x) = 1−N/productdisplay n=1/braceleftBigg 1 2+1 2erf/bracketleftBig (|x−xn|−rn)/ζ/bracketrightBig/bracerightBigg , (17) whereNis the total number of spheres/discs, {x1,...,xN}are their locations, and {r1,...,r N}their radii. For the 64% and 75% cases, the number of 2D discs was N= 426 and N= 628, respectively. The parameter ζdetermines the diﬀusivity of the solid/solid interface ( ζ→0 recovers the sharp interface limit) and was selected here to be ζ= 10µm. (φ(x) for the packed spheres is plotted in the right pop-out of Fig. 4.) 17 (a) MF = 64% ,p= 1MPa (b) MF = 78% ,p= 1MPa (c) MF = 64% ,p= 4MPa (d) MF = 78% ,p= 4MPa Figure 10: Traces illustrating observed burn front behavio r for varying mass fractions and pressures. Gray traces are p lotted at intervals of 50ms, blue traces at intervals of 250ms. The burn rate is determined as before by ﬁtting the slope of the eﬀec tive interface position using linear regression. Unlike with the pure AP or sandwich cases, a substantial amou nt of variability is observed in the leading point burn location. This is due to the high degree of complexi ty observed in the burn front, which will be discussed shortly. Simulations have been conducted for pr essure from 1–7MPa at 1MPa intervals for both volume fractions, and compared to corresponding experimental ob servations for equivalent SCPs as reported by Kohga [9] (Fig. 8b). The observed error ranges between 5% to 80% error for the 64% mass fraction, and from 1% to 100% for the 70% mass fraction. Errors are calculated based on the RMS diﬀerence between between the instantaneous velocity and the approximate ve locity, and are represented using error bars. The high degree of variance is a result of the high variability in th e instantaneous velocity due to the substantial heterogeneity in the packing. This is clearly observed i n the plots of the burn rate position and the calculated velocity vs. time (Fig. 9). Qualitatively,theresultsdemonstratemarkeddependencebothonp ressureandonconcentration,whichis well-known to exist. This indicates that the approximations made i n the model, despite their simpliﬁcation of complex burn behavior, do not introduce substantial error. It also shows the predictive power of the model, which is especially apparent here, since the parameters we re unchanged from the values obtained in the sandwich calibration tests. (Obviously, a much closer match coul d be attained if the parameters had been calibrated for this case!) We suggest, therefore, that this mode l can be used in its present form to investigate the behavior of other types of SCP conﬁgurations with reasonab le accuracy. As with the sandwich cases, we also investigate the morphological details of the solid phase regression. Unlike with the previous cases, in which the front approaches a steady -state shape, the randomness of the 18 Figure 11: Comparison of burn proﬁles at p= 5MPa for 64% AP (top) and 78% AP (bottom). Intervals are at ∆ t= 0.16s. the sphere packing induces a highly random and richly complex burn f ront. Here, we present trace plots illustrating the progression of the burn front in time (Fig. 10). The tr aces are superimposed on plots of the AP particle locations, and are plotted at 50ms intervals. The trace di agrams illustrate the qualitative diﬀerences between the burn front surface time history. At constant low pressure, we observe that the morphology of the low conce ntration AP (Fig. 10a) is substantially more complex than the high concentration case (Fig. 10b). Th is is unsurprising: low AP concentrations have lower interfacial surface area, and leave large is lands of binder that have no AP with which to react. On the other hand, high AP concentrations result in smal ler and fewer binder islands, meaning that the majority of the binder is available to react with the AP and there is less binder left behind. A similar trend is observed at high pressure (Figs. 10c and 10d). Trends are also compared for both concentrations at 5MPa (Fig. 11); here it is clear that, aside from the obvi ous diﬀerence in ﬂame speed, there is a distinct diﬀerence in morphology as well as the size and shap e of left-behind islands. 19 (a) Burn of 64% MF composition at 5kPa, 180ms intervals (b) Burn of 78% MF composition at 5kPa, 100ms intervals Figure 12: 3D simulation of packed AP spheres in HTPB using th e diﬀuse interface method 4.5. 3D AP spheres packed in HTPB matrix Packed sphere SCP conﬁgurations are inherently three-dimensional. As discussed above, 2D simulations cannot faithfully represent a 3D packed sphere structure, as the 3D e quivalent of any 2D approximation is an array of packed cylinders, rather than spheres. To determine thi s eﬀect of dimensionality, 3D simulations have been performed for the p= 5kPa case, for both 64% and 78% AP mass fractions. A simulation domain of 4mm normal to the burn direction and 0 .25mm×0.25mm in-plane is used, and the same packing conﬁguration is used as that for the corresponding 2D simulations (Fig. 8a). T he two faces on the outside of the original domain do not intersect any of the enclosed spheres; the ot her two include cross-sections of the packed spheres; this allows the simulation to resolve potential edge eﬀects that are only apparent in 3D. Though not the focus of the work, the process of output visualization is n ontrivial and warrants explana- tion. VisIt [72] is used to perform initial processing of the raw BoxLi b output, which entails the calculation of theη≥0.5 isovolume. The geometry is then imported into Blender [73]. Since the raw mesh has vary- ing levels of reﬁnement and does not include full resolution of φit is diﬃcult to visualize the composition. To render the composition over the imported geometry, then, a speci al shader has been developed by the authors using the open shader language (OSL) to eﬃciently encode the sp here packing information without 20 the computational expense of explicitly meshing individual spher es. This allows for surface and volumetric rendering of arbitrary AP/HTPB compositions and geometries using modern rendering techniques. HTPB is rendered using a non-specular semi-transluscent blue shader with subsurface scattering; AP is rendered using an orange, opaque, glossy shader. Each high-resolution 2D render requ ires approximately 10 seconds on a desktop computer to run. Burn rate is calculated similarly as for the 2D case. The η= 0.5 isosurface is calculated for each output ﬁle, and the most extreme point in the isosurface is treated as the eﬀe ctive position. It is possible that instantaneous speeds are overestimated ( vis ´ a vis experimental results) using this method, as a physical trip wireorhigh-speedimagingmaynotdetectextremaoftheburnfrontmor phologythatarecapturedhere. Itis not expected that this would aﬀect the average burn rate, however, onl y the variability in the instantaneous burn velocity. Velocity is calculated using the central diﬀere nce method, and a linear regression is ﬁt to determine average velocity. The variance of the velocity is repre sented as the error. For both the 64% and the 78% case, the calculated 3D burn rate is lower than that predicted by the 2D case, although the error between the 3D prediction and the experiment al measurement is approximately the same (Fig. 8b). We attribute this diﬀerence to two factors. The ﬁrst , and most signiﬁcant, is the diﬀerence in 2D area fraction vs 3D volume fraction for equivalent packings. As mentio ned above, the area fraction is substantially greater than the volume fraction. Therefore, the 3D si mulation equivalent to the 2D cases would have a lower interface content, reducing the speed at which the interface regresses. The second factor is the presence of two non-sphere-intersecting faces, which ha ve a lower burn rate due to the locally lower interface content. Upon inspection of the burn front morphology, howeve r, the average curvature of the burn front appears isotropic (i.e., no obvious “bowing out” of the interface b etween non-sphere-intersecting faces), which indicates that edge eﬀects are negligible. The morphology of the regressing interfaces qualitatively appears to obe y similar trends as that observed in the 2D simulation. The 3D morphology of the unburned ( η≥0.5) region is visualized at 180ms intervals for the 64% case and 100ms for the 78% case. (Diﬀerent intervals are used to all ow for better qualitative comparison of the interface.) For the 64% case, the interface moves more s lowly, as its motion is hampered by large binder peninsulas. The burn progresses via rapid jumps forw ard as the surface encounters AP and solid/solid interface regions, and otherwise moves slowly. The burn leaves behind large islands of binder that would normally be advected away by the ﬂow. The burn progresses quic kly in the 78% case and is rarely slowed substantially by binder obstacles because of the large AP concen tration. This simulation shows that moreislandsofbinderareleftbehind. Althoughtheyaremorenumerous, theislandsaresubstantiallysmaller than those in the 64% case. It should be emphasized that these islands are not “ejected” from the interface but merely left behind. The apparent trail of islands is thus the re sult of the greater speed of the interface. In both cases, since a greater concentration of islands that are located at t he non-sphere-intersecting faces, the eﬀects of edges are signiﬁcant in 3D simulations. 21 5. Conclusion In this work, we present a solid-phase model for calculating the bur n rate and morphology for solid composite propellants using a diﬀuse interface method. The model is based on classic work from phase ﬁeld theory, which allows for the thermodynamics of the system to b e captured faithfully in a geometrically arbitrary setting. We demonstrate the interpretability between t he model parameters and physical values, by showing the relationship between free energy, mobility, interf ace shape, and ﬂame speed. We then develop a scheme for reduced-order modeling of the eﬀect of interface mixi ng on the burn rate, using a heuristic combustion model and diﬀuse species ﬁeld parameter ﬁeld. The model is applied to a selection of experimentally veriﬁable SC Ps. It is ﬁrst calibrated to AP mono- propellant burn by comparison to legacy data. Sandwich structures with both PBAN and HTPB binders are then used to calibrate the binder model parameters for a range of pre ssures. Finally, the model is ap- plied predictively (i.e., without adjusting parameters) to stu dy packed AP spheres in an HTPB matrix at varying mass concentrations and varying pressures. For the sandwich c ases, against which the model was calibrated, the error was observed to be 10% or less in most cases. For the p redictive packed spheres case, the error is more substantial in some cases (up to 100%) but was less than 10% f or the 64% case and when p= 5MPa for both packings. Further development and calibration of the mod el will focus on increasing its predictive capacity for these and other SCP conﬁgurations. Future work will also require validation of the model against other experimental datasets, considering broader an d more realistic conﬁgurations, including the Miller packings [74], that have served as a benchmark f or validation in other computational modeling eﬀorts [75–77]. An important conclusion that can be drawn from this work, regarding solid phase models of this type, is that it is necessary to explicitly account f or interface content. In other words, it is an oversimpliﬁcation to simply rely on individual decomposition/bur n rates; to do so will inevitably fail to capture the increased burn rate observed in SCPs burning below 1M Pa. It is also important to reiterate the connection between pressure, concentration, and islands of left-be hind binder, which may have important implications for unburned binder and combustion eﬃciency. We have proposed a highly reduced-order model intended to capture the burn rate and morphology of SCPs, and we caution that care should be taken in the interpretation and e xtension of the model results. The primary simpliﬁcation is the approximation of the mobility, L, using a pressure and species dependent power law. In reality, Lshould be understood to depend primarily on temperature, which i s determined by the heat ﬂux, which is determined by the reaction chemistry, whi ch is dependent on the mass ﬂux, which in turn is dependent on the mobility. Though the proposed formulation includes several order-reducing sim pliﬁcations, and is apparently able to make some predictions with them, the real power of the phase ﬁeld ap proach lies in its ability to robustly integrate multiphysics phenomena. Speciﬁcally, in the solid phase , the use of continuous ηandφﬁelds allow for the “phasing in” of solid-phase governing equations to model t hermal transport and thermoelastic response. In fact, a robust theory has already been developed for imposi ng boundary conditions over diﬀuse 22 phase ﬁeld boundaries [78], and elastic solvers have already been imple mented with those capabilities [43, 45, 57] to simulate mechanics; all these lay the groundwork for future ex pansion of the model. Of greater consequence, however, is the resolution of the ﬂuid phase—essenti al to capturing behavior such as detailed ﬂame locations, heat ﬂux back to the solid, and erosive burning [79]—whic h is not present in the current model. We refer the reader to current work by some of the authors that p resents a method for robust coupling of solid and ﬂuid phases across a diﬀuse deﬂagrating boundary [51], which will provide the basis for concurrent simulations of the solid and gas phase. Acknowledgments BK and BR were supported by the Oﬃce of Naval Research [grant number #N00014-21- 1-2113]. MQ was supported in this eﬀort by UCCS Department of Mechanical and Aerospace Engineering start-up funds. This work used the INCLINE cluster at the University of Colorado Colorado Sp rings. INCLINE is sup- ported by the National Science Foundation, Grant #2017917. This work also used th e Extreme Science and Engineering Discovery Environment (XSEDE) on the STAMPEDE2 supe rcomputer at the Texas Advanced Computing Center (TACC) through allocation TG-PHY130007; XSEDE is supported by the National Sci- ence Foundation [grant number ACI-1548562]. Finally, the authors gratefully ac knowledge many discussions with Dr. Brian Bojko that substantially strengthened this work.