Concise Mechanistic Model of Phase Change Material Solidification Kinetics

1 Concise Mechanistic Model of Phase Change Material Solidification Kinetics Louis Desgrosseilliers1, Nerea Uranga2, Daniel Carbonell3, Ignacio Gurruchaga1 1 SPF Institut für Solartechnik, Ostschweizer Fachhochschule, Rapperswil, Switzerland 2 Fundación Tekniker, Eibar, Spain 3 DCarbo Energy Consulting, Sant Llorenç Savall, Spain E-mail: [email protected] Abstract A concise power law reaction kinetic expression of reversible crystallization phase change was presented for use in global modelling of crystallization and nucleation kinetics of solid-liquid phase change materials (PCMs). This was developed as a closed form, wholly mechanistic expression of phase change kinetics essential to predict coupled heat-mass latent heat evolution in PCMs. This constituted a significant departure from the semi-empirical and phenomenological formulations that have so far dominated PCM sciences. Notably, the presented formulation expresses mass supersaturation as the fundamental driving force for phase change as opposed to relying on supercooling degree or temperature rate. In contrast to practice in the industrial crystallization process industry, it was postulated that crystal size population balances could be neglected in PCMs used for thermal storage due to averaging effects of crystal size population balances since the rates of latent heat evolution and heat transport are exclusively prioritized rather than yield of a desired crystal size. A favourable formulation of the global power law crystallization kinetic model was obtained through evaluation against the 0D transient reference case derived from the cooling curve of manually seeded crystallization of 11 K supercooled 38.1 %mass NaCHO2 aqueous solution, yielding NaCHO2•3H2O(s), under development as a cold storage PCM. This formulation was based on sequential reactions of nucleation and crystallization, the latter linearly dependent on the mass fraction of nuclei. Values of apparent activation energies of nucleation and crystal growth and their respective pre-exponential rate constants were obtained by least-squares fitting of cooling curve data. Keywords: mixture PCMs, solid-liquid phase change, phase change materials, power law kinetic model, nucleation rate, crystallization rate, non-isothermal kinetics This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 2 1. Introduction A versatile and simple, though mechanistic, macroscopic treatment of solid-liquid phase change kinetics would offer distinct advantages in the dynamic analysis and simulation of phase change materials (PCMs) and the performance of thermal storage systems. Such a formulation would readily find application in simulation and analysis of coupled heat-mass latent heat processes for thermal storage, such as non-equilibrium calorimetric analysis of PCMs, study of heating/cooling rate dependent phase change hysteresis, solidification from supercooled liquids, and dynamic thermal storage system performance simulation. In recent years, semi-mechanistic [1] and phenomenological [2] approaches have been introduced to model supercooling and/or phase change hysteresis in static PCM (i.e., no convection transport), each one constituting a departure from classical moving boundary problems of idealized phase change governed only by heat balances. Barz et al. [2] have presented a convenient phenomenological framework to predict phase fraction from path dependent scaling relationships of the apparent, hysteresis enthalpy profile of a PCM. While shown computationally efficient and apparently more representative of experimental data of incomplete phase change cycling than were classical models of idealized phase change thermodynamics, phenomenological approaches cannot fundamentally elucidate knowledge of the inherent phase change processes. Also, as this was modelled using enthalpy-temperature profiles obtained through dynamic calorimetric experiments, such as those used to report enthalpy-temperature in PCM supplier datasheets, the data could be valid only for heating/cooling rate intensities (W.kg-1) near those of the calorimetric measurement conditions. In contrast, Günther et al. [1] developed a semi-mechanistic approach that, while remaining classically mechanistic in its treatment of thermodynamic properties, heat transfer, and propagation of secondary nucleation, relied on explicit empirical functions obtained from isothermal experiments of single crystal growth extrapolated on the bulk to govern the rate of phase change. In consequence, their method produced results of which the phase change solidification kinetics were fixed to prescribed profiles. Also, their model required explicit definition of a nucleation temperature derived from small scale experiments, thus not representing nucleation as a path dependent dynamic process. While their model could be used to reproduce observed phase change behaviour under identical conditions to those of their source data for informing the model equations, it would be however incapable of predicting behaviour of the same PCM under substantially different experimental conditions, namely in the prediction of nucleation from supercooling/supersaturation of the liquid phase. Furthermore, both methods [2] and [1] have been developed for congruently melting PCMs only (i.e., solid phase and liquid phase have identical compositions), thus neglecting phase change of incongruent PCM mixtures, namely aqueous solutions of salt hydrates. Both [1,2] have opted to simplify modelling by adopting expressions using exclusively temperature (e.g., supercooling degree) and temperature rate as the sole driving forces for phase change. However, especially concerning solidification processes, temperature has only an indirect effect on the driving force for phase change, whereas it is in fact a mass supersaturation that is the responsible driving force, long since recognized in the fields of industrial crystallization and metallurgy [3–6]. With temperature, however, one can compute the mass supersaturation by determining the equilibrium phase fraction that is material-system specific and can be obtained from either a heat balance (congruent PCMs) or liquidus curve data from mixture phase diagrams (incongruent PCMs). Many more analogies of phase change can be drawn between different material systems on the basis of mass supersaturation than are possible using supercooling degree alone. While both classical nucleation theory (CNT) [4,7] and law of mass action models of bulk crystallization [3,4,6] employ mass supersaturation (whether directly or indirectly) as the underpinning driving force for nucleation and crystal growth, their formulations are too esoteric and cumbersome for practical use in PCM science and macroscopic simulation of PCM thermal storage systems undergoing solid/liquid phase change. CNT is too often concerned with the dynamics of either only primary nucleation or of single crystal growth and generally requires considerable material property data generally unavailable in the domain of PCMs. Lane introduced fundamentals of CNT to PCM science in his seminal work [8], but this was not further integrated or adapted for macroscopic analysis of phase change processes for thermal storage. Law of mass action models (i.e., power law reaction kinetics) of bulk crystallization, on the other hand, use much more concise formulations with reliance on fewer physical parameters, but are focused on computation of population balance modelling of crystallites rather than just global mass balances of phase fraction yields. This is of course motivated by the economic value of controlling yield of a desired crystal size and purity in industrial production. The JMAK equation (Johnson-Mehl-Avrami-Kolmogorov), otherwise called Avrami equation [4,5], though famously used for its extreme simplicity of exact solution accounting for impingement of growing crystals, represents only isothermal crystallization processes, thus unsuitable for cyclical and/or non-isothermal phase change process modelling and simulation. In the field of PCM science and engineering for thermal storage, the prediction and control of crystal size and population balances are not generally of any great concern, whereas dynamic prediction of bulk phase change yield and its effect on the overall thermal energy balance and heat transfer constitute the primary objectives. Furthermore, in PCM science and engineering, kinetic expressions of phase change must also be able to dynamically achieve various states of 2-phase equilibrium. This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 3 In this paper, a generalized power law formulation of macroscopic, reversible solidification phase change kinetics is presented and further demonstrated for the case of non-isothermal cooling, heterogeneously triggered nucleation in a salt hydrate PCM-water solution and subsequent equilibrium-stage solidification. This specific focus was motivated by the authors’ study of T-history cooling curves for salt hydrate-water solutions manually nucleated from the supercooled state that were performed under the auspices of the EU Horizon funded BEST-Storage project (www.best-storage.eu). The material system in question has been investigated for use in a cold storage Phase Change Slurry thermal storage system. The power law phase change kinetic framework has been examined with respect to single step and two sequential step formulations and explored for the effects of apparent reaction order and apparent activation energy. These were studied using as reference case T-history cooling curve data measured for NaCHO2-water solution, manually seeded in the supercooled state. 2. Method This work focused exclusively on power law kinetic expressions with Arrhenius temperature dependence and supersaturation driving force of the form: d 𝑤 s d 𝑡 = sign ( 𝜎 ) 𝐾 exp ( ― 𝐸 a 𝑅(𝑇 + 273.15) ) | 𝜎 𝑛 | , (1) 𝜎 = 𝑤 l ― 𝑤 ∗ l 𝑤 ∗ l , and (2) 𝑤 ∗ l = 𝑓 ( 𝑇 ) , (3) where w is the mass fraction, σ the supersaturation quotient [3,4], K the Arrhenius pre-exponential rate constant, Ea the apparent activation energy of the net process, R the ideal gas constant, T the temperature in °C, and n the reaction rate order. Subscripts s and l and superscript * indicate the bulk solid phase, solute in the liquid phase, and chemical equilibrium property, respectively. The apparent reaction order n is defined 1 ≤ n ≤ 3 for crystal growth and secondary nucleation – where crystal growth ≤ secondary nucleation [3] – and as high as n ≈ 20 for primary nucleation [4]. Lastly, f (T) represents a polynomial fit to the liquidus curve of the solidifying species in the mixture’s equilibrium phase diagram. The compiled list of symbols can be found at the end of this document. This formulation was founded on four key axioms: 1. Arrhenius temperature dependence of the reaction rate constant, 2. Liquid phase supersaturation driving force solely responsible for solidification processes, 3. Single expression of reversible solidification process, therefore able to achieve the equilibrium state, and 4. Solidification rate is generally independent of the size of crystals. Although appearing deceptively simple, the second point was essential to avoid concluding in error that solidification processes exhibit anti-Arrhenius behaviour (i.e., negative activation energy), as was remarked by polymer crystallization scientists in the early 1950s and later acknowledged to be false [9]. Therefore, the observation in most solidification from liquids where solidification rates increase at temperatures below their thermodynamic equilibrium transition temperatures is owed to the increase in mass supersaturation at lower temperatures (except for mixtures with liquidus curves converging to the solidus at low temperatures). It is important to note also two importance reasons for selecting the liquid phase supersaturation quotient as promoted by Myerson [3] as driving force rather than the quotient of the solid phase conversion deficit used by others (e.g., Foubert et al. and Mazzanti et al. [4]): i) the liquid phase supersaturation represents the excess reagents for crystallization, more apt for law of mass action kinetics than the reaction product deficit, and ii) nonzero value of the denominator wl* at T*, especially important for modelling of single component, congruent PCM solidification. The third point was exploited to simplify the reaction rate expression to avoid separately expressing the reverse reaction rate in order to achieve the equilibrium state, wl*. In contrast to Foubert et al.’s formulation specifying separate forward and reverse rate expressions [4], the current formulation can accomplish reaching equilibrium in a concise formulation due to always preserving the sign of σ in Eq.(1) no matter the reaction rate order ‘n’. For this study, however, the phase change kinetics power law model was used only to model monotonic solidification, so it remained to be seen how the model would perform when representing net dissolution of the solute, especially at temperatures exceeding the saturated liquid temperature of the average composition of the mixture PCM. While item 4 contradicts practice in industrial crystallization, it is nonetheless an appropriate simplification for PCM thermal storage systems. In the development of PCM thermal storage systems, there is no effort taken to ensure control of crystal size populations, thus leaving the control of the heating and cooling rates free to respond to only the operational requirements (external heat transport boundary conditions) for heat absorption and heat release. In the case of solid-on-coil PCM heat exchange, a polycrystalline bulk layer is formed. Therefore, these typical conditions make it plausible to assume that there would exist at any time a broad and random distribution of crystal sizes and nuclei. Therefore, it was assumed that an averaging This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 4 effect dominates the solidification rate of PCMs, in that the birth of new crystals, the growth of mature crystals, impingement, and the fracture of large crystals all contribute simultaneously to bulk solidification with neither dominating. Of course, this would not always hold true (e.g., very low, constant cooling rates generating few nuclei and larger individual crystals), but was expected to be nonetheless an expedient and useful concept within the thermal energy storage field. Assumptions on the most likely anticipated reaction mechanisms (single step vs. two sequential steps) and reaction orders relevant to PCMs have been examined in the following sections for the formulation of the complete system of equations representing 0D transient solidification of solute from a mixture PCM. 2.1. Single Reaction Step The formulation remained as it was in Eqs.(1)-(3). This formulation represented a system where the nucleation rate is sufficiently strong that bulk solidification does not exhibit a distinct induction period. This could be potentially suitable to represent PCM materials exhibiting minimal supercooling and those doped with nucleation additives. In this mode, the reaction order was constrained to 1 ≤ n ≤ 3, representing fast kinetics at low values of the supersaturation quotient, σ. Fast kinetics were thought to be necessary to reproduce equilibrium solidification behaviour using only a single reaction step model, whereas multiple reaction step models could combine contributions of terms each representing fast or slow kinetics. Independent of this mode succeeding to represent the reference case, this mode would serve mainly to study separately the properties of reaction orders 1-3 and prioritize which to be assigned in the study of the two sequential reactions formulation as well as define anticipated ranges for their reaction parameters Ea and K. 2.2. Two Sequential Reactions The defining feature of this mode was the sequential nature of nucleation and crystal growth rates where the crystal growth rate was directly proportional to the mass fraction of nuclei, nucleation serving as the overall rate limiting step. This formulation was anticipated to produce an initial induction period in which, after seeding, the nucleation rate would dominate, and up to a critical value of nuclei mass fraction, crystal growth would dominate the overall solidification kinetics. The governing equations of this mode were: 𝑤 s = 𝑤 n + 𝑤 c , (4) d 𝑤 n d 𝑡 = sign ( 𝜎 ) 𝐾 n exp ( ― 𝐸 a ,n 𝑅 ( 𝑇 + 273.15) ) | 𝜎 𝑛 | , and (5) d 𝑤 c d 𝑡 = 𝑤 n sign ( 𝜎 ) 𝐾 c exp ( ― 𝐸 a ,c 𝑅 ( 𝑇 + 273.15) ) | 𝜎 𝑚 | , (6) where subscripts n and c denote nucleation and crystal growth, and n, and m represent the respective powers of nucleation reaction order, and crystal growth reaction order. As has been already indicated in Section 2, the apparent reaction order for secondary nucleation is defined within the 1 ≤ n ≤ 3 limits and the one for the crystal growth is constrained to m ≤ n. For unseeded (primary) nucleation, one could substitute Eq.(5) with a formulation of sequential/parallel reactions of primary (denoted by subscript “1”) and secondar nucleation (“2”), d 𝑤 n d 𝑡 = sign ( 𝜎 ) [ 𝐾 n,1 exp ( ― 𝐸 a ,n,1 𝑅 ( 𝑇 + 273.15) ) | 𝜎 𝑛 1 | + 𝑤 s 𝐾 n,2 exp ( ― 𝐸 a ,n,2 𝑅 ( 𝑇 + 273.15) ) | 𝜎 𝑛 2 | ] , (7) with the secondary nucleation’s dependence on a prior population of crystallites, ws and n1 > n2 [3]. This, however, was not pursued for demonstration in this study as the reference case for validation was one of manually seeded nucleation (secondary nucleation). 2.3. Materials and Reference Experiment The reference case for demonstration and partial validation of the features of the reaction mechanisms presented in Sections 2.1. and 2.2. was that of a custom prepared T-history cooling curve experiment performed using solutions of NaCHO2 in water. The T-history apparatus and one sample tube assembly are shown in Figure 1. The apparatus consisted of a precision heating/chilling bath (±0.01 K (manufacturer statement), Julabo DYNEO DD1000F) filled with 10 %mass ethanol in water operated at 60 % pump output (27 L/min at 100 %). Three sample tubes were suspended from a custom lid fixture into the bath: each consisted of 10 mm OD/8 mm ID/178 mm long borosilicate glass sample tube surrounded by polyethylene foam insulation of 33 mm OD/14 mm ID and soft polyvinyl chloride sleeve added over the glass tubes to improve the fitting. Insulating the sample tubes followed refinements of the T-history technique at ZAE Bayern [10] and helped ensure that each sample could This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 5 accurately rely on only one single temperature probe despite the strongly convective environment of the chilled bath. Lumped mass approximation of the individual samples was supported by Bi = 0.04 (Biot number) determined experimentally using water as the sample and assumptions of both rotationally-symmetric (dT/dθ ≈ 0) and axially-independent heat transport (dT/dz ≈ 0) due to very low cooling rate, high viscosity (~20 mPa.s ), and poor aspect ratio of sample volume suppressing internal natural convection. Each sample tube was also equipped with a top-mounted 3-port connector (b in Figure 1), allowing for nucleation agent to be added manually without disturbing the T-history experiment. Nucleation agent was introduced via the horizontal port to manually trigger solidification from the supercooled liquid state (closed over with tape when not in use), and a 2 mm OD T-type stainless steel sheathed thermocouple probe was inserted and positioned approximately along the tube’s centreline via the top port of the connector. The heights of the sample tube thermocouples were also adjusted for the sensor tip to lie at approx. the mid-depth of the sample PCM volume (≈ 5 cm from the bottom of the tube). The thermocouples were calibrated to ±0.17 °C uncertainty (95 % confidence) between 0 - 17 °C and logged at 0.1 Hz using a Keysight DAQ970A multichannel temperature transmitter/logger. The sample tubes were each filled with approx. 5 mL of solution prepared from 37.6±0.1 %mass NaCHO2 anh. (oven dried, 97 %mass purity reagent grade from Sigma, CAS 141-53-7, Table 1) and balance of reverse osmosis water, deemed a pseudobinary mixture of 38.1±0.1 %mass (95 % confidence) NaCHO2 in water when omitting the mass of impurities. For such a low fraction of impurities, it was believed that there should be negligible effect of higher order mixtures on the phase equilibrium relationship of water-NaCHO2, especially considering that negligible differential heat of solution has been reported for the analogous system of water-NaC2H3O2 [11]. This analogy was justifiable by the comparable aqueous dissociation constants of each compound’s conjugate acid [12]. The tube assemblies with loaded samples were initially conditioned at room temperature (~16 °C) to achieve thermal equilibrium, then inserted simultaneously into the chilled bath at 0±0.01 °C. Nucleation agent NaC2H3O2 anh. (99 % purity ACS grade from Scharlau, CAS 127-09-3, Table 1) was manually added to the sample tubes once each sample reached < 1 °C (supercooled > 10 K), at which point nucleation was induced (recorded as 5 290 s for sample 1). The samples then reached their respective equilibrium conditions (deemed > 5 460 s for sample 1), followed by a slow cooling process approaching equilibrium-stage solidification for the remainder of the experiment. See sample thermograms in Figure 2 for more detail. The nucleation and solid-liquid equilibrium domains where solidification was detected by the effect of latent heat release on the sample temperature served as the reference case for solutions to the two kinetic model formulations. Table 1: Material provenance Species CAS Supplier Purity NaCHO2 141-53-7 Sigma 97 %mass NaC2H3O2 127-09-3 Scharlau 99 %mass Reverse osmosis water - Onsite - As seen in the water-NaCHO2 equilibrium phase diagram in Figure 3 (produced from digitized data [13] and IUPAC-NIST tabulated data [14]), nucleation of 38.1 %mass NaCHO2 in water below 10.6 °C yields the trihydrate solid compound. Serving as the measured data for least-squares fitting the models proposed in Sections 2.1. and 2.2., the corresponding equilibrium mass fractions of solid trihydrate, ws*, assumed to have been yielded during the solid-liquid equilibrium cooling phase in Figure 2, were inferred using the polynomial fit to the trihydrate compound liquidus (solid-liquid) curve in Figure 3 and the Lever Rule (see Eq.(9)): 𝑤 ∗ l = 1.101 × 10 ― 4 𝑇 2 + 6.024 × 10 ― 3 𝑇 + 3.045 × 10 ― 1 ± 0.001 , and (8) 𝑤 ∗ s = 𝑤 ave ― 𝑤 ∗ l 𝑤 TH ― 𝑤 ∗ l , (9) where T is in °C, wave corresponds to the binary mixture’s average composition (38.1±0.1 %mass NaCHO2), wTH is the stoichiometric mass fraction of NaCHO2 in the trihydrate compound (55.720±0.0005 %mass), and for which Eq.(8) obtained R2 = 1.000 though using only the 4 data points featured by IUPAC-NIST Solubilities Database. However, in liquidus curve determination, 4 points should suffice to provide the two termini (peritectic and eutectic) and two points in between to establish curvature and steepness of slope for a quadratic fit. Needed to calculate the mass supersaturation quotient, σ, in Eq.(2) are: 𝑤 w = 1 ― 𝑤 ave 𝑤 TH = 0.316 ± 0.002 , and (10) 𝑤 l = 𝑤 TH ( 1 ― 𝑤 w 1 ― 𝑤 s ) , (11) This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 6 where ww corresponds to the mixture’s mass fraction of water in excess of the stoichiometric ratio needed to form the trihydrate compound. The ordinary differential equation (ODE) proposed in Section 2.1. was solved as an initial value problem for ws and the system of ODEs in Section 2.2. solved for [wn, wc], both using the ode15s stiff solver in Matlab R2023a to reduce computation effort in the equilibrium solidification regime where solution stiffness was anticipated. The relative and absolute tolerances for ode15s were specified as 10-7 and 10-10, respectively. Initial values for solid fractions were each null and the time interval of the solution was specified as the beginning of the nucleation domain to the end of the solid-liquid equilibrium domain shown in Figure 2 ([5 290, 10 540] s for sample 1). The remaining kinetic parameters, Kn, Kc, Ea,n, and Ea,c were determined by least-squares fitting of the respective ODE solutions to experimentally derived ws*,exp (Eq.(9) calculated for the sample 1 cooling curve at t ≥ 5 460 s in Figure 2). Note that ws* in Eq.(9) was valid only for equilibrium solidification and could not be used to estimate solid fraction in the nucleation domain (5 280 ≤ t ≤ 5 460 s in Figure 2) marked by fast solidification kinetics. Parameter convergence was achieved using fmincon constrained optimization solver in Matlab, set to the default interior-point algorithm. In each case, integer values of n and m 1 were specified and varied to study the properties of the models. The general form of the least-squares optimization objective function and constraints was: min ( ∑ ( 𝑤 ∗ s,exp ― 𝑤 s ( [ 𝐾 i ] , [ 𝐸 a , i ] ,𝑇,𝜎 ) ) 2 ) s.t. { 𝐾 i > 0 𝐸 a,i > 0 , (12) where [] indicates array sets of each K and Ea, index i = [n,c] (n = nucleation and c = crystal growth) and ws is the ODE timeseries solution of total solid fraction. The ODE timeseries solution, ws, was retimed in Matlab (retime function for timeseries data type) using piecewise cubic spline interpolation to obtain corresponding values at each time interval of ws*,exp (0.1 Hz). Also, in solving ws, linear interpolation was used to estimate T from the sample 1 cooling curve data in Figure 2 when the timestep used by ode15s solver did not coincide with the cooling curve’s 0.1 Hz data capture. 3. Results To simplify comparison of the two 0D kinetic solidification phase change model formualtions, the experimentally derived equilibrium solid yield, ws*,exp, was calculated using Eq.(9) for only the sample 1 equilibrium domain cooling curve in Figure 2 (t ≤ 5 460 s) as the other two temperature curves were closely overlapped with that of sample 1 and so would serve no further purpose. Therefore, all of the model derived solutions, ws, were generated for a single data set for direct comparison to oneanother. The initial time used for the onset of nucleation and crystallization in the computation of the kinetic models was inferred from the sample 1 cooling curve in Figure 2-b at which time the temperature profile first showed signs of monotonic selfheating from latent heat release due to nucleation/crystallization, deemed t = 5 290 s. The end state of the nucleation domain, deemed t = 5 460 s, was ascertained by the observation of the cessation of this monotonic self-heating phase. The comparison of converged solutions obtained using the two different phase change kinetic model formulations presented in this study was only the first phase of comparison, as the reference case could not yet serve as a precise validation case of the models, in particular since the dynamic evolution of the solid fraction during the nucleation phase could not be determined directly from the experimental results. This will instead be reported in the following study of the detailed, coupled heat-mass partial differential equation (PDE) solutions to all three T-history cooling curves in Figure 2, as only the temperature curve could serve to accurately reveal which of the phase change kinetics solutions most faithfully represented the dynamic, nonisothermal evolution of both seeded nucleation and crystallization. The results of this study would therefore generate the shortlist of kinetic parameters judged most promising for this final stage of validation. 3.1. Single Reaction Step Table 2 shows the list of converged solutions for n = 1 found from the search for local minima around initial guesses for Ea in the range [10, 104] J.mol-1 and for K in the range [10-3, 101] kg.kg-1.s-1. Initial guesses supplied to fmincon were fixed at each value of Ea = [101, 102, 103, 104] J.mol-1 with the initial guess order of magnitude varied for K until one or more local minima were found or deemed infeasible. In cases where local minima were found, initial guesses Ea and K were repeated close to those minima to further probe the neighbouring solution space. 1 At first, n and m were also obtained by least-squares optimization, but these solutions did not differ significantly from the initial guesses of n and m as there were many local minima possible for K n, K c, E a,n, and E a,c for each n and m . This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 7 Note that the converged solutions for n = 1 in Table 2 were obtained only after adding to Eq.(12) the additional minimization constraint max ( 𝑤 s ( 𝑡 ≤ 5 460 𝑠 ) ) ≤ 𝑤 ∗ s,exp ( 5 460 𝑠 ) (13) to suppress runaway crystallization in the nucleation domain as shown in Figure 4. While fundamentally unrealistic, this oscillating runaway crystallization behaviour served at least to demonstrate the inherent feature of the proposed general form of the power law kinetic model formulation to encompass both the forward and reverse reactions, producing therefore a return force to the initally observed runaway reaction and later reaching the corresponding solid-liquid equilibrium states. Table 2: Converged kinetic parameters and coefficient of determination of fit for single reaction model with n = 1. # Ea (J.mol-1) K (kg.kg-1.s-1) R2 1a 9.98×103 4.50×10-1 0.966 1b 9.98×103 4.57×10-1 0.968 1c 1.01×103 1.00×10-2 0.981 1d 1.01×103 1.56×10-2 0.994 1e 9.15×102 1.38×10-2 0.992 1f 9.15×102 1.35×10-2 0.992 1g 2.74×102 8.10×10-3 0.982 1h 2.74×102 1.10×10-2 0.993 1i 2.00×102 9.20×10-3 0.989 1j 6.07×101 9.30×10-3 0.992 Whether from the experimental body of knowledge on crystallization from liquids, CNT or JMAK equation isothermal crystallization kinetics, there was no accepted basis for this type of seeded crystallization oscillating runaway kinetic behaviour. While there were no means available to directly measure the solid mass fraction during the nucleation step, the general consensus in the literature on both isothermal and continuously cooled non-isothermal crystallization kinetics proposes that nucleation solidification invariably follows a sigmoidal progression, therefore strictly monotonic [4] and therefore nonoscillating. The corresponding lists of convergence results for n = [2,3] are not shown as their local minima followed similar patterns to solutions obtained for n = 1, thus only the solid fraction profiles of their deemed global minima solutions are shown herein. The converged solutions for n = 1 in Table 2 formed solution families grouped strongly around initial guesses of Ea in the range [60, 104] J.mol-1, meaning that for each Ea initial guess, fmincon iterated K in preference to Ea to achieve convergence, indicating a much stronger sensitivity to K than for Ea. For initial guesses of Ea ≥ 105 J.mol-1, fmincon failed to iterate both K or Ea any further from the initial guesses. However, there remained still to identify the globally optimum and physically relevant converged solution for n = 1. A brief survey of activation energies of nucleation and crystallization processes of ceramics and metals/metalloids suggested that scaled from 800 – 1 100 K to ≈ 300 K (peritectic temperature of NaCHO2•3H2O = 18 °C [14]) using the law of equivalent states (equal entropies of transition), activation energies for both nucleation and crystal growth would be estimated in the range ≤ 105 J.mol-1 [15–17]. Therefore, along with the observed convergence difficulties for the initial guess 105 J.mol-1, these reinforced a reasonable upper bound limit of Ea < 105 J.mol-1 for both nucleation and crystal growth. The survey of the n = 1 solution space shown in Table 2 revealed that the apparent global optimum was situated in vicinity to solution #1d. It became clear that those with Ea < 103 J.mol-1 (solutions #1g-1j) were unrealistic since at such low values, the exponential term would become mostly invariable in the sample 1 cooling curve range of 0-18 °C, as demonstrated by the converged values of K having remained nominally equal to 10-2 kg.kg-1.s-1. A comparison of the Arrhenius exponential term at Ea = [102, 103, 104] J.mol-1 for 0-18 °C made it clear that both Ea = [102, 104] J.mol-1 represented extreme cases of invariability (either near 0 or near 1). For the case Ea = 103 J.mol-1, the exponential term was evaluated in the range 0.6-0.7 for 0-18 °C and its slope was 10x greater than for either Ea = [102, 104] J.mol-1, both nominally 10-4 K-1. Therefore, solutions with Ea ≈ 103 J.mol-1 represented cases with most marked temperature dependence of the power law rate expression. It became clear also from values of R2 in Table 2 that those solutions with Ea ≈ 104 J.mol-1 (solutions #1a-1b) were unfavourable to faithfully represent the apparent equilibrium crystallization behaviour of NaCHO2•3H2O. Therefore, solution #1d became the deemed globally optimal solution in Table 2 as it had also superior R2 to solutions #1c and #1e-1f. Figure 5 shows the apparent globally converged solutions for the cases of n = [1, 2, 3], where the solutions for n > 1 were obtained in the same manner as was done for n = 1 (solution 1d). In the case of n = 2, Ea = 9.81×102 J.mol-1 and K = 0.134 kg.kg-1.s-1; for n = 3, Ea = 9.98×103 J.mol-1 and K = 39.9 kg.kg-1.s-1. All solutions were numerically fitted to the This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 8 experimentally derived nominal ws*,exp timeseries data, reserving a sensitivity study to the upper and lower bounds of ws*,exp for only the most promising kinetic model solutions. As seen in Figure 5 (b), the runaway crystallization kinetics in the nucleation domain were successfully abated by the optimization condition in Eq.(13) for all instances of n = [1, 2, 3]. However, that is where the similarities ceased between solutions of different reaction order. The first order apparent optimal solution (n = 1) achieved the most suitable kinetics of equilibrium solidification, ws*,exp (i.e., at low values of σ), albeit at the expense of poor convergence to ws*,exp at the transition from the nucleation domain to the equilibrium domain in proximity to t = 5 460 s. Both the second and third order optimal solutions, while poorly converged overall to ws*,exp, exhibited more rapid rise in modelled ws in the nucleation domain than did the n = 1 solution. These were the results of faster kinetics than for n = 1 when σ > 0.1, due in both cases to orders of magnitude larger values of the Arrhenius exponential term compensating for the reductions in the power law σ term. What the two higher order solutions also revealed was their inherent property to sharply dampen kinetics below a threshold value of σ: for the converged n = 2 solution this was σ < 0.024, and σ < 0.044 for n = 3. While certainly not improving model fit to experimentally derived ws*,exp in Figure 5, this feature could prove valuable in isolating nucleation kinetics when only above a σ threshold value in two sequential reactions modelling, allowing therefore crystal growth kinetics to dominate below this threshold. This type of behaviour would be consistent with observed norms in industrial crystallization from solution [3]. Therefore, while the first order (n = 1) single reaction model showcased fast kinetics at low values of σ to suitably represent equilibrium crystallization kinetics, it could not fundamentally alone capture both nucleation-dominated kinetics and a transition to crystal growth dominance. The second and third order models failed to exhibit equilibrium crystallization kinetics but showed promise as a tool to regulate the transition from nucleation dominance to crystal growth at a threshold value of σ in a mixed order, two reaction model. Furthermore, as seen in Figure 5-b, neither first nor second and third orders single reaction model solutions could exhibit sigmoidal profiles nor nuclei induction periods that are typical of crystallization processes from liquids [4,5]. Failure of the single reaction model to solely capture all aspects of nucleation and crystal growth were anticipated but served nonetheless as a simple case study of the individual merits of reaction orders n = [1, 2, 3] to inform selection of reaction orders for the two sequential reactions models. As such, a sensitivity study for the n = 1 single reaction model kinetic parameters was not performed. For the remaining work on two sequential reactions modelling, it became clear that m = 1 would be most suitable to represent crystal growth in the equilibrium solidification domain via Eq.(6) where low values of σ would be encountered. Both n = 2 or 3 were to be explored for the secondary nucleation kinetics in Eq.(5) due to their abilities to effectively dampen nucleation kinetics upon reaching a threshold value of σ. 3.2. Two Sequential Reactions In the search for local optima obtained by least-squares fitting of the two sequential reactions kinetic phase change model, Matlab’s patternsearch constrained optimizer was used to further search the solution space near local minima identified using fmincon in the same manner as was done for finding single step reaction model local minima. This was done to probe the nearby space more thoroughly using a non-gradient, grid-based optimization algorithm. The solutions obtained are shown in Table 3. Solutions were obtained with initial guesses of types Ea,n = Ea,c, Ea,n > Ea,c, and Ea,n < Ea,c to explore the implications of higher or lower nucleation activations energies for both Ea,n and Ea,c in the range [102, 103, 104] J.mol-1. Values where Ea,n > Ea,c were thought valuable to explore as they could potentially prove consistent with observations for some PCMs where semi-ordered molecular clusters were found to persist in the bulk liquid phase at T > Tm (melting temperature) and were believed to serve as seed nuclei in the cooling phase [18,19]. Higher values of net Ea,n for reversible nuclei formation would theoretically result in nuclei persisting in the bulk liquid phase at T > Tm but also that new nuclei are more difficult to form at T < Tm, resulting in deeper supercooling/supersaturation when the initial nuclei population is sufficiently low due to liquid phase superheating pretreatment to suppress bulk crystallization. The latter feature is represented in the formulation and results presented herein of the two sequential reactions model, while the former remains to be demonstrated in similar kinetic modelling of melting/dissolution processes. In Table 3, two types of solutions were obtained (sample profiles shown in Figure 6 and Figure 7): 1) very close fitting (most with R2 = 0.995-0.997) with kinetic behaviour in two stages marked by fast nucleation terminating in a plateau before the end of the nucleation domain, then near-equilibrium crystallization for the remainder (noted as “two-stage” in Table 3); and 2) close fitting (R2 = 0.930-0.990) smooth kinetic profiles without an intermediate plateaux (noted as “smooth”). Of these two solution types, the smooth profiles were the most compelling. While both solution types exhibited initial nucleation induction phases consistent with crystallization from solution (owed to wn in Eq.(6)), only the smooth profiles exhibited fully sigmoidal This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed 9 progression by also gradually approaching the equilibrium solidification domain which was believed to better conform to the reference temperature T-history temperature profile in the nucleation domain (Figure 2) than two-stage profiles due to the slower and more gradual crystallization rate (hence rate of latent heat release) of the smooth profiles. However, such questions remained inconclusive in this study for lack of suitable data for direct validation of the model solutions in the nucleation domain. Comparisons of the cooling curves in Figure 2 to the resulting simulated heat-mass coupled model temperature profiles obtained using the two-stage and smooth converged kinetic solutions in Table 3 would be further examined in a related study on the Thistory cooling curves of seeded, supersaturated 38.1 %mass NaCHO2-water solutions. Table 3: Two sequential reaction model least-squares solutions for n = [2,3] and m = 1. Grey shaded rows indicated the apparent global optima. Low scoring of solutions due to very high wn was based on the conventional notion that nuclei, while numbering many, due to their small size represent generally very little mass fraction overall (assumed < 0.01). # n m Ea,n (J.mol-1) Ea,c (J.mol-1) Kn (kg.kg-1.s-1) Kc (kg.kg-1.s-1) max(wn) R2 Comment 2a 6.14×10-4 6.32×102 8.00×10-4 1.30×101 0.001 0.9967 Two-stage, very low Ea,n and low Ea,c 2b 5.79×101 1.01×103 8.00×10-4 1.33×101 0.0017 0.9948 Smooth, very low Ea,n 2c 7.04×101 1.00×103 5.00×10-4 1.22×101 0.0013 0.9787 Smooth, very low Ea,n 2d 1.01×102 1.00×103 1.00×10-4 1.01×102 0.0003 0.9966 Two-stage, very low Ea,n 2e 7.13×102 1.00×104 7.90×10-3 1.00×102 0.012 0.9964 Two-stage, very low Ea,n 2f 9.15×102 9.96×103 4.51×10-2 7.24 0.07 0.9695 Smooth, very high ws,n 2g 9.96×102 1.01×103 7.61×10-2 1.14×10-1 0.096 0.9805 Two-stage, very high ws,n 2h 1.00×103 9.99×102 5.20×10-3 1.02 0.011 0.9195 Smooth, lower fit than 2m 2i 1.01×103 1.00×104 8.10×10-3 9.92×101 0.01 0.9955 Two-stage 2j 3.78×103 1.00×104 4.30×10-3 7.47×102 0.0018 0.9969 Two-stage 2k 9.98×103 1.12×103 5.13×10-2 2.00×101 0.0015 0.9967 Two-stage 2l 9.99×103 9.99×103 3.60×10-1 1.16×102 0.01 0.9957 Two-stage 2m 1.00×104 1.34×103 1.95×10-2 1.90×101 0.0008 0.9402 Smooth, best fit n = 2 2n 2 1 9.83×104 9.64×104 1.00×1017 1.00×1017 0.07 0.9327 Smooth, Kn & Kc very high, very high ws,n 3a 1.58×10-2 6.60×102 3.00×10-2 1.60 0.013 0.9955 Two-stage, very low Ea,n 3b 1.37×102 1.03×103 2.68×10-2 1.59 0.01 0.99 Smooth, low Ea,n 3c 2.39×102 9.98×103 3.06×10-2 1.00×102 0.01 0.9958 Two-stage 3d 1.00×103 1.00×104 1.51×10-2 9.92×101 0.005 0.9337 Smooth, lower fit than 3g 3e 9.98×103 1.33×103 9.80×10-2 3.28×101 0.0006 0.9861 Smooth, similar to 3g 3f 9.98×103 8.98×102 9.00×10-2 4.90×101 0.0005 0.9962 Two-stage 3g 9.99×103 9.63×102 9.33×10-2 3.28×101 0.0006 0.9897 Smooth, best fit n = 3 3h 3 1 1.00×105 9.98×102 5.00×1015 1.00×102 0.0002 0.9933 Smooth, Kn too high Until such time, four solutions shown in Table 3 emerged as the most promising (see rows shaded in grey): solutions #2j and #3f for two-stage profiles; #2m and #3g for smooth profiles. Among the two smooth solutions, the third order solution (#3g) exhibited a notable improvement in fit to the reference data (R2), while the two-stage solution #2j showed the best fit overall (followed very closely by #3f). Figure 6 shows the characteristic differences between two-stage and smooth solutions, using the two best fit solutions of each kind – #2j and #3g – as representative examples and comparing them to the single reaction step best fit solution, #1d. Figure 7 compares the two smooth solutions of differing reaction orders, #2m and #3g where the differences in R2 in the same solution family were largest. Other than their wn profiles differing by max(wn), the differences between both the crystal growth fraction, wc, and the total solid fractions, ws, of #2j and #3f were insignificant. In Figure 6-a can be seen the two-stage profile plateau for #2j and the smoothness of the #3g profile in comparison. Figure 6-b further shows that although #3g has ultimately an attenuated overall crystallization rate in comparison to #2j, both profiles followed similarly rapid growths after the initial stages of nucleation induction, the latter clearly absent in the #1d profile. The high initial values of σ were surely responsible for this despite their differences in total nucleation yield, wn. However, the crystal growth rate for #2j remained higher at lower values of σ nearing the end of the nucleation domain, while that of #3g became slowed in comparison. In the equilibrium-solidification domain in Figure 6-a, both profiles are seen to converge for t > 5 800 s (including #1d). Therefore, their individual distinctions were strongest in the initial phase only, principally in the nucleation domain, therefore where conditions were conducive to fast reaction rates. This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=5378887 Preprint not peer reviewed