We combine a thermoacoustic experiment with a thermoa- coustic reduced order model using Bayesian inference to accu- rately learn the parameters of the model, rendering it predictive. The experiment is a vertical Rijke tube containing an electric heater. The heater drives a base ﬂow via natural convection, and thermoacoustic oscillations via velocity-driven heat release ﬂuc- tuations. The decay rates and frequencies of these oscillations are measured every few seconds by acoustically forcing the sys- tem via a loudspeaker placed at the bottom of the tube. More than 320,000 temperature measurements are used to compute state and parameters of the base ﬂow model using the Ensem- ble Kalman Filter. A wave-based network model is then used to describe the acoustics inside the tube. We balance momen- tum and energy at the boundary between two adjacent elements, and model the viscous and thermal dissipation mechanisms in the boundary layer and at the heater and thermocouple loca- tions. Finally, we tune the parameters of two different thermoa- coustic models on an experimental dataset that comprises more than 40,000 experiments. This study shows that, with thorough Bayesian inference, a qualitative model can become quantita- tively accurate, without overﬁtting, as long as it contains the most inﬂuencial physical phenomena. Address all correspondence to this author.NOMENCLATURE Ag Cross-sectional area of gas [m2],Ag=pD2=4 As Cross-sectional area of solid wall [m2],As=pe(D+e) c Speed of sound propagation [m s1] cp;g Speciﬁc heat capacity of gas at constant pressure [m2s2K1] cs Speciﬁc heat capacity of solid [m2s2K1] df Diameter of one of the heater ﬁlaments [m] D Internal diameter of tube [m] Dp Diameter of single prong [m] Dw Diameter of single wire [m] DS Total diameter including two prongs and two wires (if present) [m], DS=D+2(Dp+Dw)orDS=D F Friction force at the solid wall [N] g Acceleration of gravity [m s2] hi Convective heat transfer coefﬁcient from gas to solid [W m2K1] ho Convective heat transfer coefﬁcient from solid to ambient [W m2K1] ki Inviscid pressure loss coefﬁcient at the heater location, ki=Dp=1 2rhU2 h L Length of the tube [m] n Interaction index NuiNusselt number for convective heat transfer from gas to solid, Nu i=hiD=lg NuoNusselt number for convective heat transfer from solid to ambient, Nu o=hoL=la 1 p Gas pressure [Pa] Pr Prandtl number, Pr =rgncp;g=lg ˙q=˙Q Thermal power [W] ˙Qi Thermal power from gas to solid [W], ˙Qi=hipi(TgTs)dx ˙Qo Thermal power from solid to ambient [W], ˙Qo=hopo(TsTa)dx R Reﬂection coefﬁcient R g Gas constant of air [m2s2K1],R g=p=(rT) s Complex frequency of the thermoacoustic system [s1], s=sr+i si si Frequency of oscillation [s1] sr Growth/Decay rate of oscillation [s1] t Time [s] T Temperature [K] u=U Gas velocity [m s1] x Axial coordinate [m] g Heat capacity ratio of air, g=1:4 dbl Boundary layer thickness [m] Dx Length of small acoustic element [m] e Tube thickness [m] QNon-dimensional temperature, Q= (TT1)=T1 l Thermal conductivity [W m1K1] n Momentum diffusivity [m2s1] pi Inner perimeter of solid wall [m], pi=pD po Outer perimeter of solid wall [m], po=p(D+2e) r Density [kg m3] t Time delay [s] twall Shear stress at the tube wall [kg m1s2] () Average gas quantity ˜() Quantity per unit length ()Non-dimensional quantity ()0First-order perturbed quantity ()1 Geometric or gas quantity evaluated at tube inlet section ()2 Geometric or gas quantity evaluated at tube outlet section ()a Quantity at ambient conditions ()d Quantity downstream of heater ()dd Quantity referring to drag device ()e Equivalent quantity ()g Quantity that refers to gas ()h Quantity that refers to heater or evaluated at heater loca- tion ()i Quantity that refers to acoustic element i ()m Quantity that refers to microphone ()p Quantity that refers to a single heater prong ()s Quantity that refers to solid wall ()t Quantity that refers to thermocouple ()th Thermal quantity ()u Quantity upstream of heater()vis Viscous quantity ()w Quantity that refers to a single heater wire 1 INTRODUCTION Thermoacoustic oscillations are nonlinear phenomena that occur whenever heat release rate ﬂuctuations are sufﬁciently in phase with pressure oscillations. They can cause structural vibra- tion of the system or even complete destruction. Predicting and eliminating thermoacoustic oscillations is a signiﬁcant challenge in gas turbine design. Due to the high sensitivity of the growth rate of oscillations to geometry, boundary conditions, and system parameters, models often fail at accurately predicting the insta- bility of a thermoacoustic system. This means that a stable lab- scale engine does not guarantee that its full-scale version will be stable as well. An example of this is the Saturn V engine, which required more than 2000 full-scale tests before a stable conﬁg- uration was achieved [1]. Recently introduced legislation re- garding NO xemission [2], will require gas turbines to operate at lean premixed conditions. These conditions, however, are more conducive to thermoacoustic instability because lean ﬂames are more sensitive to equivalence ratio perturbations [2–4]. This mo- tivates research into quantitative modelling of these phenomena. The goal of this study is to combine statistical methods with a vast amount of experimental data in order to extract quantitatively-accurate reduced order models based on physi- cal principles. This technique is called data assimilation, and was ﬁrst developed for meteorology problems [5]. The prob- lem of having a computational model that depends on parame- ters with unknown value and uncertainty represents a classical inverse problem. Many engineering problems fall into this cat- egory, and therefore many mathematical tools have been devel- oped for this over the past decades. In the ﬁrst part of this study the unknown parameters are inferred using an Ensemble Kalman Filter (EnKF), which is a data assimilation technique, whereas in the second part we infer the model parameters by simple regres- sion. The paper is organized as follows. Sec. 2 provides a de- scription of the physical mechanism that triggers thermoacoustic instabilities in a Rijke tube. Sec. 3 describes the experimental apparatus and the way data are acquired and processed. Sec. 4 presents a model of the base ﬂow that arises by natural convec- tion inside the Rijke tube, and describes how this is combined with the experimental observations. Sec. 5 contains details on the acoustic model, which is built using the information coming from the base ﬂow model. Sec. 6 tests and compares two differ- ent physics-based thermoacoustic models. 2 PHYSICAL DESCRIPTION OF THE RIJKE TUBE In this study we use a Rijke tube, which is a simple ther- moacoustic device that supports thermoacoustic oscillations. It 2 FIGURE 1 : Sketch of the experimental apparatus including 6 mi- crophones, 8 thermocouples to measure the gas temperature at the tube centerline, a thermocouple to measure the ambient tem- perature, a loudspeaker, and an electric heater held in place by a pair of prongs and wires. Note that the loudspeaker acoustically forces the system at a frequency, f, close to the natural frequency of the tube, fn. is a vertical tube, generally open at both ends, which contains a heat source. The heat source sets up a natural convection ﬂow through the tube. On top of the mean ﬂow, small acoustic per- turbations can propagate. For a tube open at both ends, the ﬁrst acoustic mode has a pressure node and a velocity antinode at the two ends of the tube, and a pressure antinode and a velocity node in the middle. Air ﬂows towards the centre of the tube dur- ing the compression phase of an acoustic cycle and away from the centre during the expansion phase. When the heat source is placed in the bottom half of the tube, it experiences a higher- than-average velocity during the compression phase, which in- creases the heat transfer to the ﬂow, and a lower-than-average velocity during the expansion phase, which decreases the heat transfer to the ﬂow. Typically there is also a small time delay be- tween the velocity perturbation and the subsequent heat transfer perturbation. The time delay is such that moments of increased heat transfer occur during moments of increased pressure, and moments of decreased heat transfer occur during moments of de- creased pressure. In the absence of dissipation, this causes the acoustic energy to grow [6]. In the presence of dissipation, typi- cally through visco-thermal mechanisms and acoustic radiation, the acoustic energy grows only if the above driving mechanism FIGURE 2 : Electric heater made up of 0.559 mm diameter Ni-Cr ﬁlaments wound between two parallel ceramic rings. exceeds the dissipation. Ref. [7] found that the optimal posi- tion for self-excited thermoacoustic oscillations is when the heat source is placed at xh=L=0:25. When the heat source is instead placed in the top half of the tube, moments of increased heat transfer coincide with moments of lower pressure, and moments of decreased heat transfer coincide with moments of higher pres- sure. This causes the acoustic energy to decrease [6]. 3 EXPERIMENTAL SETUP 3.1 Apparatus A sketch of the experimental apparatus is shown in Fig. 1. The rig consists of a 1 m long stainless steel vertical tube with an internal diameter of 47.4 mm and a wall thickness of 1.7 mm. An electric heater (Fig. 2) is attached to two rods and held in place at four different positions from the bottom end of the tube: xh=L= [0:55;0:50;0:45;0:40]. The heater is powered by an Elektro-Automatik EA-PSI 5080-20 A DC programmable power supply with maximum power 640 W, controlled through National Instruments LabVIEW. The experimental apparatus, as described so far, is that in [8]. Six G.R.A.S. 40SA probe microphones are used to record the pressure near the inner surface of the tube at different ax- ial locations measured from the bottom of the tube: xm=L= [0:95;0:85;0:75;0:65;0:55;0:45]. Each microphone is equipped with a 20 mm long probe, on top of which a temperature shield is mounted. The raw pressure signal is sampled at 10 kHz, which is much higher than the anticipated frequencies of the thermoa- coustic oscillations, 170 to 190 Hz. All data is acquired through a National Instruments BNC-2110 DAQ device using LabVIEW. A 4WVisaton FRS 8 loudspeaker is ﬁxed at the base of the tube. The loudspeaker is connected to a Stage Line STA-500 Pro 3 0 0.5 1 1.5 2 2.5-150-100-50050100150 0 0.5 100.51FIGURE 3 : Pressure oscillations recorded by the six micro- phones after sinusoidal forcing at 170 Hz. The inset ﬁgure shows the corresponding pressure eigenmode. Power ampliﬁer with maximum power 600 W controlled through National Instruments LabVIEW. Eight type-K thermocouples are installed along the cen- treline of the tube through small holes at positions xt=L= [0:90;0:80;0:70;0:60;0:50;0:40;0:30;0:20]. An additional ther- mocouple is placed near the inlet section to record the ambi- ent temperature. All thermocouples are logged with four TC-08 USB DAQ boxes. In the following discussion, we will assume that the temperature does not vary in the radial direction, so that a 1D model can be used. This assumption could be relaxed, but the modelling would be signiﬁcantly more expensive. 3.2 Data acquisition The experiment is automated through National Instruments LabVIEW. The input power is varied in steps of 10 W, from 10 W to 180 W. Each step lasts 70 minutes. For every heater po- sition, each experiment lasts 21 hours. The nine thermocouples measure the temperature simultaneously every 7 seconds. Ev- ery 7 seconds, a sinusoidal signal at 170 Hz is provided to the loudspeaker. At these heater powers, the system is thermoacous- tically stable, so the acoustic oscillations decay. Their decay rate and frequency depend, however, on the thermoacoustic driving from the heater. The experimental campaign is carried out in the stable regime because the power that can be provided to the heater is limited to 300 W, which is not enough to acquire a large set of experimental data in the unstable regime. Investigation of the unstable regime would require an additional microphone to be placed near the top end of the tube, and a phase-shift ampli- ﬁer for active control of the system. FIGURE 4 : Sketch of a small element of tube of length dxwith gas speed U(x), gas temperature Tg(x), solid temperature Ts(x), ambient temperature Ta, and inner and outer heat ﬂow rates ˙Qi and ˙Qo, respectively. 3.3 Data processing Fig. 3 shows a plot of the pressure signal recorded by the six microphones. To extract the decay rate and frequency, ﬁrstly, the decaying part of the signal is isolated. Secondly, the Fourier transform is applied. Thirdly, the logarithm of the resulting sig- nal is taken. Finally, by means of a lifted cosine weighted func- tion, a weighted least squares regression is performed in order to ﬁt a straight line. The slope of the straight line corresponds to the decay rate, whereas the frequency of the oscillations is given by the frequency where the peak in the Fourier spectrum occurs. This procedure is applied to each microphone, and the six values are then averaged to obtain a more robust estimate, although no signiﬁcant scatter exists between the six measurements. The av- eraged decay rate and frequency represent the real and imaginary parts of the complex frequency s, respectively. 4 BASE FLOW MODEL AND ANALYSIS 4.1 Model We consider a small element of tube of length dx(see Fig. 4). The energy balance in the solid reads (rsAsdx)cs¶Ts ¶t=¶˙Qs ¶xdx˙Qo+˙Qi (1) By assuming that the diffusive heat transfer ˙Qs=ls¶Ts ¶x, and that ˙Qoand ˙Qican be modelled with convective heat transfer coefﬁcients hoandhi, respectively, and by further assuming con- stant thermal conductivity ls, Eq. (1) can be re-written as rsAscs¶Ts ¶t=lsAs¶2Ts ¶x2hopo(TsTa)+hipi(TgTs)(2) 4 By re-arranging and replacing the convective heat transfer coef- ﬁcients with Nusselt numbers, one obtains (3)¶Ts ¶t=ls rscs¶2Ts ¶x2Nuola rscspo AsL(TsTa) +Nuilg rscspi AsD(TgTs) The same procedure can be applied to a control volume con- taining gas moving at local speed U, thus giving (rgAgdx)cp;g¶Tg ¶t+U¶Tg ¶x =¶˙Qg ¶xdx˙Qi+˜˙Qhdx(4) In Eq. (4), the source term˜˙Qhaccounts for the electric heater. The inlet gas quantities ()1are assumed to be ambient ()a. One can now apply the ideal gas law, as well as Fourier’s law, and use mass conservation in its integral form¶(rAU) ¶x=0 , to obtain (5)r1Agcp;gT1 Tg¶Tg ¶t+r1Agcp;gU1¶Tg ¶x =lgAg¶2Tg ¶x2Nuilgpi D(TgTs) +˜˙Qh The gas thermal conductivity is assumed to be equal to that of air at ambient conditions, i.e. lg=la=l1. Radiation is not modelled. We also solve the integral momentum equation, in which we assume that (i) the unsteady term is negligible, i.e. the iner- tia forces are negligible compared with the buoyancy and drag forces, and (ii) the total pressure losses Dpare concentrated at the heater location, and modelled through an inviscid pressure loss coefﬁcient ki, such that Dp=ki(rhU2 h)=2. The integral mo- mentum equation is r2U2 2r1U2 1+kirhU2 h 2=ZL 0(r1rg)gdx (6) By using mass conservation and the ideal gas law, Eq. (6) can be re-arranged to U2 1ki 2A1 Ah2Th T1+A1 A22T2 T11 =ZL 0(r1rg) r1gdx(7) In our case A1=A2=Ah, but in general these areas could be different.4.2 Non-dimensionalization We non-dimensionalize Eqs. (3), (5) and (7) with the ref- erence scalesfL;g;T1g, which naturally give a time scale (L=g)1=2, and a velocity scale (gL)1=2. Temperatures are mea- sured relative to T1and then divided by T1. For example, T becomes Q=TT1 T1. Non-dimensional quantities are denoted by()to distinguish them from the corresponding dimensional quantities. The energy equation for the solid becomes (8)T1 (L=g)1=2¶Q s ¶t=ls rscsT1 L2¶2Q s ¶x2 Nuola rscspo AsLT1Q s+Nuila rscspi AsDT1(Q gQ s) We deﬁne the following non-dimensional parameters h 1=1 L(gL)1=2ls rscs(9) h 2=po As(gL)1=2la rscs(10) h 3=pi As(gL)1=2L Dla rscs(11) and re-write Eq. (8) compactly as ¶Q s ¶t=h 1¶2Q s ¶x2Nuoh 2Q s+Nuih 3(Q gQ s) (12) The energy equation for the gas becomes (13)raAgcp;gT1 (L=g)1=21 (Qg+1)¶Q g ¶t +raAgcp;gT1 (L=g)1=2U 1¶Q g ¶x=laAgT1 L2¶2Q g ¶x2 Nuilapi DT1(Q gQ s) +˜˙Qh We deﬁne the following non-dimensional parameters h 4=1 L(gL)1=2la racp;g(14) 5 h 5=pi Ag(gL)1=2L Dla racp;g(15) ˙Q h=(L=g)1=2 T1˜˙Qh raAgcp;g(16) and re-write Eq. (13) compactly as 1 (Qg+1)¶Q g ¶t+U 1¶Q g ¶x=h 4¶2Q g ¶x2Nuih 5(Q gQ s) + ˙Q h (17) The integral momentum equation becomes (18)Z1 0Q g Qg+1dx=U 12ki 2A1 Ah2 (Q h+1) +A1 A22 (Q 2+1)1 4.3 Numerical implementation A 1D Finite-Difference scheme with 101 points is used for spatial discretization of the base ﬂow model. The scheme is fourth-order accurate and is centred everywhere except at the boundaries, where it is biased. Zero temperature gradients are imposed at the inlet and outlet boundaries. An explicit fourth- order accurate Runge-Kutta scheme is used for time integration. The source term in Eq. (17), ˙Q h, is modelled with a Gaussian distribution centred at the heater location, with variance set arbi- trarily to 0 :0002 m2, normalized such that its integral equals the input power. The second-order accurate midpoint rule is used for the integral in Eq. (18). 4.4 Bayesian inference The 1D unsteady model is marched forward in time to re- produce the evolution of the experiment. This, however, requires knowledge of the three parameters Nu i, Nu o, and ki. We wish to exploit the information from both the model and the experi- ments in order to obtain the best possible estimate of the gas and solid temperature proﬁles. We must consider that our base ﬂow model is nonlinear in the temperature Q g. The Ensemble Kalman Filter is a suitable tool to address these problems [5]. This repre- sents the statistically optimal method when working with Gaus- sian distributions of state and parameters, in the sense that it pro- vides the minimum variance estimate or, equivalently, the max- imum likelihood estimate of f(yjd)[5]. The idea, which relieson Bayes’ theorem, is to ﬁnd the optimal estimate for state and parameters, given the marginal probability of the model often called the prior ,f(y) , which is the probability of the model with no knowledge of the data, and the probability of the data given the model often called the likelihood function, f(djy) . The product of these two probability density functions, once properly normalized, gives the probability of the model given the data often called the posterior ,f(yjd) , which represents the statistically best estimate that one can compute, given the model and the experimental data with the associated uncertainties. This is summarized by f(yjd) =f(y)f(djy) f(d), which represents Bayes’ theorem. In the present case, the state is a vector containing the solid and gas temperatures at every grid point, and the parameters are Nu i, Nu o, and ki. The EnKF randomly samples from the prior distributions and simulates every member of this ensemble until experimen- tal observations become available. The prior distributions of the state and the parameters are deﬁned by their mean and standard deviations. We set the initial state to be the ambient state with zero variance. We set the initial parameters [Nui;Nuo;ki]Tarbi- trarily to have mean [7;60;9:5]Tand variance [0:2;0:3;0:1]T. We assume that the experimental temperature measurements have uncertainty4 K. We use an ensemble made up of 30 members. Once the temperature measurements become available (every 7 seconds), state and parameters are updated using Bayes’ theo- rem. Because the model is nonlinear, the initial Gaussian shape of the prior distribution deforms into a non-Gaussian distribu- tion. Before the update step, the mean and covariance of this new distribution are computed from the ensemble members and a new Gaussian distribution is assumed, with the same mean and covariance. The update state is performed and the optimal state and parameters are estimated with their uncertainties, based on the likelihood function. The optimal distributions then become the new prior distributions to sample from, and the process re- peats. With this process, one is able to (i) infer the evolution of the model parameters with their uncertainties, and (ii) improve the estimate of the state, which becomes more robust by learning from the experimental measurements. 4.5 Analysis Figs. 5a to 5c show the evolution of the model parame- ters from low powers (10 W) to high powers (180 W), with an uncertainty of 2s. The results for xh=L= [0:40;0:50] have higher uncertainties and less regular trends than those for xh=L= [0:45;0:55]especially as far as Nu iis concerned. This is because for xh=L= [0:40;0:50]the thermocouple measure- ments at the heater location are affected by radiation and con- duction, and are therefore discarded, being these mechanisms not included in our model. This is the most informative loca- tion (see Fig. 5e). Fig. 5d shows the inlet bulk velocity, which increases to a plateau at high powers. Despite the imperfection 6 (a) Inner Nusselt number. (b) Outer Nusselt number. (c) Pressure loss coefﬁcient. (d) Inlet bulk velocity. (e) Gas temperature at t=21 h. (f) Outlet gas temperature. FIGURE 5 : (a), (b), (c): Base ﬂow model parameters with uncertainty 2sas functions of time ([h], bottom axis) or heater power ([W], top axis). (d): Bulk velocity of the gas at the inlet of the tube. (e) Gas temperature distribution inside the tube at t=21 h ( ˙Qh=180 W) with corresponding experimental measurements. (f) Outlet gas temperature evolution. 7 FIGURE 6 : Sketch of a network model with ﬁve acoustic ele- ments separated by four jump conditions. The travelling waves in each acoustic element are denoted by fiandgi.R1andR2are the inlet and outlet reﬂection coefﬁcients. in the model, the gas temperature (Fig. 5e) seems to capture well the experiments even for the largest discrepancy case ( ˙Qh=180 W). Finally, Fig. 5f shows the evolution of the outlet gas temper- ature, which is particularly inﬂuential for the phase of the outlet reﬂection coefﬁcient. 5 ACOUSTIC MODEL AND ANALYSIS 5.1 Model To model the acoustics, we split the domain into 50 acoustic elements, each with constant density, pressure, and temperature, and connect these with a network model [9]. Forces (shear stress at the wall) and heat transfer (to and from the wall) are integrated across adjacent acoustic elements and then modelled by jump conditions between the elements. The choice of using 50 acous- tic elements is simply motivated by the fact that a smaller num- ber does not signiﬁcantly reduce the computational cost. Fig. 6 shows a sketch of a network model with ﬁve acoustic elements separated by four drag devices , denoted by ()dd, representing four jump conditions. We assume that entropy waves have no inﬂuence and, knowing that the mean ﬂow is much slower than the speed of sound, assume zero mean ﬂow. The linearized mo- mentum and energy equations in each acoustic element are r¶u0 ¶t+¶p0 ¶x=0 (19) ¶p0 ¶t+gp¶u0 ¶x=0 (20) The solution of the above equations is obtained by performing a Riemann decomposition, which results in a f-wave travelling downstream and a g-wave travelling upstream. The pressure andthe velocity ﬂuctuations in the generic acoustic element iare p0 i(x;t) =fi tx ci +gi t+x ci (21) u0 i(x;t) =1 rici fi tx ci gi t+x ci (22) We take Laplace transforms of the travelling waves, hence fi tx c =Fi(s)estesx ciandgi t+x c =Gi(s)este+sx ci. The jump conditions between the elements model the visco-thermal phenomena in the system. The ﬁrst is the heat and momentum transfer that occurs in the boundary layer due to velocity and temperature gradients at the tube walls and around the heater prongs and wires. The second is the heat and momentum trans- fer where the heater and thermocouples are inserted. The third is the thermoacoustic mechanism, which we deal with separately in Sec. 6. The radiation of acoustic energy through the tube ends is accounted for by a reﬂection coefﬁcient, experimentally calcu- lated and assumed to be constant for all the experiments. The visco-thermal drag in the boundary layer is distributed inside each acoustic element. From [10], the ﬂuctuating wall shear stress is twall=(rn=dbl)u0, where dbl=2p(2n=si)1=2. Consider now an element of tube with length Dxand total perime- terpDS, with DS=D+2(Dp+Dw). This diameter allows us to account for the presence of the two prongs and wires to which the heater is attached. The total ﬂuctuating force on the ﬂuid element is DF=DxpDStwall=Dx DSr 23=2(sin)1=2u0(23) We now deﬁne an equivalent area, and corresponding equivalent diameter, to account for a reduced ﬂow area in the acoustic ele- ments that contain heater prongs and wires: Ae=AApAw=p 4(D22D2 p2D2 w) =p 4D2 e (24) (If an acoustic element does not contain heater prongs and wires, then we set DS=De=DandAe=A.) Integrating the momen- tum equation across adjacent acoustic elements gives the pres- sure jump as a known linear function of the acoustic velocity u0 (25)p0 i+1(t)p0 i(t) =DF(t) Ae =DxDS pD2er(2nsi)1=2u0 i(t) kvis;blu0 i(t) 8 Similarly, the heat transfer from the wall into the gas is D˙q=DxpDSl dblT0(26) where dblis assumed to be the same as for the viscous boundary layer because the Prandtl number for air is close to 1. The gas is assumed isentropic and ideal, soT0 T=p0 pg1 gandp=rR gT. By combining the previous relationships, we obtain the heat transfer as a known linear function of the acoustic pressure p0: D˙q=Dx DS(sin)1=2 23=2Prp0(27) This can be written as a jump condition for the acoustic velocity. Integrating the energy equation across adjacent acoustic elements gives (28)u0 i+1(t)u0 i(t) =g1 g1 p4 pD2eD˙q(t) =Dxg1 g1 p(sin)1=2 23=2Prp0 i(t) kth;blp0 i(t) By combining Eqs. (21), (22), (25) and (28), the momentum and energy jump conditions at xddcan be written in terms of FandG as (29)Fi+1esxdd ci+1+Gi+1e+sxdd ci+1Fi 1kvis;bl rici esxdd ci Gi 1+kvis;bl rici e+sxdd ci=0 (30)Fi+1esxdd ci+1 ri+1ci+1Gi+1e+sxdd ci+1 ri+1ci+1 +Fi kth;bl1 rici esxdd ci +Gi kth;bl+1 rici e+sxdd ci=0 For the ﬁrst and last acoustic elements F1esxu c1R1G1e+sxu c1=0 (31)R2FNesxd cNGNe+sxd cN=0 (32) Eqs. (29) to (32) can be written in matrix form as A(s)W=0 (33) where sis the complex frequency of the system, whose real and imaginary parts represent growth rate and frequency of the ther- moacoustic oscillations. Wis the vector of complex amplitudes [Fi;Gi]T, with i=1; ::: ; N, and Nthe total number of acoustic el- ements. Further modelling is described in the following section, where some results are presented. 5.2 Analysis The decay rate and frequency of oscillations in the empty tube, averaged over 1100 identical experiments, are s=7:61 rads1+i168:10 Hz. We assume that the two ends are acoustically identical, i.e. R1=R2=R, and solve Eq. (33) for R, given s. This gives a reﬂection co- efﬁcient for the empty tube of R=0:9758 +i0:1003. Once the reﬂection coefﬁcient is known, the heater is placed inside the tube at different positions, without switching it on. The following 11 positions are investigated: xh=L= [0:55;0:50;0:45;0:40;0:35;0:30;0:25;0:20;0:15;0:10;0:01]. For each case, 100 identical experiments are performed. The heater is modelled as a visco-thermal drag element with a time delay, following the theoretical analysis of [11] and knowledge of the qualitative form of the feedback sensitivities in [10]. Hence the jump conditions Eqs. (25) and (28) at the heater location are p0 i+1(t)p0 i(t) =kvis;hu0 i(ttvis;h) (34) u0 i+1(t)u0 i(t) =kth;hp0 i(ttth;h) (35) The four parameters kvis;h,kth;h,tvis;h, and tth;hare obtained by ﬁnding the best ﬁt to Eqs. (33) to (35), given 100 measurements ofsat each of the 11 heater positions. Fig. 7 shows the predictions of the network model compared to the experimental measurements for all heater positions. The calibrated parameters are reported in the caption. The time of the viscous drag, tvis;h, is negative, while the thermal time delay, tth;h, is positive, in agreement with [11]. The pressure loss coef- ﬁcient due to viscous drag at the heater location, kvis;h, is several orders of magnitude larger than the thermal loss coefﬁcient, kth;h. Nevertheless, setting kth;hto zero results in relatively large errors in Fig. 7, showing that the thermal transfer to the heater cannot be neglected, even when the heater is switched off. 9 0 0.2 0.4 0.6-12-10-8 Exp Model 0 0.2 0.4 0.6164166168170FIGURE 7 : Decay rate (top) and frequency (bottom) of acoustic oscillations when the heater is switched off and placed at differ- entxh=L. The calibrated heater parameters are: kvis;h=27:3 kg m2s1,kth;h=1:53105kg1m2s,tvis;h=1:36103s, andtth;h=8:98104s. We tried building several other models, which all failed to predict the correct trend. For example, we tested different com- binations of models in which the pressure jump is proportional top0or the velocity jump is proportional to u0. This is in per- fect agreement with the predictions of the feedback sensitivities in [10] (see Fig 7a in [10]). The 8 thermocouples are then placed in the tube and the decay rate and frequency of acoustic oscillations are measured with a further 100 experiments at each of the 11 heater posi- tions. There is a small but noticeable shift in the decay rate (top of Figs. 7 and 8). We introduce 8 new jump conditions into the model, one at each thermocouple location. We then ﬁx the ratio kvis;t=kth;tto be the same as that of the heater, kvis;h=kth;h, and set the time delays to zero, for simplicity, so that the visco-thermal drag from the thermocouples is quantiﬁed by a single parame- ter. As before, we perform regression on Eq. (33) to calculate the optimal value of this parameter. This leads to accurate results, shown in Fig. 8, with the resulting best-ﬁt parameters reported in the ﬁgure caption. 6 THERMOACOUSTIC MODEL COMPARISON When the heater is switched on, in addition to accounting for the visco-thermal drag losses that occur in the boundary layer and across the heater and thermocouples (Sec. 5.2), one has to ac- count for the thermoacoustic mechanism. By integrating Eq. (20) 0 0.2 0.4 0.6-12-10-8 Exp Model 0 0.2 0.4 0.6164166168170FIGURE 8 : Decay rate (top) and frequency (bottom) of acoustic oscillations when the heater is switched off and placed at differ- entxh=L, and thermocouples present inside the tube. The cali- brated thermocouple parameters are: kvis;t=0:257 kg m2s1, kth;t=1:44107kg1m2s. across the heater, the jump condition is u0 i+1(t)u0 i(t) =kth;hp0(ttth;h)+g1 g1 p4 pD2˙q0 h(t)(36) A common model for ˙ q0 h(t)is the ntmodel [12], according to which the normalized heat release rate ﬂuctuations are pro- portional to the normalized velocity ﬂuctuations upstream of the heating element with a time delay: ˙q0 h(t) ˙Qh=nu0 h(tt) U1(37) Assuming that the parameters nandtin Eq. (37) are constant, one can do regression using the remaining 43200 experimental decay rates and frequencies at different heater powers in order to infer their values. One can then compare the performance of the model to the experiments. Results are shown in Fig. 9 and the best-ﬁt parameters are reported in the ﬁgure caption. The agreement between the lines is fairly good and the main source of discrepancy occurs when predicting the frequency. The error becomes larger as the heater power increases. This is mainly because of the constant reﬂection coefﬁcient assumed for this study. In reality, the phase of the reﬂection coefﬁcient changes with the temperature, therefore, as the heater power increases, it 10 0 5 10 15 20 25-20-15-10-5 0 5 10 15 20 25170180190FIGURE 9 : Decay rate (top) and frequency (bottom) of the oscil- lations for the case of heater switched on and placed at different xh=L. The heater power is increased by 10 W every 70 minutes (bars at the bottom). Red lines: model. Black lines: experiments. Optimal thermoacoustic parameters: n=0:271,t=0:517 ms. is less and less valid to assume constant R2(Fig. 5f). On the other hand, the decay rate trend is correctly captured with fairly good accuracy for all the conﬁgurations. The second model we propose is similar to Eq. (37) but tis modelled as t=kt0:2df U1(38) where dfis the diameter of one of the many heater ﬁlaments (see Fig. 2), which is 0.559 mm. The choice of using this model takes inspiration from [11], who states that ”For frequencies small compared with w0=20U1=df, the departure of the heat-transfer ﬂuctuations from their quasi-steady form consists essentially of a time lag of the order of 0:2df=U1”. By doing regression over nandkt, one ﬁnds that the best-ﬁt ndoes not change, whereas kt=2:04, which is close to the value of 1 implied by Lighthill’s statement. With this value of kt, one obtains a time delay tin the range [0:451:20]ms. The results are plotted in Fig. 10. The frequencies predicted by the network model are similar to those in Fig. 9. The modelled decay rates, however, seem to be negatively affected by the dependence on velocity, especially at low powers, where the velocity shows signiﬁcant transients (see Fig. 5d). Overall, the ﬁrst model, where nandtare assumed to be constant, seems to perform slightly better. 0 5 10 15 20 25-20-15-10-5 0 5 10 15 20 25170180190FIGURE 10 : Decay rate (top) and frequency (bottom) of the os- cillations for the case of heater switched on and placed at dif- ferent xh=L. The heater power is increased by 10 W every 70 minutes (bars at the bottom). Red lines: model. Black lines: experiments. Optimal thermoacoustic parameters: n=0:271, kt=2:04. 7 CONCLUSIONS In this paper we use Bayesian inference on a thermoacoustic system and assimilate 324,000 temperature measurements and 46,500 decay rates/frequencies to infer the parameters of a re- duced order model. The experiment consists of an electrically- heated Rijke tube in the thermoacoustically stable regime. Every 7 seconds we acoustically pulse the system with a loudspeaker at the bottom of the tube, and measure the subsequent decay rate and frequency with six probe microphones, sampling at 10 kHz. We measure the temperature of the air inside the tube at 8 axial locations, and use these data with their uncertainty to improve the predictions of an unsteady 1D conjugate heat transfer model for the base ﬂow. For this, an Ensemble Kalman Filter with 30 members is used to provide statistically optimal estimates of state and parameters. This allows us to compute the inlet bulk veloc- ity and the gas temperature at every grid point. The base ﬂow model informs an acoustic network model, made up of 40 drag devices to account for visco-thermal losses in the boundary layer, and 9 drag devices to account for visco-thermal losses across the heater and the 8 thermocouples. We perform 1100 identical ex- periments with the empty tube in order to compute the reﬂec- tion coefﬁcients at the top and bottom ends, assuming them to be acoustically identical. We then place the heater inside the tube at 11 different locations without switching it on, and perform 100 11 identical experiments for every heater position. This allows us to infer the model parameters of the visco-thermal losses across the heater. We ﬁnd that the viscous loss is characterized by a negative time delay, whereas the thermal loss has a positive time delay, in agreement with what found by [11]. We repeat the same pro- cedure after placing the thermocouples inside the tube and work out through regression the inﬂuence of the thermocouples on the acoustic ﬂow and thermal ﬁelds. Finally, we perform experi- ments with the heater switched on and investigate four different heater positions. We vary the heater power by 10 W every 70 minutes. We use the parameters learned in the cold tube cases to learn the parameters of two thermoacoustic models we propose: a classical ntmodel with constant parameters, and a modiﬁed version from [11], in which t=kt(0:2df)=U1. This paper shows that an Ensemble Kalman Filter of many conjugate heat transfer simulations can assimilate experimental data from a rig, even if it never reaches steady state. Without this assimilation, the sound speed in a corresponding thermoacous- tic network model is too inaccurate for the model to be quan- titatively accurate. With this assimilation, combined with re- gression to infer another 7 parameters from 46,500 datapoints, a quantitatively accurate thermoacoustic model can be created. Although not shown here, several other models were tried, all of which failed to model the experimental data qualitatively, let alone quantitatively. Both heat release models perform well. Al- though we ﬁnd that the ﬁrst model performs slightly better than the second, we conclude that, as long as a model is qualitatively correct, it can be made quantitatively accurate over the range ex- amined. Further improvement of this method will be to infer the re- ﬂection coefﬁcients with the multi-microphone method. This will allow us to account for the inﬂuence of the high tempera- ture at the outlet end, which affects the phase of the reﬂection coefﬁcient and, consequently, the model predictions of the os- cillation frequency. We will also infer the parameters of more complex thermoacoustic models derived via system identiﬁca- tion based on high-ﬁdelity numerical simulations [13] in order to assess systematic model error (epistemic error). Another im- provement will come from replacing regression with a more ro- bust and quantitatively informative technique in order to provide uncertainty estimates also for the acoustic and thermoacoustic parameters. For this part of the work, a more natural choice than EnKF would be to use Markov Chain Monte Carlo, a Bayesian inference method that would allow us to relax the Gaussian dis- tribution assumption albeit at a higher computational cost. In addition to this, in our future work we will train a purely data- driven machine learning algorithm and compare its performance to our physics-based statistical learning technique. Finally, in- vestigating different conﬁgurations and acquiring new data will allow us to test whether a physics-based model can extrapolate in a reliable manner.ACKNOWLEDGMENT This project has received funding from the European Union’s Horizon 2020 research and innovation programme un- der Grant Agreement No. 766264.