S Schropff
email: solene.schropff@univ-amu.fr

F Petitpas
email: fabien.petitpas@univ-amu.fr

E Daniel
email: eric.daniel@univ-amu.fr

Compressible single-phase flows in heated and cooled ducts

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

This paper proposes analytical or semi-analytical solutions of one-dimensional, single phase, steady compressible flows in a non-adiabatic duct. We propose the following simple question: What are the possible solutions for a compressible flow in a duct of length L receiving a certain amount of power Q out ?

This problem has had quite established solutions, but they may appear too restrictive for actual applications. Simple studies can be found in [START_REF] Bruce | Addition of heat to a compressible fluid in motion[END_REF] or in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] : they exhibit analytical solutions for so-called "T 0 -change". Indeed, the heat power applied to the flow will increase (heating) or decrease (cooling) its stagnation temperature, denoted by T 0 . The previously mentioned analyses have shown a specific solution called "choked solution", which corresponds to a Mach number equal to unity, obtained for a maximum heat. Author's analysis by [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] (p.203) suggests that if the heat is larger than this authorized value, the change of flow conditions is done upstream of the duct.

To obtain those solutions, the flow conditions (such as velocity, pressure, temperature) must be specified at the inlet. However, in subsonic regime, it is not possible to keep the quantities fixed as waves will travel up and down the duct, which will change the inlet state ; only the inlet stagnation state conditions are constant. Therefore, in realistic context, the previously proposed analysis cannot be applied to recover solutions of actual flows (without questioning the rigor of the analysis).

In the following study, upstream variations are taken into account to determine solutions, by connecting the inlet of the duct to a tank with fixed stagnation conditions. We examine the case where the heat is not uniformly distributed along the duct, but instead obeys a position-dependent function. In doing so, the analysis can be reduced to a solely heated or cooled flow. The solutions are extended to the "Stiffened-Gas" equation of state, allowing to consider both liquid flows and ideal gases.

Beyond furthering the knowledge on compressible flows, these analytical solutions will help to give preliminary insights on the behavior of a more complex system than what has been developed thus far. Indeed, compressible fluids with high-speed motion have been studied for military, aeronautical and spatial purposes [START_REF] Bruce | Addition of heat to a compressible fluid in motion[END_REF]: as example, cryogenic fluids are used into spatial systems, for propulsion purposes, but also to cool engine and nozzles walls that are subjected to several MW/m 2 of thermal heat fluxes.

Realistic and quantitative numerical methods and models have therefore been developed to perform more complex simulations of non-adiabatic flows [START_REF] Bermúdez | Numerical solution of non-isothermal nonadiabatic flow of real gases in pipelines[END_REF], [START_REF] Schmidmayer | ECOGEN, an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF]. However, those simulations often remain confidential given the fields of research, or too restrictive (limited to a specific geometry, element, regime). Hence, the solutions proposed in this study could be used as future reference solutions for the validation of numerical tools (examples will be presented in Section 5,applied to ECOGEN open-source code by [START_REF] Schmidmayer | ECOGEN, an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF]).

In this paper, we first recall the theoretical elements for flows moving in non-adiabatic ducts. Basic local relations are exhibited in the extended case of "Stiffened-Gas" equation of state. In the second part, a more global solution is developed and provides a full description of the flow behavior in a heated/cooled duct at steady state. The third part is dedicated to the study of some limit cases that leads to the plotting of abacus useful to discriminate subsonic and supersonic flow regimes. In the final part, typical results are presented and analyzed for flows in heated and cooled ducts. Finally, the usefulness of these reference solutions is proved by validating a numerical tool on two practical cases (for subsonic and supersonic behaviors).

2 Model for non-adiabatic single-phase flows

Unsteady model

Although unsteady solutions are not addressed in this study, the governing equations are presented down below as they lay the ground for the problem that will be dealt with in steady conditions. For non-adiabatic compressible inviscid single phase flows, the fluid state at each point is governed by the Euler equations, including heat source term:

     ∂ρ ∂t + div(ρu) = 0 ∂ρu ∂t + div(ρuu + P I) = 0 ∂ρE ∂t + div((ρE + P )u) = δ q (1)
where ρ represents the density, u the velocity vector, P the pressure and E = e + 1 2 u 2 the total specific energy, with e the internal specific energy. δ q represents the received (or lost) thermal power per volume unit (W/m

3 ) of fluid.

The entropy equation is:

∂ρs ∂t + div(ρsu) = δ q T (2)
where T represents the temperature of the fluid and s the specific entropy. System (1) is closed by any convex equation of state (EOS) e = e(P, v). In the present work, the "Stiffened Gas" (SG) EOS is used:

e(P, v) = P + γP ∞ γ -1 v + e ref (3)
This EOS allows to describe the thermodynamic properties of both liquids and gases, as it takes into account both the repulsive effects and the attractive molecular properties that binds the matter together. The EOS parameters γ, P ∞ and e ref are obtained from reference thermodynamic curves, characteristics of the material and transformation under study. See [START_REF] Le Métayer | Élaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF] for details. When P ∞ is null, the SG EOS can be simplified to the Ideal Gas (IG) EOS.

Steady-state analysis

General problem

Let us consider a one-dimensional flow in a rectilinear pipe, defined by a length L, a constant cross-section S and a circumference C, receiving or losing a quantity of heat flux φ over its external surface S ext = C × L. Its inlet state quantities are denoted by subscript (•) in .

We will consider two cases to obtain mono-dimensional steady-state solutions of system [START_REF] Bermúdez | Numerical solution of non-isothermal nonadiabatic flow of real gases in pipelines[END_REF]. The first one (one can call the local solution) is presented in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF]. The inlet quantities are imposed. Hence, we may then ask ourselves what the flow solution is at a distance x from the inlet, denoted by subscript (•) x , as the duct has received a certain amount of heat flux φ. The configuration is depicted in figure 1.

The second case is related to a more global solution, which is supposed to be more realistic. At the inlet, only the stagnation conditions are prescribed. At the outlet, denoted by subscript (•) out , the pressure is prescribed: in the case of the subsonic flow, it is equal to the external one. In the supersonic case, no quantities can be imposed. This configuration will be described in depth later on.

Heat flux 𝜑 𝑃

!" 𝑣 !" 𝑢 !" 0 Inlet 𝑥 ? 𝑃(𝑥) 𝑣(𝑥) 𝑢(𝑥)

Local solution

The following relations are local because they can be used whatever the location of x. Corresponding to the configuration presented in figure 1, integration of system (1) at steady-state between the two points (•) in and (•) x along axis x, results in the following system:

ρ in u in = ρ x u x = ṁs ( 4
)
ρ in u 2 in + P in = ρ x u 2 x + P x (5) ṁs (H x -H in ) = qs,x (6)
where ṁs = ρu = ṁ/S is the mass flow rate per unit area of cross-section, qs = Q/S is the power received per unit area of cross-section and H = h + 1 2 u 2 is the total specific enthalpy. Combining equations ( 4) and ( 5) leads to the equation of the Rayleigh line in the (P, v) plane:

P x = ṁ2 s (v in -v x ) + P in (7)
For a given steady flow, ṁs is constant. Thus for a known inlet state, -ṁ2 s is the slope of its Rayleigh line in the (P, v) plane.

The definition of internal specific enthalpy h = e + P v combined to equations ( 4) and ( 5) leads to equation (6) being rewritten as:

e x (P x , v x ) -e in (P in , v in ) + P in + P x 2 (v x -v in ) = qs,x ṁs (8)
The thermodynamic closure is achieved thanks to a convex EOS e(P, v) (3). Using the SG EOS, an explicit relation between the pressures and specific volumes between states (•) in and (•) x is obtained from (8):

P x = (P in + P ∞ ) (γ + 1)v in -(γ -1)v x (γ + 1)v x -(γ -1)v in + 2(γ -1) (γ + 1)v x -(γ -1)v in qs,x ṁs -P ∞ (9)
This relation describes the Crussard curve in the (P, v) plane. It can be rewritten to express the specific volume at point (•) x :

v x (P in ) = (γ + 1)(P in + P ∞ ) + (γ -1)(P x + P ∞ ) (γ + 1)(P x + P ∞ ) + (γ -1)(P in + P ∞ ) v in + 2(γ -1) (γ + 1)(P x + P ∞ ) + (γ -1)(P in + P ∞ ) qs,x ṁs (10)
The Hugoniot curve is obtained in the limit case for an adiabatic transformation, where qs,x = 0:

P x = (P in + P ∞ ) (γ + 1)v in -(γ -1)v x (γ + 1)v x -(γ -1)v in -P ∞ (11)
The thermodynamic path of fluid between states (•) in and (•) x is represented in figure 2. Each state is constrained by both the Rayleigh line and its respective Crussard curve. This figure clearly shows that for the imposed inlet flow conditions, only one solution is reachable at the location x according the sign of the power received, given by the intersection of the Rayleigh line and the corresponding Crussard curve.

Pressure

Specific volume

Hugoniot curve (in) Crussard curve (x+) Crussard curve (x-) Rayleigh line

) & > 0 ) & < 0 % !" , ' !" % # $ , ' # $ % # % , ' # % Figure 2:
Thermodynamic path for a fluid flowing in a heated/cooled duct in the (P, v) plane. Starting from state (•) in , the fluid can be heated to reach (•) + x state or cooled down to reach (•) - x state.

Local solution in Mach number

Using the SG EOS, Euler equations between (•) in and (•) x can be also rewritten in terms of Mach number. Indeed, analysis of the solution in Mach number allows to obtain information about the global dynamics of the flow. Moreover, this number is quite pertinent and usual for the analysis of compressible flows. The Mach number is defined by M = u/c, with c the speed of sound defined as c 2 = ∂P ∂ρ s . In the case of the SG EOS, the speed of sound reads:

c = γ(P + P ∞ ) ρ (12)
Equation ( 5) is thus expressed as a ratio of pressure:

P in + P ∞ P x + P ∞ = 1 + γM 2 x 1 + γM 2 in ( 13
)
The equation (13) coupled to (4) allows to express the ratio of temperature:

T in T x = M 2 in (1 + γM 2 x ) 2 M 2 x (1 + γM 2 in ) 2 (14)
Taking into account for heat variation, the so-called T 0 -change (change in stagnation temperature) between the inlet (•) in and point (•) x is a direct measure of the amount of heat input. The energy equation ( 6) becomes:

h 0,x -h 0,in = C p (T 0,x -T 0,in ) = qs,x ṁs (15)
The inlet stagnation state, which is retrieved by bringing the flow at rest through an adiabatic process can be expressed as a function of Mach number. This yields the following isentropic relations:

P 0,in + P ∞ P in + P ∞ = 1 + γ -1 2 M 2 in γ γ-1 (16) T 0,in T in = 1 + γ -1 2 M 2 in ( 17
)
The combination of equations ( 14) and (17) allows to find a relation between the Mach number at the inlet (•) in and at point x, through a ratio of stagnation temperatures:

T 0,in T 0,x = M 2 in (1 + γM 2 x ) 2 (1 + γ-1 2 M 2 in ) M 2 x (1 + γM 2 in ) 2 (1 + γ-1 2 M 2 x ) (18
! !,# /! !,$% = 1.1 ! !,# /! !,$% = 1 ! !,# /! !,$% = 0.9 ! !" ! # ! !,# /! !,$% = 1.2 ! !,# /! !,$% = 0.8 0.9 0.8 1.1 1.2
Figure 3: Local solution for M x versus M in : when the inlet is subsonic (M in < 1), the only physical solution is a subsonic outlet (M x < 1). When the inlet is supersonic (M in > 1), the solution can either be supersonic (M x > 1), or subsonic through a shock (M x < 1).

The analysis of these solutions is detailed in literature [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] and show how complex a compressible flow can be. Important results are summarized in figure 3, which illustrates the behavior of the Mach number at point x M x versus the inlet Mach number M in , for fixed values of ratios of stagnation temperatures T 0,x /T 0,in . Obviously, T 0,x /T 0,in = 1 represents adiabatic solutions, where either no change of flow regime happens or a shock wave is a possible solution. For T 0,x /T 0,in ̸ = 1, we are either in the case of heating (T 0,x /T 0,in > 1) or cooling (T 0,x /T 0,in < 1). On these branches, there are two possible values of M x for one value M in . When the inlet is supersonic, the outlet can also be directly supersonic, or subsonic through a shock wave. However, for initially subsonic inlet cases, the only possible physical outcome is also a subsonic flow [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF].

Limitations of local solutions

The previous configuration that corresponds to the local solution is too far of a realistic system and is not selfsufficient to describe the full reality of the flow. Indeed, in a subsonic flow, the inlet continuously re-adapts to the external conditions (such as outlet pressure and heat flux) and cannot be fixed. It means that one cannot impose the inlet conditions and then apply a power heat to calculate the solution at point x: doing so, travelling acoustics waves will change the inlet conditions. Hence, the conditions at the inlet of the duct are no longer known, meaning that the flow along the duct cannot be determined. Notice that for a supersonic flow, the prescription of the inlet state is direct as acoustic perturbations cannot go back upstream : the local solution is then valid. Besides, in actual system, the stagnation conditions are generally known and the flow is then governed by the external pressure.

So, we propose a study that is closer to this configuration, by considering imposed tank conditions instead of inlet conditions. In the tank, the velocity is equal to zero while the pressure P 0 and temperature T 0 are given data. It also means that at the flow regime in the section connected to the tank is necessarily subsonic, that is to say M in < 1. The possible solutions for Mach number are then depicted in figure 4, which is a restriction of the solutions plotted in figure 3: on one hand, locally heating increases the flow velocity (curve slope over 1) up until the sonic point ; on the other hand, cooling decreases the flow velocity (curve slope under 1).

In the following analysis, a complete description of the flow is proposed. Stemming from the previous part, steady relations will be established in two parts: the first set of relations are the local equations adapted in the case of interest to obtain the flow solution along the duct. The second set of relations is between the inlet of the duct and the tank at the stagnation conditions; notice that the inlet flow conditions are unknown in this configuration. Coupling both systems will allow to define the flow characteristics along the duct, for given stagnation conditions, an outlet pressure and an applied heat flux function.

Min

! !,# /! !,$% = 0.9

! !" ! # ! !,# /! !,$% = 0.8 ! !,# /! !,$% = 1 ! !,# /! !,$% = 1.1 ! !,# /! !,$% = 1.2
Figure 4: Local solution for M x versus M in : focus on the subsonic inlet area (M in < 1), which can only lead to a subsonic or sonic outlet (M x ≤ 1), whatever the heat-flux applied.

3 Global solution in a duct

Observations

Before delving into a mathematical description of the considered flows, we present thereafter some typical qualitative results of flows in a duct. Some flow behaviors are easily deduced from isentropic or adiabatic flow analysis, or from the basic results [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] depicted in figure 3, while the more complex ones are presented in light of the subsequent analysis. For now, the aim is simply to highlight the key parameters of the flow, in order to familiarize the reader with its behavior in this new configuration.

We describe a flow in a duct subjected to a total amount of power Q = Q i , where Q i is a local applied amount of heating/cooling power, which can be positive or negative. This study is inspired by the study of adiabatic flows in nozzles: in such flows, the area of the cross section of the nozzle is of prime importance, through the ratio of the area over the critical area. In the present work, the local amount of heat will play this role. Nevertheless, the characteristic of the configurations (area for nozzle, heating for the duct) is not sufficient to determine flows conditions: they are driven by the external pressure P out and more accurately by the pressure ratio defined as Π = P out /P 0 . Because the tank conditions are given fixed data, the only way to create a flow through the nozzle or the duct is to vary the outlet pressure P out , like one would do opening or closing a valve of a vacuum pump. To generate a flow from the tank to the outlet of the duct, the pressure is chosen accordingly and Π is necessarily less than 1.

Isentropic flow

When the duct is not subjected to any heat flux, the flow is defined by only one state (constant along the duct) which can be recovered through equations ( 16) and (17): its pressure is equal to the external pressure P out .

Decreasing Π = P out /P 0 leads to an increase of the flow velocity. As the inlet of the duct is connected to the tank, isentropic relations govern the flow between both, hence why the inlet state can only be subsonic. Therefore, the velocity increases to reach the sonic value, as the pressure ratio is decreased: then, a sonic choking of the flow would occur at the outlet and the downstream pressure would no longer be able to affect the flow inside the duct, remaining at the sonic state.

Heated flow

Now the duct is subjected to a fixed heating power Q heat > 0. Because of the heat exchange, the flow cannot be constant in the duct: indeed, the amount of heat received at a location x is less than the heat received at the location

Mach number Distance M = 1 Isentropic (Q heat = 0W) Heated (Q heat > 0W) Pressure Distance Isentropic (Q heat = 0W) Heated (Q heat > 0W)
Figure 5: Qualitative behavior of a flow within a duct connected to a tank and an outlet pressure with a subsonic inlet, for different heating powers: Mach number (left) and pressure (right).

x + dx. Therefore, the inlet will adapt to the outlet pressure in the way described by figure 4: the inlet Mach number M in (always subsonic) is necessary less than M x at any position x in the duct. So, for a heating power, the flow will accelerate all along the duct. Moreover, the Mach number cannot become supersonic to maintain the physicality of the system, as figure 4 shows that with a subsonic inlet, the outlet must also be subsonic. Considering the pressure, its evolution along the duct is on the opposite of the Mach number, meaning that a decrease will be observed between the inlet and the outlet cross sections. The behavior of these quantities are illustrated in figure 5, where Mach and pressure profile are plotted for a given heating powers.

Cooled flow

The flow is now cooled down by a fixed cooling power Q cool < 0. The pressure and Mach profiles evolve in the opposite way to the heated case: pressure is increasing between inlet and exit section, while the flow is decelerated.

Hence, when Π is decreased, the inlet pressure decreases and the inlet Mach number rises to become supersonic: this is not part of our study as it will not respect the hypothesis of a subsonic inlet.

Obviously, the upper value of the ratio Π is bounded by unity ; this short analysis indicates that some lower limits may exist and could be exhibited by deeper calculations.

Heated and cooled flow

In this work and in order to complete our observations, we combine both configurations: after the heated subsection receiving Q heat , a subsection subjected to a fixed cooling power Q cool is added. For a rather high pressure ratio Π (but still bounded by 1), the flow remains subsonic: the heated part accelerates, whereas the cooled part decelerates. When the outlet pressure is decreased, the flow overall accelerates as the end part of the heated subsection is brought closer to the sonic point.

Once Π gets below a certain value, the flow reaches the sonic point at the end of the heated part. Then, the flow in the cooled part can become subsonic or supersonic because both solutions are possible. It suggests that this bifurcation is governed by specific pressure ratios called critical pressure ratios.

Cooled and heated flow

This configuration has no theoretical interest and the whole solution is provided by the analysis of figure 4 that shows the flow will decelerate and accelerate while the outlet section will remain subsonic.

This short description of flow receiving a power shows the outlines of the study that will be presented hereafter: the solutions for either heated or cooled flows are well documented in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF]. Cooled and heated duct flows solutions are easily deduced from the previous observations. On the contrary, the solutions for heated and cooled duct flows seem to exhibit non trivial solutions. The next part of the study will focus on this configuration.

System description

We propose to investigate an uneven heating distribution on the duct, which consists in an extension of the work done in [START_REF] Fathalli | Liquid-vapor flows in heated ducts: Reference solutions[END_REF]. This will allow to describe as many cases as possible, as we have seen that it is now possible to reach either a subsonic outlet or a supersonic outlet with a subsonic inlet in such configuration.

! 0 Inlet ! !"#$ Point of heat flux variation # Outlet Imposed outlet state $ %&$ $ '( % '( & '( Tank state $ ) , ( ) Heated subsection ) !"#$ Cooled subsection ) %&$ 𝑥 !"#$ 𝜑 %&&'
Figure 6: Tank flowing in a heated and cooled pipe To study the flow, we consider a duct connected to a tank (denoted by subscript (•) 0 ) whose state is supposed known (pressure and temperature), and an outlet pressure imposed by external conditions. The whole duct is governed by non-adiabatic relationships as it receives or loses a certain amount of heat flux, whereas the tank and inlet states are coupled by isentropic relationships. The system is depicted in figure 6.

ϕ heat ϕ cool x heat Linear heat flux ϕ (W/m 2 ) Distance Q heat Q out x heat Power received Q (W) Distance

Duct relations

The duct under study is divided in two sub-sections, separated by the point of variation of flux distribution (defining the subscript (•) heat ): a heated subsection between x = 0 and x = x heat , followed by a cooled subsection between x = x heat and x = L. The local relations exposed previously are suitable to describe the behavior of the flow in the duct.

Along axis x, a function of heat flux φ(x) (W.m -2 ) is distributed on the external surface of the duct:

φ(x) = φ heat > 0 if x ≤ x heat φ cool < 0 otherwise
The power Q(x) =

x 0 φ(η) C dη (W ) applied to the flowing fluid is determined by the distribution of φ(x) (figure 7):

Q(x) = C x φ heat if x ≤ x heat C [x heat φ heat + (x -x heat ) φ cool ] otherwise
The power received by the flow on the entire heated subsection is Q heat = C x heat φ heat , the power lost by the flow on the entire cooled subsection is Q cool = C (L -x heat ) φ cool and the power received by the flow on the entire duct is

Q out = Q heat + Q cool . If Q out > 0,
the duct has received more heat than it has lost. Conversely, if Q out < 0, the duct has lost more heat than it has received. Finally, if Q out = 0, we can consider the outlet flow adiabatic.

The power received per unit area of cross-section is defined as : qs (x) = Q(x)/S (W.m -2 ). As we apply a function of heat flux qs (x) ∝ Q(x) on each point x ≡ (•) x of the duct, equations ( 7) and (9) govern the flow between the inlet and any point of the duct. The applied heat flux can be directly correlated to the notion of T 0 -change through equation ( 15):

T 0,x T 0,in = qs,x ṁs C p T 0,in + 1 (19)

Tank-inlet relations

The closure of the system is thus done by connecting the inlet of the duct to a tank. The inlet stagnation state is recovered by slowing down the flow to rest through an isentropic process and is equivalent to the tank state. This means that P 0,in ≡ P 0 and T 0,in ≡ T 0 . On the tank-inlet side, the only available relations are: H in = H 0 , s in = s 0 where H is the total specific enthalpy, and s is the specific entropy. Even if the supersonic solution is mathematically acceptable, only the subsonic solution is sought to be consistent with an actual tank. The first condition H in = H 0 allows us to establish a relation linking the mass flow rate and the inlet pressure:

ṁs (P in ) = 2(h 0 -h is in ) v is in (20)
v is in (P in ) and h is in (P in ) are respectively the isentropic specific volume and internal enthalpy provided using the second condition s in = s 0 . Both are a function of pressure and depend on the equation of state. The isentropic relationships are written using the SG EOS:

v is in (P in ) = v 0 P 0 + P ∞ P in + P ∞ 1 γ (21) h is in (P in ) = γv 0 (P 0 + P ∞ ) γ -1 P 0 + P ∞ P in + P ∞ 1-γ γ + e ref (22)
In these equations, the quantities related to the inlet cross section, such as P in , are unknown.

Reference solutions : determining flow behavior

The studies done on heated ducts in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] and [START_REF] Fathalli | Liquid-vapor flows in heated ducts: Reference solutions[END_REF] completed by the previous observations have shown that once the duct is connected to a tank, as its inlet must be subsonic due to the coupling, it can only reach the sonic point at most and cannot become supersonic. However, when appending a cooled subsection to the duct, we have seen that under a specific pressure ratio, the sonic state is in fact reachable and the cooled subsection can lead to a supersonic state at the outlet flow.

In this part, we find out by analytical calculations what are the conditions required to determine the flow regime.

As introduced in the previous observations, this type of description may remind adiabatic nozzle flows, with the following analogy done:

• The converging part of the nozzle would be the heated subsection of the rectilinear duct.

• The diverging part would be the cooled subsection.

• The fixed geometry of the nozzle (ie: the area ratio between the throat and the outlet A * /A(x)) is viewed similarly to the fixed power function Q(x) applied to the duct. • The throat of the nozzle would be the point of heat flux change (x heat ). We will then apply this analogy and we first hypothesize the existence of the same theoretical critical pressure ratios presented for example in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF]: Π 1 , Π 2 and Π 3 (table 1). They divide the flow into different specific regimes: fully subsonic, presence of a sonic point (choked flow), supersonic, shocked... Identifying those critical pressure ratios will allow us to categorize and properly determine the flow behavior. It must be noted that despite not yet having noticed a steady shock wave in the cooled subsection (Π 1 > Π ≥ Π 2 ), its existence is alleged and may be proven true or untrue later on.

The different values of critical Π are now sought by solving the system of equations composed of the local solution and the tank-inlet relations.

Fully subsonic flow

The flow is supposed to be fully subsonic when Π > Π 1 (M < 1 throughout the whole duct). Given that we do not reach sonic state at the end of the heated subsection, the duct also remains subsonic in the cooled subsection. Therefore, information coming from further than the point of heat flux variation (such as the outlet state) can impact the inlet.

Π = 1 No flow 1 > Π > Π 1
Fully subsonic flow Π = Π 1 Subsonic inlet, choked flow at the end of the heated subsection, subsonic outlet

Π 1 > Π > Π 2
Subsonic inlet, choked flow at the end of the heated subsection, steady shock wave in cooled subsection, subsonic outlet Π = Π 2 Subsonic inlet, choked flow at the end of the heated subsection, steady shock wave at the outlet cross section

Π 2 > Π > Π 3
Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Shock waves outside of the duct) Π = Π 3 Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Adapted flow outside of the duct) Π 3 > Π Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Expansion waves outside of the duct)

Table 1: Nozzle pressure ratios definition applied to heated and cooled duct.

Inlet state : Firstly, we seek for the determination of the inlet state, as it will constraint the rest of the flow. In the subsonic case, the outlet pressure P out is a known constant and pressure information coming from downstream has an impact on the inlet flow of the duct. Equation ( 7) is evaluated between (•) in and (•) out to numerically determine the inlet pressure P in under the following implicit form:

f (P in ) = ṁ2 s (v out -v is in ) + P out -P in = 0 ( 23
)
On one hand, the isentropic relations between the tank and the inlet allow to calculate the specific volume at the inlet v is in (P in ) with equation ( 21) and the mass flow rate ṁs (P in ) with equation (20). On the other hand, the local relations evaluated between the inlet and the outlet allow the determination of the specific volume at the outlet v out (P in ), with equation (10), as the heat flux Q out ∝ qs,out is also known.

The value of the inlet Mach number M in (P in ) is determined by the isentropic relation between the tank and the inlet, through equation ( 16).

Duct state : Once the inlet state is known, we now seek to characterize the flow at any point x in the duct. This will be done by using the local solutions of the duct. The ratio of stagnation temperature at point x, denoted T 0,x /T 0,in , is obtained through equation ( 19). This equation is coupled with equation (18), and both express the so-called T 0 -change. This yields the following equation, which is a biquadratic function of M x :

f (M x ) = M 4 x βB -γ 2 A T 0,x T 0,in + M 2 x B -2γ T 0,x T 0,in -A T 0,x T 0,in = 0,      β = γ-1 2 A(P in ) = M 2 in 1 + βM 2 in B(P in ) = (1 + γM 2 in ) 2 (24)
See appendix Solving M x for solutions of f (M x ) = 0: there are 4 solutions, 2 of which are negative and thus non-physical, and 2 other which are positive, thus physical, but either subsonic (< 1) or supersonic (> 1). As we are studying the fully subsonic case, the lowest positive root of the equation (0 < M sub x < 1) represents the subsonic Mach number of this case's flow.

Given the Mach number at any point x ∈ [0; L], the pressure P x is computed from equation ( 13), with M in determined from equation (16):

P x (P in ) = (P in + P ∞ ) 1 + γM 2 in 1 + γM 2 x -P ∞ (25)

Choked flow

The flow is choked once Π ≤ Π 1 : a sonic state is reached in the duct. This state can only be observed at the x heat location, which is the end of the heated subsection. Given that the flow is sonic at x heat , no information can come from downstream of this location: it acts a barrier for any backwards travelling wave. Then, the behavior of the flow within the cooled subsection can follow different paths, depending on the critical pressure ratios presented in table 1. Our goal now is to determine these hypothetical values.

Inlet state : Given prior considerations, to determine the inlet state in the case of a choked flow, we focus on the the coupled system tank-inlet and inlet-x heat instead of inlet-outlet. Indeed, the flow is sonic at x heat and that will not change. Equation ( 7) is therefore evaluated between (•) * in and x heat ≡ (•) heat to numerically determine P * in , where (•) * in denotes the inlet state for which the flow ends up choked:

f (P * in ) = ṁ2 s (v heat -v is in ) + P heat -P * in = 0 (26)
In the same fashion as the subsonic flow, the relations between the tank and the inlet allow to compute the mass flow rate ṁs (P * in ) with equation ( 20) and the specific volume v is in (P * in ) with equation ( 21). The local relations evaluated between (•) * in and (•) heat lead to the determination of the specific volume v heat (P * in ), through equation ( 10), given that Q heat ∝ qs,heat is known and fixed.

The remaining unknown of equation ( 26) is the pressure at the end of the heated subsection, P heat . Our hypothesis that the flow is choked means that it reaches the sonic state and thus, M heat = 1. The expression connecting the Mach number and the pressure is given by equation (13), evaluated locally between (•) * in and (•) heat :

P heat (P * in ) = (P * in + P ∞ ) 1 + γM 2 in 1 + γ -P ∞ (27)
The value of the inlet Mach number M in (P * in ) is determined by the isentropic relation between the tank and the inlet, via equation (16).

Duct state : Now, the inlet state is known and from it we seek to characterize the flow at any point x of the duct. The determination of the stagnation temperature T 0,x through equation (15), coupled with the ratio of stagnation temperature in the duct T 0,x /T 0,in with equation (18), yields the same polynomial of order 4 presented by equation (24), function of Mach number at point x denoted M x (P * in ). See appendix for the resolution of the polynomial.

In the heated subsection (x ∈ [0; x heat [), the Mach number is necessarily subsonic. Thus, its value the lowest positive root of equation (24) (0 < M sub x (P * in ) ≤ 1). We also know that the sonic point is reached at (x = x heat ). In the cooled subsection however (x ∈]x heat ; L]), once the flow is choked, there are two possible regimes depending on the pressure ratio that is imposed: subsonic (therefore still 0 < M sub x (P * in ) < 1), or supersonic (therefore M sup x (P * in ) > 1). They respectively correspond to the lowest positive root and highest positive root of equation ( 24). This behavior can be seen in figure 3: once the beginning of the cooled subsection is sonic, there are two possible states for T 0,x /T 0,in < 1: one where M x > 1 and one where M x < 1.

Depending on Mach number value M x (P * in ), the pressure for each possible paths is evaluated from equation (25). Again, the value of Mach number at inlet M in (P * in ) is determined from equation (16). Figure 8 depicts the possible critical paths described above. Critical outlet states : When the flow is choked (sonic point at the end of the heated subsection), one can see that the outlet pressure P out is not needed anymore. Indeed, P out is absent from the system formed by equations (26), ( 27) and (24) to determine M x . In this configuration, only the heat flux function and the tank parameters condition the flow through the variation of stagnation temperature described with equation (19).

The critical pressures for choked flow and then subsonic outlet and supersonic outlet are evaluated from equation (25) at point (•) out . Relating back to table 1, they respectively allow to determine Π 1 = P sub out /P 0 and Π 3 = P sup out /P 0 , both functions of P * in . Hugoniot curve "in" Crussard curve "heat" Crussard curve "out" Rayleigh line

𝑀 !" < 1 𝑀 #$% < 1 𝑀 #$% > 1 𝑀 &'(% = 1 (a) (b) (c)'' (c)'
Figure 9: Thermodynamic path for a heated/cooled choked flow in the (P, v) plane, with corresponding Mach number values. Starting at a subsonic state on the Hugoniot curve "in" (a), the fluid receives a certain amount of heat.

If the pressure ratio is critical in the given configuration, the flow accelerates until it reaches sonic point at "heat" state. The corresponding Crussard curve "heat" is tangent to the Rayleigh line of the flow at the sonic point (b).

Then, as the fluid is cooled down to reach outlet state, it moves to another Crussard curve "out" (c): the outlet state winds up at the intersection of the Rayleigh line and the outlet Crussard curve. Depending on the Π, the outlet state is either the subsonic state (c'), or the supersonic state (c") ; there can be no intermediate steady shocked state between the latters.

To determine Π 2 = P shock out /P 0 , we assume an adiabatic shock wave located in the outlet section of the duct, starting from the supersonic outlet state. Being adiabatic, the shock wave state remains on the Hugoniot curve of the supersonic outlet state and must end up at the intersection with the Rayleigh line, as in a choked flow, its slope -ṁ2

s is fixed by heated duct properties. Constrained by these two curves, the outlet shock wave can only bring the supersonic outlet state to the subsonic outlet state. This means that the shocked outlet state is equivalent to the subsonic outlet state determined previously with Π 1 (P shock out = P sub out ⇒ Π 2 = Π 1 ). The thermodynamic path is displayed in figure 9. Thus recalling table 1, given that Π 2 = Π 1 , we have demonstrated that there can be no steady shock wave in the cooled subsection. If any part of this subsection has reached a supersonic state, in a steady state the whole subsection will become supersonic. In the choked flow conditions, given a specific Q heat and inlet state, only two outlet states are possible.

The updated pressure ratio table is displayed in table 2. At this point, knowing the heating and cooling powers Q heat and Q cool , the critical pressure ratios are known. We can then determine the flow behavior within the duct for a given pressure ratio Π.

Circling back to the nozzle studies analogy, only knowing the ratio of exit and throat areas is necessary (therefore only one data), whereas here a couple of data is needed (any combination of Q heat and Q cool ).

Π = 1

No flow 1 > Π > Π 1

Fully subsonic flow Π = Π 1 Subsonic inlet, choked flow at the end of the heated subsection, subsonic outlet

Π 1 > Π > Π 3
Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Shock waves outside of the duct) Π = Π 3 Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Adapted flow outside of the duct) Π 3 > Π Subsonic inlet, choked flow at the end of the heated subsection, supersonic outlet (Expansion waves outside of the duct) 4 Critical pressure ratios behavior

In the following section, we study the critical pressure ratios: the goal is to determine the limits of the possible regimes, depending on the power function Q(x) that the duct is subjected to, and provide analytical and physical limitations of our study framework. We especially take a look at the behavior of Π 1 , as it marks the threshold between a subsonic outlet and a supersonic outlet.

Asymptotic analysis -Limit cases

Isentropic flow

To begin with, when the flow in the duct is neither heated nor cooled (thus isentropic), we deduce that the state within the duct is constant. As P (x) = P out ∀x, the Mach number through the whole duct is determined by equation ( 16).

If the flow actually reaches sonic state, then P * in = P out and M in = M out = 1. Equation ( 16) evaluated in this specific configuration yields the unheated isentropic sonic limit case:

Π(Q(x) = 0, ∀x) = P out P 0 = P 0 + P ∞ P 0 2 γ + 1 γ γ-1 - P ∞ P 0 ≡ Π is (28) In the case of ideal gases, Π is = 2 γ + 1 γ γ-1 .

Heated flow

If the duct is now heated, the flow state is no longer constant. If it reaches sonic point at its outlet, it means that M out = 1. The outlet pressure P out can be expressed by evaluating equation ( 13) between (•) in and (•) out :

P out = (P * in + P ∞ ) 1 + γM 2 in 1 + γ -P ∞
If the heating power is largely increased, the inlet pressure P * in rises to reach tank pressure P 0 and the pressure along the duct decreases as the flow is heated, to reach the outlet state determined by the pressure ratio. This can be seen through figure 5. A quasi zero pressure gradient between the tank and the inlet implies that the flow is almost at rest at the inlet. The inlet pressure cannot get higher than the tank pressure, as the direction of the flow would then be reversed and our study hypotheses would no longer be valid. Through isentropic equation ( 16) between the tank and the inlet, we deduce the inlet Mach number M in → 0. In that configuration, we can express the theoretical value of lowest critical pressure ratio Π h :

Π(Q → +∞) = P out P 0 = P 0 -γP ∞ (1 + γ)P 0 ≡ Π h (29)
In the case of ideal gases, Π h = 1/(1 + γ). However, for SG, this limit depends on the stagnation pressure P 0 (note for both cases, the limit case is Π h < 1). For a given heating power, once the duct reaches sonic state, it has been shown that the mass flow rate is fixed. Increasing the heating power leads to a decrease of critical mass flow rate ( ṁ * s ≡ ṁs (P * in )), as shown by figure 10. Figure 11 displays the behavior of Π 1 versus the heat flux qs,heat ∝ Q heat in order to remove geometrical dependency, in the following conditions:

• Gas: γ = 1.358, C v = 1247, e ref = 1.97 • 10 7 J/kg, P ∞ = 0 Pa • Tank state: P 0 = 1.5 Bar, T 0 = 394 K When the heating power decreases, Π 1 converges towards the singular value Π is presented in equation ( 28). It is the pressure ratio needed to obtain a choked flow without any heating in a given configuration : if Π > Π is , the flow can never reach the sonic point, as heating will only decrease the Mach number. On the other hand, if the heating power increases, the value of Π 1 decreases and tends towards the theoretical pressure ratio Π h introduced with equation (29). This corresponds to the lowest pressure that can be reached at the sonic point when increasing the heating power. At a certain point, increasing heating becomes pointless, as the critical pressure ratio Π 1 is bounded from below by the value of Π h .

Heated and cooled flow

Heated and cooled flows are here considered. We assume that the flow is choked and therefore reaches sonic point at the end of the heated subsection, thus M heat = 1 and P heat is expressed by equation ( 27).

On one hand, if the cooling power is negligible in front of the heating power, the flow behavior does not vary much after entering the cooled subsection, as one could consider the cooled subsection close to adiabatic. The pressure

Q heat = 1 kW Q heat = 10 kW Q heat = 100 kW Q heat = 1 MW Q heat = 5 MW Q heat = 10 MW π is π h
Figure 12: Critical pressure ratio Π 1 evolution depending on the value of ratio of heating over cooling power, for different values of heating power. For a given heating power and power ratio, if Π > Π 1 , then the flow is fully subsonic. If Π < Π 1 , then it is choked and supersonic in the cooled subsection.

ratio needed to bring the end of the heated subsection to the sonic point follows the curve presented by figure 11. If the heating power is high enough, the pressure ratio tends towards the value of the critical pressure ratio of a duct that is solely heated: in that case, P out → P heat and Π 1 → Π h . On the other hand, when the cooling power applied to the flow is much larger than the heating power, and thus the flow is overall cooled down at the outlet, the outlet state pressure P out caused by a sonic flow increases compared to the inlet pressure (see figure 2). However, to keep our hypothesis of flow going from a tank to an outlet, the pressure ratio Π = P out /P 0 must be less than 1 ; otherwise, the direction of the flow will be reversed. Therefore, we pose the following inequality, which will yield a condition on the lowest applicable cooling power Q min cool : Π 1 = P out /P 0 < 1. To solve it, we must determine the outlet pressure P out (P * in ) by evaluating equation (25) at the outlet. It depends on the inlet state, which is known and fixed by the properties of the heated subsection, as we suppose a choked flow. The inlet pressure P * in is determined from equation (26) and the inlet Mach number M in (P * in ) is determined from equation ( 16).

The outlet pressure also depends on the outlet Mach number M sub out (P * in ), which can be evaluated by finding the lowest positive root of equation (24) at point (•) out . M sub out depends on the outlet stagnation temperature T 0,out , which depends on the power applied Q out = Q heat + Q cool . Therefore, after retrieving T min 0,out , the lowest possible cooling power that keeps a flow going from the tank to the outlet can be easily determined from equation (19):

Q min cool > T min 0,out T 0,in -1 ṁs C p T 0,in S -Q heat (30)
Beyond the apparent simplicity of this equation, the calculation of T min 0,out is not trivial and equation details can be found in appendix Solving Q min cool . As we have seen with the above limit cases, when Q heat ≫ |Q cool |, the outlet pressure tends towards P heat < P * in . When Q heat ≪ |Q cool |, the outlet pressure rises towards P 0 (at most, given our study hypothesis). Thus, when the system under study is both heated and then cooled down, Π 1 ∈]Π h ; 1[. 4.2 Critical pressure ratios dependence regarding Q heat and Q cool (abacus) Now, we consider a nondescript combination of heating and cooling power. Figure 12 shows the evolution of the critical pressure ratio Π 1 against the ratio of heating over cooling power |Q heat /Q cool |. For a fixed power ratio, there are multiple possible critical pressure ratios. Therefore, either the heating power or cooling power must be a fixed parameter as well. For example, for a power ratio

|Q heat /Q cool | = 1, if Q heat = 5 M W , then the outlet flow becomes supersonic for Π < Π 1 = 0.75.
Let us consider a fixed power ratio

|Q heat /Q cool |, where either Q heat ≫ |Q cool | or Q heat ≪ |Q cool | .
On one hand, if the heating power Q heat decreases, the flow critical pressure ratios converge towards the singular value Π is presented in equation (28). It is the pressure ratio needed to obtain a choked flow without any heating in a given configuration : if Π > Π is , the flow can never reach the sonic point, as heating will only decrease the Mach number.

On the other hand, if the heating power Q heat increases, the flow critical pressure ratios convergence point decreases and tends towards the theoretical pressure ratio Π h introduced with equation (29). This is the lowest pressure that can be reached at the sonic point when increasing the heating power. At a certain point, increasing heating becomes pointless, as the critical pressure ratio Π 1 is bounded from below by the value of Π h .

Results

We present some results that illustrate the flow behaviour in various conditions: these results are obtained by the analytical solutions proposed in this study. They are also the answer to the question presented in the introduction: What are the possible solutions for a flow in a duct of length L receiving a certain amount of power Q out ? Then, we also present comparisons with numerical results provided by the numerical solution of the Euler equations system (1), thanks to a code based on finite volume method.

Parameters study

For every test case that is presented, we consider a duct of length L = 1 m heated for x ∈ [0; x heat ] and cooled for x ∈]x heat ; L], with a total volume of V = 1/(4π) m 3 . For a cylindrical duct, that leads to a radius of R = 1/(2π) m, a cross-section S = πR 2 = 1/(4π) m 2 and a perimeter P = 2πR = 1 m.

The tank conditions are also held constant in the different cases : P 0 = 1.5 Bar, T 0 = 394 K.

The fluid under study is a gas obeying the IG EOS with the parameters γ = 1.358, C v = 1247, e ref = 1.97 • 10 7 J/kg, P ∞ = 0 (which is particular case of the SG EOS). The heat capacity at constant volume C v is required to extract the temperature given by relationship e = C v T .

The value of the heating power is Q heat = 200 kW. The cooling power would be a parametric term to get some various conditions through

Q out = Q heat + Q cool .
The corresponding value of Q min cool is -5.45 MW. All those parameters define the system under consideration. To calculate various flows, a dynamic parameter must be changed to differentiate these cases, which is what we will study thereafter.

Variation of pressure ratio Π

The dynamic parameter of this study is obviously the outlet pressure through the ratio Π: we take a look at the impact of the pressure ratio on the behavior of the flow.

In this test case, we pose Q cool = -250 kW, meaning that Q out = -50 kW: the flow is overall more cooled than heated. The critical pressure ratios can first be calculated by the method presented in Section 3.3, which yields the following values: Π 1 = 0.6313, Π 3 = 0.4303. The whole solution is then calculated between the inlet and the outlet. The particular solutions corresponding to Π 1 and Π 3 are labelled in figure 13. In this figure, the evolution of the pressure, the Mach number, the temperature and the flow velocity is depicted, for various values of Π ∈ [0.60; 0.73].

Several solutions are represented for Π ≥ Π 1 : the mass flow rate ṁs (P in ) varies as the state of the heated subsection still depends on the variations of P out (as depicted on figure 14). In this range, the pressure and the Mach number follow a coherent pattern. The Mach number increases up to x heat and then decreases in the cooled part of the duct, and the pressure follows the reverse pattern. Extremum values are reached at x heat for Π = Π 1 .

Once Π < Π 1 and the flow is choked, ṁs (P * in ) is fixed by heated subsection, which is determined and cannot change (see figure 14). There is no steady shock wave in the cooled subsection and every supersonic solution is overlapping (see red curve in figure 13).

Let's take a look at the evolution of the temperature, which is quite interesting and relates an example of possible particular behaviors within compressible flows. For greater values of Π, in the heated part, the temperature increase is monotonous up until x heat and then, in the same manner, decreases in the cooled part. This result seems to be obvious: when a fluid is heated, its temperature increases. However, for lower values of Π, the temperature decreases in the heated part of the duct and increases in the cooled part (for Π > Π 1 ). This result is a known particularity of non-adiabatic compressible flows: the flow follows this pattern if 1/ √ γ < M < 1 [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF] (see for example purple curve in figure 13).

Variation of outlet heat power Q out

We now fix the pressure ratio Π and take a look at the impact of variable power.

Q out = 0W Q out > 0W Q out < 0W
Figure 15: Pressure (top row) and Mach number (bottom row) evolution in a heated and cooled duct for different outlet states, with a fixed heating power. For Π > Π 1 (left), the outlet state is constrained but the whole duct re-adapts. For Π < Π 1 (right), the heated subsection state is fixed but the cooled subsection re-adapts.

As Q heat is fixed, modifying Q out is equivalent to modifying the cooling power Q cool . When Q out > 0, then |Q cool | < Q heat : therefore, there is overall more heating than cooling power on the duct. When Q out < 0, then |Q cool | > Q heat : then, there is more cooling power.

Figure 15 underlines the difference in flow behaviour for different values of outlet power Q out , with a fixed Q heat , depending on the flow regime. Overall, in the case of a non-adiabatic outlet, the pressure varies compared to an adiabatic outlet. Only two representative flow cases are considered: the first case is for a pressure ratio Π > Π 1 , ensuring a fully subsonic flow, and the second case is for Π < Π 1 .

For a subsonic flow (Π > Π 1 ), in the case of an adiabatic outlet, there is symmetry between inlet and outlet, which is lost for a non-adiabatic outlet. If the outlet state is no longer adiabatic, as the outlet pressure P out conditions the flow, the full duct re-adapts depending on the value of Q out . Within the duct, as the flow is heated and cooled down but remains subsonic, it follows the behavior described by figure 4. Between the inlet and the outlet, when Q out > 0, the pressure generally decreases (P in > P out ) and the Mach number generally increases (M in < M out ). On the opposite, for Q out < 0, the pressure generally increases (P in < P out ) and the Mach number generally decreases (M in > M out ).

For a supersonic outlet flow (Π < Π 1 ), the heated subsection is conditioned by P heat , which itself is determined by the fixed heating power Q heat . Thus, only modifying Q out cannot impact the heated subsection and the inlet state, but does change the cooled subsection state. When Q out > 0, P out decreases less than when Q out < 0, and the Mach number increases less. This can be explained by figure 3: now that the supersonic state is reached in the cooled subsection, heating leads to a decrease of Mach number and cooling leads to an increase of it.

There would be little information in studying the effects of the value of Q heat , which would only change the levels without changing the qualitative evolution of quantities, compared to the case presented in this section. We can simply underline the following points: depending on if the flow is choked or not, an increase of heating power will vary the inlet Mach number and the inlet pressure, as described in figure 4. It will also decrease the value of Π 1 , as figure 12 shows.

P Π > Π 1 Π < Π 1 Π > Π 1 Π < Π 1 NX = 100
0.106% 0.266% 0.061% 0.739% NX = 1000 0.057% 0.085% 0.006% 0.222% Table 3: Maximum values of relative pressure error depending on the level of meshing NX and the pressure ratio Π, for numerical simulations of a gas and a liquid.

Example of validation for CFD codes

These analytical results are used to validate CFD codes simulations of compressible single-phase flows in non-adiabatic ducts, whether they are being solely heated, cooled down, or both. In this subsection, we use the CFD code described in [START_REF] Schmidmayer | ECOGEN, an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF], which uses the [START_REF] Godounov | Résolution numérique des problemes multidimensionnels de la dynamique des gaz[END_REF] method, in the following configuration:

• Tank state: For different levels of meshing (NX being the number of mesh cells), figure 16 and figure 18 compare the numerical code results with analytical fully subsonic results, respectively for the gas and the liquid. Figure 17 and figure 19 compare them with analytical supersonic outlet results, also respectively for the gas and the liquid.

P 0 = 1.5 Bar, T 0 = 394 K • Geometry: V = 1/4π m 3 , x heat = L/2 • Power: Q heat = 200 kW, Q cool = -
We define the maximum relative error on variable y as a function of meshing as follow: ε y = max |ynum-y analytical | y analytical . Table 3 displays the different maximum values of relative pressure error ε P : corresponding to the previously mentioned figures, the finer the mesh, the closer the numerical solution gets to the analytical one. Moreover, when the solution is fully subsonic, the numerical solution is closer to the analytical one.

Conclusion

Reference analytical or quasi-analytical solutions in heated and cooled rectilinear ducts have been addressed for compressible inviscid single-phase flows. They were developed by using the "Stiffened Gas" equation of state, allowing to describe gases or liquids.

Following the study done in a heated duct, a cooled subsection was added downstream to study its impact on the behavior of the flow. An analogy to nozzle behaviors was done through critical pressure ratios definition: the method to determine them was developed, depending on the value of the heat flux function applied to the duct and the tank conditions used to prescribe the inlet state.

The definition of two possible states was provided. When the pressure ratio between the outlet of the duct and the tank is greater than a specific critical pressure ratio, the mixture remains fully subsonic. When the pressure ratio becomes lower, the heated subsection reaches a sonic point at its end and the cooled subsection triggers a unique supersonic flow. The thermodynamic analysis showed that there cannot be a steady shock-wave in such flows.

The value of the critical pressure ratio demarcating the shift between the fully subsonic regime or the supersonic outlet regime was also studied: departing from the nozzle analogy where only a single ratio of areas is needed to determine the regime of the flow, both values of powers must be provided for heated and cooled ducts. Limit cases for this critical pressure ratios were developed, for cases without heat, extreme heating power or extreme cooling power.

Finally, results were provided, in good agreement with the analytical study. Varying the pressure ratio when it remains higher than the critical one leads to varying fully subsonic solutions with different mass flow rates. When it becomes lower, the cooled duct is fully supersonic and the solution is unique: the mass flow rate is fixed by the heated subsection and no more variation of solution happens. Locally, even with a more complex two-step function of power, varying the power applied results in predictable change of behavior described in [START_REF] Shapiro | The dynamics and thermodynamics of compressible fluid flow[END_REF].

The solutions stemming from this analytical study are expected to be used for validation of non-adiabatic compressible flows simulation tools. A short validation of a numerical code using the [START_REF] Godounov | Résolution numérique des problemes multidimensionnels de la dynamique des gaz[END_REF] method was provided, yielding excellent results and very small relative errors.

To go further in this research work, one can expect that the numerical methods will be the favored method because of the complexity of the flow. Considering a non-adiabatic liquid flow, one can reasonably expect a mass transfer leading to a two-phase (liquid/vapor) flow. As we mentioned, some analytical solutions were also developed for two-phase flows in [START_REF] Fathalli | Liquid-vapor flows in heated ducts: Reference solutions[END_REF]: like the present paper, these studies provide help in analysing the flow as well allowing to check the quality of the codes by comparison with analytical solutions. Thus we have the following positive Mach numbers:

M sub x (P in ) = -b + √ ∆ 2a ≤ 1 M sup x (P in ) = -b - √ ∆ 2a ≥ 1
To be physical, there are two criteria to respect:

∆ ≥ 0 -b ± √ ∆ 2a ≥ 0
B Solving Q min cool from equation (30)

The limit of our computation is Π 1 = P out /P 0 = 1. Equation ( 25) is evaluated at the outlet for outlet pressure P out value:

Π 1 = P out P 0 = P * in + P ∞ P 0 1 + γM 2 in 1 + γM sub out 2 - P ∞ P 0 = 1
The isentropic relation (16) yields the following expression for inlet Mach number M in :

M 2 in = 2 γ -1 P 0 + P ∞ P * in + P ∞ γ-1 γ - 1
The lowest value of equation ( 24) for a choked flow is evaluated at the outlet for M sub

Figure 1 :

 1 Figure 1: Flow in a duct receiving or losing a heat flux

Figure 7 :

 7 Figure 7: Relationship between heat flux variables: φ(x) (left) and Q(x) (right)

Figure 8 :

 8 Figure8: Mach number (left) and pressure (right) evolution for a choked flow within a duct. The heated subsection (red line) starts subsonically and reaches the sonic point (M x ≤ 1, dotted line), but two possible paths exist within the cooled subsection (blue lines), depending on Π: subsonic (M x < 1, dotted line) or supersonic (M x > 1, dashed line).

2 )Figure 11

 211 Figure 10: ṁ * s evolution depending on heating power. When qs,heat increases, ṁ * s decreases.

Figure 16 :

 16 Figure16: Analytical and numerical results comparison of a gas being heated and cooled down, for Π > Π 1 , with mesh refinement.

  250 kW We compute two different fluids: a gas (P ∞ = 0) and a liquid (P ∞ ̸ = 0): • Gas: γ = 1.358, C v = 1247, e ref = 1.97 • 10 7 J/kg, P ∞ = 0 Pa • Liquid: γ = 3.423, C v = 1231.2, e ref = -1.15 • 10 6 J/kg, P ∞ = 10 4 Pa

Figure 17 :Figure 18 :

 1718 Figure 17: Analytical and numerical results comparison of a gas being heated and cooled down, for Π < Π 1 , with mesh refinement.

Figure 19 :

 19 Figure 19: Analytical and numerical results comparison of a liquid being heated and cooled down, for Π < Π 1 , with mesh refinement.

AppendicesA 2 A 2 A

 22 Solving M x from equation (24)f (M x ) = M 4 x βB -γ (P in ) = M 2 in 1 + βM 2 in B(P in ) = (1 + γM 2 in ) 2Solving the previous equation yields 4 roots, 2 negatives which are not physical and 2 positives:M x (P in ) =

2 A

 2 (P * in ) = M 2 in 1 + βM 2 in B(P * in ) = (1 + γM 2 in ) 2(31)

Table 2 :

 2 Pressure ratios definition for heated and cooled duct.

  Figure 13: Behavior of the flow variables depending on different pressure ratios, for a fixed couple {Q heat = 200 kW, Q cool = -250 kW}. As Π decreases, the flow gets progressively closer to choking. Once Π < Π 1 , the steady flow remains supersonic in the cooled subsection (red solution).Figure14: ṁs evolution depending on pressure ratio Π for a fixed couple {Q heat = 200 kW, Q cool = -250 kW}. As long as Π > Π 1 , when the pressure ratio decreases, the flow accelerates and the mass flow rate increases. When Π < Π 1 , ṁs ≡ ṁ * s no longer varies as the flow is choked.

		1.05 1.1 1.02 1 1.01							Q out = 0W Q out > 0W Q out < 0W		1.15 0.95 1.2 0.9	π 1 π 3					Q out = 0W Q out > 0W Q out < 0W
	Pressure (Bar) Pressure (Bar)	0.7 0.94 0.75 0.8 0.85 0.9 0.95 1 0.95 0.96 0.97 0.98 0.99								Mach number Pressure (Bar)	0.75 0.65 0.8 0.85 0.9 0.95 1 1.05 1.1 0.7 0.75 0.8 0.85
		0.6 0.65 0.93	0 0	π 1 π 3	0.2 0.2	0.4 0.4	x (m) x (m)	0.6 0.6	0.8 0.8	1 1	0.7 0.6 0.65	0 0	0.2 0.2	0.4 0.4	x (m) x (m)	0.6 0.6	0.8 0.8	1 1
		370 0.85 0.86							Q out < 0W Q out > 0W Q out = 0W		500 1.15 520 1.2	π 1 π 3
		360 0.84									480 1.1
	Temperature (K) Mach number	330 340 350 0.8 0.81 0.82 0.83								Velocity (m/s) Mach number	380 400 420 440 460 0.95 1 1.05
		320 0.79									360 0.9
		310 0.78	0 0	π 1 π 3	0.2 0.2	0.4 0.4	x (m)	0.6 0.6	0.8 0.8	1 1	340 0.85 320	0 0	0.2 0.2	0.4 0.4	x (m)	0.6 0.6	0.8 0.8	1 1
							x (m)								x (m)
					240											π 1
					238
					236
				. 2 m s (kg/s/m )	230 232 234
					228
					226
					224
					222 0.6		0.62	0.64	0.66		0.68		0.7		0.72
										π = P out /P 0

The function f T0,out T0,in = Π 1 -1 = 0 is solved by a numerical method. Knowing the value of T0,out T0,in , one can then evaluate Q min cool .


Naoufal El Bekri

Lucas Drumetz

Franck Vermet

TIME CHANGED NORMALIZING FLOWS FOR ACCURATE SDE MODELING

Keywords: Stochastic differential equations, deep generative models, normalizing flows, time series, convex neural network

The generative paradigm has become increasingly important in machine learning and deep learning models. Among popular generative models are normalizing flows, which enable exact likelihood estimation by transforming a base distribution through diffeomorphic transformations. Extending the normalizing flow framework to handle time-indexed flows provided dynamic normalizing flows, a powerful tool to model time series, stochastic processes, and neural stochastic differential equations (SDEs). In this work, we propose a novel variant of dynamic normalizing flows, a Time-Changed Normalizing Flow (TCNF), based on time deformation of a Brownian motion which constitutes a versatile and extensive family of Gaussian processes. This approach enables us to effectively model some SDEs that cannot be modeled otherwise, including standard ones such as the well-known Ornstein-Uhlenbeck process, generalizes prior methodologies, and leads to improved results and better inference and prediction capability.

INTRODUCTION

Dynamical systems have widespread use in various scientific areas such as finance, geosciences, and physics. The representation of these systems usually involves Ordinary Differential Equations (ODEs), or Stochastic Differential Equations (SDEs) [START_REF] Oksendal | Stochastic differential equations: an introduction with applications[END_REF] when noise and perturbations, on top of the deterministic component, are considered. Crucial applications include modeling volatility in financial data, or uncertainty quantification and propagation in geosciences. Tackling such systems through time series modeling and machine learning is an approach that experienced a surge in popularity, particularly recently, thanks to the generative paradigm, for forecasting applications, filtering, or interpolation with a notion of uncertainty in the generated sequences.

Popular generative models include Generative Adversarial Networks (GANs) [START_REF] Goodfellow | Generative adversarial nets[END_REF] and variational autoencoders (VAEs) [START_REF] Diederik | Auto-encoding variational bayes[END_REF], but also more recently normalizing flows (NFs) [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF][START_REF] Kobyzev | Normalizing flows: An introduction and review of current methods[END_REF] and diffusion/score based models [START_REF] Song | Generative modeling by estimating gradients of the data distribution[END_REF]. Though these models can be applied to generate time series, they are not well-suited for the task because they treat such data as vectors in R T , with T the number of time steps, without accounting for the causal structure. Adaptations of GANs, VAEs, and NFs to time series data have been carried out in [START_REF] Yoon | Time-series generative adversarial networks[END_REF][START_REF] Kidger | Neural sdes as infinite-dimensional gans[END_REF], [START_REF] Li | Scalable gradients and variational inference for stochastic differential equations[END_REF][START_REF] Zeng | Latent sdes on homogeneous spaces[END_REF], [START_REF] Mehrasa | Point process flows[END_REF][START_REF] Shchur | Intensity-free learning of temporal point processes[END_REF], respectively. In this work, we focus on NFs for their capacity to access explicit likelihoods, which is crucial for applications when uncertainty quantification or anomaly detection is required.

NFs are based on the celebrated change of variables formula that provides an expression of the probability density function of diffeomorphic transformations of a random variable. With carefully chosen transformations (or compositions thereof), if the initial density is tractable (explicit likelihood and easy sampling, in most cases Gaussian), the transformed density can be easily manipulated and sampled as well, provided that the Jacobian of the transformation can be computed efficiently. By considering the theoretical limit where an infinite number of transformations is applied, we can derive the Continuous Normalizing Flow (CNF) [START_REF] Grathwohl | FFJORD: freeform continuous dynamics for scalable reversible generative models[END_REF]. In this instance, the NF is described by an ODE that can be integrated to obtain the resulting density. This approach further improves the computational efficiency of this class of models by replacing the Jacobian determinant with the integration of its trace.

NFs have been extended to the dynamic setting by replacing the tractable base distribution with a tractable stochastic process, i.e. a Brownian motion [START_REF] Deng | Modeling continuous stochastic processes with dynamic normalizing flows[END_REF], making this type of model much more efficient for time series generation. However, it was noted in [START_REF] Deng | Continuous latent process flows[END_REF] that these models are theoretically unable to handle some of the most basic and common processes, such as the classical Ornstein-Uhlenbeck process.

Thus, in this paper, we propose a generalization of these approaches by using a large family of Gaussian processes as the underlying base process instead of the conventional Brownian motion. The Gaussian processes are constructed by transforming through time the standard Brownian motion, giving rise to Time Changed Normalizing Flow (TCNF), a model that possesses mathematical properties that enable it to describe dynamics and SDEs that cannot be captured by previous flow-based models, while maintaining the expressivity of dynamic NFs. We corroborate these findings by numerical experiments on several well-known processes.

The remainder of this paper is organized as follows: we first provide, in section 2, an overview of the neural SDEs, wherein both the drift and diffusion are neural networks, and the dynamic normalizing flow approach and a discussion on the inherent limitations of such models. Afterward, in section 3, we introduce our model and describe its characteristics and the training algorithm. Finally, quantitative results are presented in section 4 and compared to other flow-based models, and provide concluding remarks in section 5. Fig. 1: A time-change ϕ is applied to the Wiener process (red) to create a new Gaussian process, which is then mapped through a bijection f θ to the observed process (blue).

BACKGROUND

Neural Stochastic Differential Equations

We consider a filtered probability space (Ω, F, P ) and a time horizon T . A diffusion process X = {X t } t∈[0,T ] is defined by the Itô Stochastic Differential Equation (SDE):

dX t = µ(X t , t)dt + σ(X t , t)dW t , t ∈ [0, T ] (1)
where W = {W t } t∈[0,T ] is the m-dimensional adapted standard Wiener process (or Brownian motion). Functions µ :

R d × [0, T ] -→ R d and σ : R d × [0, T ] -→ R d×m
are the drift and diffusion coefficients, respectively. When µ and σ are implemented by neural networks, the SDE is designated as a neural SDE [START_REF] Tzen | Neural stochastic differential equations: Deep latent gaussian models in the diffusion limit[END_REF][START_REF] Liu | Neural sde: Stabilizing neural ode networks with stochastic noise[END_REF].

Several works have been proposed to learn neural SDEs using different generative modeling frameworks including Variational Autoencoders (VAEs) [START_REF] Li | Scalable gradients and variational inference for stochastic differential equations[END_REF][START_REF] Zeng | Latent sdes on homogeneous spaces[END_REF] and Generative Adversarial Networks (GANs) [START_REF] Kidger | Neural sdes as infinite-dimensional gans[END_REF]. In this paper, we specifically focus on the normalizing flow paradigm.

Normalizing flows

A normalizing flow [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF][START_REF] Kobyzev | Normalizing flows: An introduction and review of current methods[END_REF][START_REF] Grathwohl | FFJORD: freeform continuous dynamics for scalable reversible generative models[END_REF] is a transformation designed to model a random variable X and its corresponding complex distribution p X through a base distribution p Z and a differentiable bijective function f : R d -→ R d . Such modeling allows for both exact density estimation and efficient sampling, by using the change of variable formula for X = f (Z):

log p X (x) = log p Z (z) -log |det J f (z)| (2)
where the Jacobian J f (z) = ∂fi ∂zj 1≤i,j≤d is the d × d matrix of all partial derivatives of f . Previous works have extended this framework to model time series and stochastic processes by employing a bijection that is continuously indexed by time F (., t), along with a Brownian motion as the base process, giving rise to a Continuous Time Flow Process (CTFP) [START_REF] Deng | Modeling continuous stochastic processes with dynamic normalizing flows[END_REF]:

X t = F (W t , t).
Another approach proposed by [START_REF] Deng | Continuous latent process flows[END_REF] involves incorporating latent dynamics from an Ornstein-Uhlenbeck process combined with a normalizing flow to effectively model SDEs.

These models have demonstrated the effectiveness of dynamic NFs in capturing the complex behavior of various types of stochastic processes and SDEs. However, it is important to highlight that these models have inherent limitations. One constraint arises when applying Itô's lemma to CTFP to derive the univariate Ornstein-Uhlenbeck (OU) process described by the equation:

dY t = -a(Y t -b)dt + σdW t (3)
Indeed, by applying Itô's lemma to F (W t , t) we obtain:

dF (W t , t) = ∂F ∂t (W t , t)dt + ∂F ∂x (W t , t)dW t + 1 2 ∂ 2 F ∂x 2 (W t , t)dt (4)
By comparing both Eq. (3) and Eq. ( 4), we deduce that to model the OU process, we need to have ∂F ∂x (W t , t) = σ, implying that F (W t , t) = σW t + g(t), where g is a given differentiable function. However, differentiating this relation wrt t and plugging in (3) leads to the following condition:

dg dt (t) + ag(t) -ab = -aσW t (5)
Eq. ( 5) is not feasible as the left-hand side is a deterministic function of t while the right-hand side is a stochastic one depending on W t . Thus, CTFP shows limitations and falls short in its ability to model various stochastic processes effectively.

In the following section, we propose a model that can address this limitation and achieve improved results.

DYNAMIC NORMALIZING FLOW WITH TIME-CHANGE

Time-changed Normalizing Flow

We propose to model an observed stochastic process, denoted as X = {X t } t∈[0,T ] , by combining a normalizing flow and a time-changed Wiener process to accurately capture the dynamic behavior of X t based on a realized time series {(x ti , t i )} n i=1 . In this paper, we explicitly address the univariate case, with ongoing development of the general case which requires appropriate time-change for each dimension. We introduce the concept of a time-changed normalizing flow (TCNF), defined as follows:

X t = f θ W ϕ(t) , ϕ(t) , ∀t ∈ [0, T ],
where f θ (., t) : R -→ R is a differentiable bijection parameterized by θ, while W ϕ(t) denotes a Brownian motion with a time-change [18, Sect. 1 Chap. 0]. The time-change is given by ϕ : R + -→ R + , which is a measurable, positive and increasing function. The measurable and positive properties ensure the correct definition of W ϕ(t) , while the increasing property ensures the existence of its moments. Consequently, the neural network modeling the time-change must have intrinsically positive and increasing attributes. The property of time change has important applications as it produces a family of Gaussian processes that are more general than Brownian motion. The Dubins-Schwarz theorem [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF]Theorem 5.1.6] further emphasizes this property as it states that every local martingale is simply a time-changed Brownian motion.

Therefore, by making the base process of our model a time-changed Brownian motion, we can accurately capture all instances of local martingales and semimartingales, thus generalizing the CTFP setting. In fact, the solution to Eq. ( 3) can be expressed as follows:

Y t = Y 0 e -at + b(1 -e -at ) + σe -at √ 2a W e 2at -1 (6)
which can be properly modeled by TCNF. More general cases like processes with time-dependent volatility can also be expressed via a time-change and therefore modeled by TCNF. Finally, for ϕ(t) = t we recover the CTFP setting, which is suitable for modeling SDEs that do not require a time-change as the Geometric Brownian Motion (GBM) [START_REF] Oksendal | Stochastic differential equations: an introduction with applications[END_REF].

Time-change function

To handle the time-change function, we employ a convex neural network that ensures a positive gradient, thereby guaranteeing a monotone output. Specifically, we utilize the M-MGN architecture [START_REF] Chaudhari | Learning gradients of convex functions with monotone gradient networks[END_REF] based on K network modules defined as follows:

tk = W k × t + b k , M-MGN(t) = a + V ⊤ V t + K k=1 s k ( tk ) × W ⊤ k σ k ( tk ) (7)
where W k , b k ∈ R l×1 are respectively weight and bias vectors of the k th layer, σ k : R l -→ R l is an activation function and s k : R l -→ R its antiderivative. a ∈ R and V ∈ R l×1 are additional network parameters. As the result of M-MGN is not necessarily positive, we apply a translation of the output to ensure that the time-change is positive.

Training algorithm

The purpose is to train the TCNF in order to maximize the log-likelihood of the observed dataset {(x ti , t i )} n i=1 :

L = log p Xt 1 ,...,Xt n (x t1 , ..., x tn ) (8)
To compute Eq.( 8), we use the change of variable formula and leverage the independence of increments W ϕ(ti) -W ϕ(ti-1) . Thus, the log-likelihood is expressed as follows:

L = n i=1 log p W ϕ(t i ) |W ϕ(t i-1 ) w ϕ(ti) -log det ∂f θ w ϕ(ti) , ϕ(t i ) ∂W ϕ(ti) , (9)
where w ϕ(ti) = f -1 θ (x ti ; ϕ(t i )) and p W ϕ(t i ) |W ϕ(t i-1 ) denotes the conditional Gaussian distribution with mean W ϕ(ti-1) and variance ϕ(t i ) -ϕ(t i-1 ). This constitutes one notable difference from the log-likelihood of CTFP which uses a Gaussian distribution with the same mean but variance t i -t i-1 .

EXPERIMENTS

Toy datasets

To evaluate the performance of our proposed model, we conducted experiments on three toy datasets comprising univariate unitless time series. These datasets were generated by sampling from three different stochastic processes. Also, in our experiments we used an identical architecture as CTFP, leveraging CNFs.

The first dataset (Toy-SDE1) was generated by discretizing the OU process, given by the equation: dX t = -θ(X tµ)dt + σdW t , where µ and σ represent constant parameters for the drift and volatility terms, respectively. Parameter θ captures the speed at which a given sample path of the process converges towards the drift term. This dataset aims to assess the model's ability to capture the dynamics of time changes.

The second dataset (Toy-SDE2) was generated based on the equation:

dX t = -θ(X t -µ)dt + σ √ tdW t
, describing an OU with a time-dependent diffusion coefficient, and is used to test the model's capability to capture time transformations with increasing complexity. Notably, this SDE is of interest as it is commonly used in score-based models [START_REF] Yang | Diffusion models: A comprehensive survey of methods and applications[END_REF], where noise is gradually introduced during the training process.

Finally, The third dataset (Toy-SDE3) involved the geometric Brownian motion, described by the equation: dX t = µX t dt + σX t dW t , where µ and σ represent constant parameters for the drift and volatility terms, respectively. This dataset was designed to showcase the capacity of TCNF to handle SDEs where no time-change is required, effectively learning the simple function ϕ(t) = t. This shows that our approach can encompass the CTFP framework. The quantitative comparison is carried by comparing estimations of the mean m Xt , the standard-deviation σ Xt , the interquartile range IQR= Q 3 -Q 1 , and the density p Xt respectively. For each model, we compute the mean absolute errors (MAE) against the ground-truth values. The mean, standard deviation, and quartiles are estimated based on 1000 sample paths over 1000 iterations, while the density is estimated by the change of variable formula on a grid consisting of 1000 spatial points and 500 temporal points within the time interval [0, T = 1.5]. The quantitative results reported in Table 1 and Figure 2 show first that TCNF shows no loss in generality as it handles cases where no time-change is required, and second that our model exhibits superior estimation capability, as it can capture the behavior of time-changed ground truth solutions.

Real-world datasets

To further assess our model's capacity to capture complex dynamics, we train it on two real-world datasets: Crypto- forecasting (Crypto) [START_REF] Ticchi | G-research crypto forecasting[END_REF] and Electric Consumption Load (ECL) [START_REF] Zhou | Informer: Beyond efficient transformer for long sequence time-series forecasting[END_REF]. The Crypto dataset contains historical prices of various cryptocurrencies. We focused our analysis on modeling Ethereum log returns over the 2020 period. The ECL dataset comprises electricity consumption data from multiple clients within a 15-minute interval. We choose to model the consumption of client '200' for its extended time series.

The results comprise mean absolute errors (MAE) for estimating the mean (m Xt ) and standard deviation (σ Xt ) for the Crypto dataset. We employed mean relative errors (MRE) for the ECL dataset to scale the results appropriately. These results are reported in Table 2 and compared against CTFP.

CONCLUSION

We introduced a generalized approach to modeling SDEs via dynamic NFs and time-change. By transforming the Wiener process through time, we generate various Gaussian processes, which are then mapped to the observed process via a bijection. The time-change combined with the NF enabled us to model processes that are otherwise challenging to derive due to calculus constraints. Importantly, this extension retains the advantages of dynamic NFs, such as exact density estimation and efficient sampling.

Experiments showed that our model exhibits better performance and a generalization capacity. We believe incorporating dimension-specific time-changes enables us to extend the method to higher dimensions. Additionally, improvement in the calibration of the time-change can be achieved by linking it to either the quadratic variation of the process or its moments.

  This work was supported by Agence Nationale de la Recherche under grant ANR-21-CE48-0005 LEMONADE, and by France 2030 framework program, Centre Henri Lebesgue, under grant No ANR-11-LABX-0020-01.

Fig. 2 :

 2 Fig. 2: Comparison between TCNF and CTFP on toy-SDE1 (top) and toy-SDE2 (bottom). In each figure we depict density MAE between the flow models and the ground truth (GT), along the mean (continuous line) and the IQR (dashed line).

Table 1 :

 1 Quantitative analysis: We display mean absolute errors (MAE) for estimating respectively the mean, standard deviation, interquartile range IQR = Q 3 -Q 1 , and density. The estimations are based on 1000 sample paths over 1000 iterations.

	Dataset	Model	m Xt	σ Xt	IQR	p Xt
	Toy-SDE1	CTFP TCNF 0.019 ± 0.007 0.109 ± 0.009 0.700 ± 0.014 0.033 0.022 ± 0.009 0.241 ± 0.012 0.886 ± 0.021 0.062
	Toy-SDE2	CTFP TCNF 0.024 ± 0.012 0.322 ± 0.016 1.183 ± 0.024 0.061 0.059 ± 0.011 0.523 ± 0.018 1.517 ± 0.037 0.087
	Toy-SDE3	CTFP 0.035 ± 0.020 0.086 ± 0.056 TCNF 0.036 ± 0.018 0.099 ± 0.072	1.103 ± 0.037 1.083 ± 0.038 0.007 0.005

Table 2 :

 2 Quantitative analysis of real-world dataset: We display mean and standard deviation estimation errors.

	Dataset Model	m Xt	σ Xt
	Crypto	CTFP TCNF 0.025 ± 0.001 0.131 ± 0.005 0.083 ± 0.007 0.232 ± 0.009
	ECL	CTFP TCNF 0.299 ± 0.007 2.182 ± 2.081 0.771 ± 0.108 26.972 ± 8.634