- Email: [email protected]
Hyperbolicity of the heat equation
3rd IFAC IFAC Workshop Workshop on on Thermodynamic Thermodynamic Foundations Foundations for for a a 3rd Mathematical Systems 3rd IFAC Workshop on Theory Thermodynamic Foundations Available onlinefor at awww.sciencedirect.com Mathematical Systems Theory 3rd IFAC Workshop on Theory Thermodynamic Foundations for a Louvain-la-Neuve, Belgium, July Mathematical Systems Louvain-la-Neuve, Belgium, July 3-5, 3-5, 2019 2019 Mathematical Systems Theory Louvain-la-Neuve, Belgium, July 3-5, 2019 Louvain-la-Neuve, Belgium, July 3-5, 2019
ScienceDirect
IFAC PapersOnLine 52-7 (2019) 63–67
Hyperbolicity Hyperbolicity Hyperbolicity Hyperbolicity
of of of of
the the the the
heat heat heat heat
equation equation equation equation∗∗
Guilherme Guilherme Ozorio Ozorio Cassol Cassol and and Stevan Stevan Dubljevic Dubljevic ∗ Guilherme Ozorio Cassol and Stevan Dubljevic ∗∗ Guilherme Ozorio Cassol and Stevan Dubljevic ∗ ∗ Department of Chemical and Materials Engineering ∗ Department of Chemical and Materials Engineering of Chemical and Materials Engineering University of Edmonton ∗ Department University of Alberta, Alberta, Edmonton Department of2V4 Chemical [email protected]) Materials Engineering University of Alberta, Edmonton Alberta, Canada T6G (e-mail: Alberta, Canada University T6G 2V4 (e-mail: [email protected]) of Alberta, Edmonton Alberta, Canada T6G 2V4 (e-mail: [email protected]) Alberta, Canada T6G 2V4 (e-mail: [email protected]) Abstract: In In this this manuscript, manuscript, aa comparison comparison between between the the parabolic parabolic and and hyperbolic hyperbolic partial partial Abstract: Abstract: In this manuscript, a comparison between the parabolic and hyperbolic partial differential equations for heat diffusion is studied. First, numerical results and an eigenvalue differential equations for heat diffusion is studied. First,the numerical results an eigenvalue Abstract: In thistwo manuscript, a comparison between parabolic and and hyperbolic partial differential for heatofdiffusion is are studied. First, numerical results and an eigenvalue analysis for for equations these types equations shown, which help to understand understand the difference analysis these twofor types ofdiffusion equations are shown, which help to the difference differential equations heat is studied. First, numerical results and an eigenvalue analysis for these two types of equations are shown, which help to understand the difference in the system dynamics. Then, both equations are also considered in a Stefan problem for the in the system dynamics. Then, equations are alsowhich considered a Stefan problem for the analysis foranthese two in types of both equations are shown, help toin understand the difference in the system dynamics. Then, both equations are also considered in a Stefan problem for the melting of ice block a one-dimensional setting. The results show that the hyperbolic partial melting of an ice block in aThen, one-dimensional setting. The results showinthat the hyperbolic partial in the system dynamics. bothspeed equations are also considered a Stefan problem for the melting of an ice blockshows in a one-dimensional The results show the hyperbolic partial differential equation shows a finite finite ofsetting. propagation of heat heat andthat can represent the system system differential equation a speed of propagation of and can represent the melting of an ice block in a one-dimensional setting. The results show that the hyperbolic partial differential equation shows a finite speed of propagation of heat and can represent the system as properly properly as as the the parabolic parabolic equation. equation. as differential shows aequation. finite speed of propagation of heat and can represent the system as properly equation as the parabolic as properly the parabolic equation. © 2019, IFACas(International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Distributed-parameter Distributed-parameter systems; systems; Moving Moving boundary boundary conditions; conditions; Heat Heat flows; flows; Linear Linear Keywords: Keywords: Distributed-parameter systems; Moving boundary conditions; Heat flows; Linear systems; Eigenvalue Eigenvalue problems; systems; problems; Keywords: Distributed-parameter systems; Moving boundary conditions; Heat flows; Linear systems; Eigenvalue problems; systems; Eigenvalue problems; 1. INTRODUCTION INTRODUCTION This manuscript manuscript considers considers both both hyperbolic hyperbolic and and parabolic parabolic 1. This 1. INTRODUCTION This manuscriptThe considers bothand hyperbolic and heat equations. equations. The derivation and analysis of of theparabolic numeriheat derivation analysis the numeri1. INTRODUCTION This manuscript considers hyperbolic and heat equations. The derivation and analysis theparabolic numerical results results for both both cases areboth considered in the theof next section. cal for cases are considered in next section. heat equations. Thecases derivation and analysis ofnext the numeriTransport-reaction distributed parameter system (DPS) cal results for both are considered in the section. After that, the analysis of the eigenvalues associated with Transport-reaction distributed parameter system (DPS) cal After that,for theboth analysis of the eigenvalues associated with results cases of are considered in the nextare section. Transport-reaction distributed parameter manufacturing system (DPS) After models present in chemical, petrochemical, that, the analysis the eigenvalues associated with each equation is shown. Finally, both equations conmodels present in chemical, petrochemical, manufacturing each equation is shown.ofFinally, both equations are conTransport-reaction distributed parameter manufacturing system (DPS) that, the eigenvalues associated models present in chemical, petrochemical, and process process industry take the the mathematical form given given by After each equation isanalysis shown. Finally, both equations are with considered in aathe moving boundary problem, which numerical sidered in moving boundary problem, which numerical and industry take mathematical form by models present in chemical, petrochemical, manufacturing equation is shown. Finally,problem, both equations are conand process industry take the mathematical formThe given by each hyperbolic partial differential equations (PDEs). The main sidered in ashown moving boundary which numerical results are shown in the the subsequent section. hyperbolic partial differential equations (PDEs). main results are in subsequent section. and process industry take the mathematical form given by sidered in a moving boundary problem, which numerical hyperbolic partial equations The main results are shown in the subsequent section. conservation laws differential are embedded embedded in aa (PDEs). modelling variety conservation laws are in modelling variety hyperbolic partial differential equations (PDEs). The main results are shown in the subsequent section. conservation lawshyperbolicity are embedded modelling variety provided by by the hyperbolicity of in thea transport transport systems 2. THE THE HEAT HEAT EQUATION EQUATION provided the of the systems conservation lawshyperbolicity are embedded in a transport modellingassystems variety 2. provided by the of the which is physically relevant and desired property action 2. THE HEAT EQUATION which is physically relevant and desired property as action provided byisthe hyperbolicity of the transport THE HEAT EQUATION which is physically relevant and desired property assystems action at distance distance precluded and physically meaningful finite Considering a 2. at is precluded and physically meaningful finite fluid at at rest rest with with constant constant density density and and nenewhich is physically relevant and desired property as action Considering a fluid at distance is precluded and physically meaningful finite speed of phenomena propagation is ensured. On the other Considering a fluid at rest with constant density and nespeed of phenomena propagation is ensured. On the other glecting non-linear terms in gradients and time-derivatives, at distance is precluded and physically meaningful finite glecting non-linear in gradients and time-derivatives, speed of propagation is ensured. On the hand, thephenomena hyperbolicity mathematically ensures the other well- Considering a fluidterms at with constant density and nehand, the hyperbolicity mathematically ensures the wellglecting non-linear terms in gradients and time-derivatives, the energy energy balance forrest the system leads to the following following speed of phenomena propagation is ensured. On the other the balance for the system leads to the hand, the hyperbolicity mathematically ensures the wellposedness of local Cauchy problems (Fischer and Marsden glecting non-linear terms in gradients and time-derivatives, the energy balance for the system leads to the following posedness of local Cauchy problems (Fischer and Marsden parabolic equation: hand, the of hyperbolicity ensures well- the parabolic equation: posedness local Cauchymathematically problems (Fischer and the Marsden (1972)). energyequation: balance for the system leads to the following parabolic (1972)). posedness of local Cauchy problems (Fischer and Marsden (1972)). parabolic equation: ∂T Diffusive transport transport of of heat heat across across macroscopic macroscopic length length ∂T (1972)). = α∆T α∆T (1) Diffusive = (1) ∂T ∂t = Diffusive transport of heat across macroscopic length scales is well described by Fourier’s law. However, any ∂t α∆T (1) ∂T scales is well described by Fourier’s law. However, any Diffusive transport heat across body macroscopic length ∂t = α∆T (1) scales well described by material Fourier’s law. is However, any initial is disturbance inof the the material is propagated initial disturbance in body propagated The heat heat flux flux can be given ∂t scales is well described by Fourier’s law. However, any The can be given by by the the Fourier’s Fourier’s law law and and initial disturbance in the material body is propagated instantly due due to to the the parabolic parabolic nature nature of of the the partial partial differdifferinstantly The heat flux given by Fourier’s considering onlycan one be dimension it is isthe defined as: law and initial disturbance the material propagated considering only one dimension it defined as: instantly due toobtained the in parabolic nature ofJordan theispartial differential equation equation obtained (Christov andbody Jordan (2005)). To The heat flux given by the Fourier’s law and ential (Christov and (2005)). To considering onlycan one be dimension it is defined as: instantly due to the parabolic nature of the partial differT is defined as: (2) q= = −k∂ −k∂ζitT ential equation obtained (Christov Jordanmodified (2005)).the To considering only one dimension eliminate this unphysical unphysical feature, and Cattaneo modified the (2) q ζ eliminate this feature, Cattaneo ential equation obtained (Christov and Jordan (2005)). To q = −k∂ζ T (2) eliminate this unphysical feature, Cattaneo modified the Fourier law to take into account the thermal inertia, which Fourier lawthis to take into account theCattaneo thermal inertia, which = −k∂ζ conductivity. T (2) eliminate unphysical feature, modified the where where k k is is the the material materialqthermal thermal The energy energy Fourier law to take into account the thermal inertia, which avoids the phenomenon of infinite propagation (Cattaneo conductivity. The avoids the phenomenon of infinite propagation (Cattaneo Fourier law to take into account the thermal inertia, which where k is the material thermal conductivity. The energy balance can be rewritten as: avoids phenomenon infinite (Cattaneo balance (1958);the Vernotte (1958); of Jou et al. al. propagation (2001)). can be material rewrittenthermal as: is the conductivity. The energy (1958); Vernotte (1958); Jou et (2001)). avoids phenomenon infinite (Cattaneo where balancek can be rewritten as: (1958);the Vernotte (1958); of Jou et al. propagation (2001)). ∂ T + ∂ q = 00 → ∂ T (3) t ζ tT balance can be rewritten as: ∂ T + ∂ q = → ∂ T = = α∂ α∂ζζ (3) t ζ t ζζ T Moreover, a Stefan problem is also considered, which is (1958); Vernotte (1958); Jou et Moreover, a Stefan problem is al. also(2001)). considered, which is ∂t T + ∂ζ q = 0 → ∂t T = α∂ζζ T (3) Moreover, a Stefan problem is also considered, which is a specific type of boundary value problem for a partial ∂ T + ∂ q = 0 → ∂ T = α∂ T (3) k k t= k , ζρ is the material t ζζ a specific type of boundary value problem for a partial k with α = density, c is the Moreover, aequation, Stefan problem is focusing alsoproblem considered, which is with α = ρc adifferential specific type of boundary value for a partial differential generally in the heat disk = C k , ρ is the material density, c is the ρc C equation, generally focusing in the heat dis- material with α =heat = Ck , ρ is material density, is the capacity pertheunit unit of mass mass and C Cc is is the adifferential specific in type of boundary value problem for a moving partial k capacity ρc equation, generally focusing in the dis- material per of and the tribution, phase changing medium. Since theheat with α =heat = C , ρ is the material density, c ρc capacity tribution, inequation, aa phase changing medium. Since the moving material heat capacity per unit of mass and C is the heat per unit of volume. For an initial differential generally focusing in the heat dismaterial heat capacity per unit of volume. For an initial tribution, in a phase changing medium. Since the moving interface is unknown a priori, the solution needs to take heat capacity capacity per=unit unit ofvolume. mass and C is the interface isinunknown a priori, medium. the solution needs to take material material heat per of For an initial value problem (T (ζ, t = 0) T (ζ)) in an one dimensional tribution, a phase changing Since the moving 0 problem (T (ζ, t = per 0) =unit T0 (ζ)) in an oneFor dimensional interface is unknown a priori, theof needs to take value into account account for the the determination determination ofsolution the moving moving boundary material heat capacity of volume. an initial into for the boundary value problem (T (ζ, t = 0) = in an dimensional space, solution of this PDE is given as: interface is unknown a of priori, to take infinite 0 (ζ)) into account forexample the determination ofsolution the moving position. One example thesethe problems isneeds theboundary diffusion infinite space, the the solution of T PDE is one given as: value problem (T (ζ, t =∞ 0) = Tthis in one dimensional an 0 (ζ)) position. One of these problems is the diffusion infinite space, the solution of this PDE is given 2as: into account for the determination of the moving boundary ∞ 2 1 (z − ζ) position. One example of these problems is the diffusion of heat in the melting of ice: as the melting occurs the infinite space, the solution of this PDE is given as: 1 (z − ζ) of heat inOne the example melting of of these ice: asproblems the melting occurs the ∞ T (z, 0)exp − T (ζ, t) = 2 position. is the diffusion T (ζ, t) = (4παt) dz (4) (4) 1 3/2 − ζ)2 dz −∞ of heat in (interface the melting of ice:the as solid the melting occurs the boundary (interface between the solid and liquid liquid phase) ∞ T (z, 0)exp − (z 4αt 3/2 boundary between and phase) 4αt T (ζ, t) = (4παt) dz (4) T (z, 0)exp − 1 (z − ζ) of heat in the melting of ice: as the melting occurs the −∞ 3/2 boundary (interface between the solid and liquid phase) will be changing position. Although the natural occurrence 4αt (4παt) T (ζ, t) = dz (4) T (z, 0)exp − will be changing position. Although the natural occurrence 3/2 −∞ boundary (interface between the solid and liquid phase) As long as (4παt) 4αtfrom 0, this will be changing position. Although the natural occurrence −∞ the initial condition is different of a Stefan problem is mostly associated with the melting of a be Stefan problem is mostly associated with the melting As long as the initial condition is different from 0, this will changing position. Although the occurrence long as the initial condition different fromwhich 0, this of problem is mostly associated with the melting As solution predicts instant heatispropagation, propagation, which is anda Stefan solidification problems, there arenatural some Stefan-like solution predicts aa instant heat is and solidification problems, there are some Stefan-like As longthe as the initial condition ispropagation, different fromwhich 0, this of a Stefan problem is mostly associated with the melting solution predicts a instant heat is and solidification problems, there are some Stefan-like called Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous called the Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous solution predicts a instant heat propagation, which is and solidification problems, there are some Stefan-like called the Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous wrote a paper in which he addressed the question of the media or even shock waves in gas dynamics (Rubenstein wrote athe paper in Conduction which he addressed theCattaneo question (1958) of the media or even shock waves in gas dynamics (Rubenstein called Heat Paradox. problems related, for instance, to the fluid flow in porous wrote a paper which he addressed theFourier’s questionlaw of and the media of He (1971)).or even shock waves in gas dynamics (Rubenstein paradox paradox of heat heatinconduction. conduction. He modified modified (1971)). a paper which he addressed theFourier’s questionlaw of and the media paradox of heatinconduction. He modified Fourier’s law and (1971)).or even shock waves in gas dynamics (Rubenstein wrote paradox of heat conduction. He modified Fourier’s law and (1971)). 2405-8963 © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.
Copyright © 2019 85 Copyright 2019 IFAC IFAC 85 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 2019 IFAC 85 10.1016/j.ifacol.2019.07.011 Copyright © 2019 IFAC 85
IFAC TFMST 2019 64 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
The parabolic equation predicts an instant heat propagation (the initial slopes show that for t > 0, T (ζ, t) > 0 ∀ ζ), while the initial response for the 2nd order hyperbolic consider that there is some time-lag for the initial propagation (there is a finite propagation velocity in the system). These results can be interesting for chemical engineering process (for instance in tubular reactors Dochain (2016); Fogler (2016)), as generally the models of reaction-diffusion nature are represented as parabolic equations. Therefore, it would be interesting to study these models represented as a hyperbolic equation instead, as it would not characterize a instantaneous transport shown in parabolic partial differential equations. 3. EIGENVALUE ANALYSIS In this section, a comparison between the eigenvalues of the parabolic and second order hyperbolic is made to understand how the system dynamics are different for the two scenarios. The eigenvalue problem is defined as Aφ = λφ, where A is the operator considered. For the ∂2 parabolic scenario, A = ∂ζ 2 , with the boundary conditions considered ∂ζ x(ζ = 0) = ∂ζ x(ζ = 1) = 0. The solution for the eigenvalue problem for the parabolic equations is well known and given by: λn = −αn2 π 2 , n = 1, ..., ∞ (8) This result indicates that the eigenvalues grow unbounded. The partial differential equation shown in Eq. (7) can be written as: 0 1 T T ∂ = (9) α ∂ 2 −1 ∂T ∂t ∂T 2 τ ∂ζ τ ∂t ∂t And an operator A can be defined as:
Fig. 1. Comparison between the results from the Heat Equations: (Upper) Parabolic; (Bottom) 2nd order Hyperbolic
based his argument on the elementary kinetic theory of gases and argued that there is a time-lag between the start of the particles at their point of departure and the time of passage through the middle layer. Therefore, if the temperature changes in time, it is the heat flux at a certain time depends on the temperature gradient at an earlier time. This assumption leads to the following definition of a modified heat flux: (5) q = −k(1 − τ ∂t )∂ζ T If τ is small:
(1 − τ ∂t )−1 ≈ (1 + τ ∂t )
T ∂x A = α ∂ 2 −1 → = Ax, x = ∂T ∂t τ ∂ζ 2 τ ∂t 0
1
(10)
The eigenvalue problem is defined as Aφ = λφ, which gives the following system of equations: φ2 = λφ1 α ∂ 2 φ1 λφ1 φ2 α ∂ 2 φ1 = λ2 φ1 → − 2 = λφ − τ ∂ζ τ 2 τ ∂ζ 2 τ
(11)
Solving this second order ordinary differential equation with the boundary conditions ∂ζ φ1 (ζ = 0) = ∂ζ φ1 (ζ = 1) = 0 results in the following condition: τ λ2n + λn + αn2 π 2 = 0, n = 1, ..., ∞ (12) and for each value of n: −1 ± 1 − 4τ (αn2 π 2 ) λn = (13) 2τ
(6)
And the energy balance can be rewritten as a second order hyperbolic PDE: (7) ∂t T + ∂ζ q = 0 → τ ∂tt T + ∂t T = α∂ζζ T If τ = 0, the original parabolic PDE is obtained. Considering the heat conduction equations shown in (3) and (7) with ∂ζ T (ζ = 0) = ∂ζ T (ζ = 1) = 0, the results shown in Figure 1 are obtained for the same initial condition 2 (T0 (ζ) = 1 − (ζ−0.5) 0.252 , for 0.25 ≤ ζ ≤ 0.75; 0, otherwise) and time interval.
And this results in a very different set of eigenvalues when compared to the parabolic equation. This system does not have a infinite set of real fast eigenvalues as shown before (for the parabolic equation, λn → −∞ as n → ∞), in fact, the whole set of eigenvalues needs to be in the region 86
IFAC TFMST 2019 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
Liquid TI (0 ≤ ζ ≤s, t)
1 where −1 τ < Re(λn ) < 0. Furthermore, when n > 4τ απ , the eigenvalues start to present imaginary parts, which is different from the eigenvalues obtained for the parabolic equation that have just real parts. 4. HEAT EQUATION EQUATION AND THE STEFAN PROBLEM To analyze the difference between the results obtained by a hyperbolic and a parabolic partial differential equation, a heat diffusion problem with phase change is considered. One of the simplest mathematical model of the change of phase phenomenon is called Stefan problem. If a phase change of a material occurs at a given point latent heat is released or absorbed and the temperature of the material at that point remains constant. It is assumed that the temperature of phase change is constant and that the temperature across the material is continuous. For simplicity, an one dimensional case is considered in a material with constant cross sectional area. If the phase change happens at a point ζ = s(t), then: (14) TI (ζ = s(t) − δ, t) = TII (ζ = s(t) + δ, t) = TC where δ is a small spatial variation such that δ → 0, TC is the phase change temperature and the subscripts I and II are used to represent two different phases. Therefore, this condition implies that the temperature close to the phase transition point is continuous and is the same in both phases. The Stefan condition is obtained by applying a energy balance between two instants t0 and t1 in the volume of a material that is going through phase change. If the phase change point moved from s(t0 ) to s(t1 ) and that s(t1 ) > s(t0 ), the energy involved in the phase change is given by the volume of the material that went through phase change between this two instants and the latent heat required for the transition: (15) Q = Cl A[s(t1 ) − s(t0 )] where Q is the energy involved in the phase transition, Cl is the latent heat (per unit of volume) and A is the constant cross sectional area. This amount of energy should be equivalent to the heat provided by heat diffusion in both sides of the phase change interface in this time interval: t1 t1 (qI − qII )dAdt = A (qI − qII )dt (16) Q= t0
A
t0
Cl [s(t1 ) − s(t0 )] =
t0
TII (s ≤ ζ ≤1, t) = 0
ζ =1 T =0
Fig. 2. Representation of the system with phase transition ature of T0 (ζ) = 0 is considered throughout the system. As the phase change temperature considered is the same as the temperature at the boundary ζ = 1 and the initial temperature in the solid phase, there is no heat flux that flows through this phase (as ∂ζ TII (s≤ζ ≤1, t) = 0 and no heat sinks or sources are considered in 0<ζ <1). This simplifies the model to a so called one phase problem, as only one phase (in this case the liquid phase) has a temperature that is time and space dependent that needs to be found (TI (0≤ζ ≤s, t)). Therefore, the Stefan condition for the moving boundary for the original Fourier’s law can be simplified to: dt s = qI = −k ∂ζ T |ζ = s(t) → dt s = −β ∂ζ T |ζ=s(t)
(19)
where, for simplicity of the notation, T is now used to represent TI and β is defined as β = Ckl , with k as the material thermal conductivity and Cl the latent heat per unit of volume, defined previously. If the modified heat flux shown in Eq. (5) is used, the Stefan condition obtained is: dt s = qI = −k(1 − τ ∂t ) ∂ζ T |ζ=s(t) → τ dtt s + dt s = −β ∂ζ T |ζ=s(t)
(20)
This results in a second order ordinary differential equation for the modified heat flux, while the original Fourier’s law results in a first order ordinary differential equation. These equations need to be solved along with the corresponding hyperbolic and parabolic diffusion equations (Eq. (7) and Eq. (3), respectively).
t0
t1
Solid
ζ =0 ζ =s T =1 T =0
where qI and qII represent the heat flux from both phases and the fact that the cross sectional area does not change was used to simplify the integral. Then, t1 (qI − qII )dt → Q = Cl A[s(t1 ) − s(t0 )] = A
65
As these partial differential equations have a moving boundary, the following change of variables is used to change from a time-varying domain to a fixed domain: ζ → ∈ [0, 1] (21) = s(t)
(17)
(qI − qII )dt
Dividing the equation by t1 −t0 and setting t1 → t0 , finally gives the Stefan condition: (18) Cl dt s(t) = qI (t, ζ = s − δ) − qII (t, ζ = s + δ)
And the derivatives can be rewritten for this new spatial variable as: 1 ∂ ∂T ( , t) = ∂ T ( , t) (22) ∂ζ T (ζ, t) = ∂ζ ∂ s
The system considered here is shown in Figure 2, which represents melting ice with fixed temperatures at ζ = 0 and ζ = 1 ( T (ζ = 0) = 1 and T (ζ = 1) = 0).
∂ ∂ζζ T (ζ, t) = ∂
The transition temperature (melting temperature) of the solid phase is considered to be TC = 0. A initial temper87
∂ ∂T ( , t) ∂ζ ∂
1 ∂ = 2 ∂ T ( , t) ∂ζ s
(23)
IFAC TFMST 2019 66 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
∂t T (ζ, t) = ∂t T (, t) +
∂ ∂T (, t) = ∂t ∂
∂t T (, t) − dt s∂ T (, t) s
(24) 1 0.8
2 ∂tt T (, t) − dt s∂ ∂t T (, t)+ s 2 ∂ T (, t) 2 (dt s)2 − dtt s s s
Temperature
∂tt T (ζ, t) = ∂ ∂T (, t) ∂ ∂T (, t) ∂ 2 ∂tt T (, t) + 2 + = ∂ ∂t ∂t ∂ ∂t2
0.6 0.4 0.2 0 0.6 0 0.4
0.2 0.4 0.2
0.6 0.8
Time
(25)
Finally, the parabolic partial differential equation can be written as: 1 (26) α 2 ∂ T (, t) = ∂t T (, t) − dt s∂ T (, t) s s And for the hyperbolic: 1 α 2 ∂ T (, t) = s 2 τ ∂tt T (, t) − dt s∂ ∂t T (, t)+ s (27) 2 2 ∂ T (, t) 2 (dt s) − dtt s + s s ∂t T (, t) − dt s∂ T (, t) s
0
1
Fig. 3. Spatial profile obtained for the Stefan problem using a parabolic heat equation (s(t = 0) = 0.5)
1
Temperature
0.8
0.6
0.4
0.2
0 0.6 0.4
The Stefan conditions also need to be properly modified. The first order is changed to: 1 dt s = −β ∂ T (, t) (28) s While the second order is given by: 1 τ dtt s + dt s = −β ∂ T (, t) (29) s The boundary conditions also need to be redefined: u(ζ = 0, t) = 1 → u( = 0, t) = 1 (30) u(ζ = s(t), t) = 0 → u( = 1, t) = 0
0.2
Time
0
1
0.8
0.6
0.4
0.2
0
Fig. 4. Spatial profile obtained for the Stefan problem using a hyperbolic heat equation interface position at s(t = 0) = 0.5 was used to show the difference on the initial diffusion between the two types of partial differential equation. The results for α = β = 1 and τ = 10−3 are shown in Figure 3 and Figure 4 for the parabolic and hyperbolic diffusion equation, respectively. Both results show that initially there is a time for which the boundary does not move. This happens because the initial condition in the liquid is the same as the transition temperature and it takes some time for the heat to flow through the liquid phase. As the diffusion starts to happen and the heat reaches the interface between the two phases, the boundary starts to move (the solid phase starts to melt) and this happens until s(t) = 1. When this condition is reached, the boundary stops to move and the system is represented just by diffusion in a fixed-domain. After some time, the system reaches the expected steady-state.
For the second order hyperbolic equation, boundary conditions are also necessary for the first order time derivatives, which can be obtained directly from the boundary conditions: u( = 0, t) = 1 → ∂t u( = 0, t) = 0 (31) u( = 1, t) = 0 → ∂t u( = 1, t) = 0 To guarantee that s(t) ≤ 1, a condition is imposed to stop the growth of the moving boundary. If s(t) > 1, the first and second derivative of s(t) are set to zero. This makes Eq. (26) and Eq. (27) convert back to the original diffusion problem (Eq. (3) and Eq. (7), respectively), as expected if there is no more change in the domain.
The main difference between the results for the parabolic and the second order hyperbolic can be seen at the initial diffusion. As shown before, in the analysis of the results of the heat equation (Figure 1), the hyperbolic equation shows a finite velocity of propagation, while the parabolic equation represents a situation of instant diffusion. These results can be seen again in Figure 3 and Figure 4. For the parabolic equation, if t > 0 all points between 0 ≤ ζ ≤ s(t) present a temperature different from zero. For the
5. RESULTS This section presents the results for the hyperbolic and parabolic partial differential equations (and the respective second order and first order ordinary differential equations) with the conditions previously chosen. A initial 88
IFAC TFMST 2019 Louvain-la-Neuve, Belgium, July 3-5, 2019 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67
1
1
0.9
0.9 0.8 Par( =1, =1)
0.7
-3
Cat( =1, =1, =10 )
s (Moving boundary)
s (Moving boundary)
0.8
67
0.6 0.5 0.4 0.3
0.7 0.6 0.5
Par( =1, =1)
0.4
Cat( =1, =1, =10 -3 ) Par( =1, =10)
0.3
Cat( =1, =10, =10 -3 ) Par( =0.1, =10)
0.2
Cat( =0.1, =10, =10 -3 ) Par( =0.1, =1)
0.2 0.1
Cat( =0.1, =1, =10 -3 )
0.1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0
Time
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time
Fig. 5. Moving boundary dynamics for the Stefan problem using a parabolic and hyperbolic heat equation (s(t = 0) = 0.5)
Fig. 6. Moving boundary dynamics for the Stefan problem using a parabolic and hyperbolic heat equation for different values of α and β (s(t = 0) = 0)
hyperbolic, the initial slope after t > 0 shows the finite speed of propagation. Even though the initial diffusion is different, both hyperbolic and parabolic present a similar growth of the boundary, as it is shown in Figure 5 for s(t = 0) = 0.5. The difference between the partial differential equations is barely noticeable and the hyperbolic equation seems to be also adequate to model the system.
ACKNOWLEDGEMENTS The support for this work is provided by CAPES 88881.128514/2016-01 (Brazil) and support for Guilherme Ozorio Cassol is gratefully acknowledged. REFERENCES
Figure 6 shows the results for s(t = 0) = 0 using different values of α and β. As the value of α increases, faster is the growth of the moving boundary, as the heat gets to the interface of the two phases in a faster pace. This can be seen if the red line is compared to the black line, and the blue line with the green line. A faster growth of the boundary is also seem if β. The higher the value of this parameter, the smaller is the heat necessary for phase transition, which makes the boundary moves faster. This can be seen if the red line is compared to the green line and the black line with the blue line. Also, a decrease in one of these parameters is not compensated if the other is increased, as it can be conclude if the green and blue lines are compared. For all cases it is possible to see that the difference between the hyperbolic and parabolic equations is barely noticeable. Once again, the hyperbolic is adequate to represent the system.
Cattaneo, C. (1958). On a form of heat equation which eliminates the paradox of instantaneous propagation. C. R. Acad. Sci. Paris, 431–433. Christov, C.I. and Jordan, P.M. (2005). Heat conduction paradox involving second-sound propagation in moving media. Physical review letters, 94 15, 154–301. Dochain, D. (2016). Analysis of the multiplicity of equilibrium profiles in tubular reactor models. IFACPapersOnLine, 49(18), 903–908. Fischer, A.E. and Marsden, J.E. (1972). The Einstein evolution equations as a first-order quasi-linear symmetric hyperbolic system, I. Communications in Mathematical Physics, 28(1), 1–38. Fogler, H.S. (2016). Elements of chemical reaction engineering. Prentice Hall, 5 edition edition. Jou, D., Casas-Vazquez, J., and Lebon, G. (2001). Extended irreversible thermodynamics (3rd Edition). Springer. Rubenstein, L. (1971). The Stefan problem, by L.I. Rubenstein. Translations of mathematical monographs. American Mathematical Society. Vernotte, P. (1958). Les paradoxes de la theorie continue de l’equation de la chaleur. C. R. Acad. Sci. Paris, 246, 3154–3155.
6. CONCLUSION This manuscript compared the results obtained for the heat diffusion considering hyperbolic and parabolic equations. First, the simple case of one-dimensional heat diffusion was considered and the eigenvalue analysis was performed to access the difference between these two equations, which shows the finite speed of propagation of the hyperbolic equation. Finally, the Stefan problem was considered and the results show that the interface dynamic for the hyperbolic equation are similar to the parabolic in the cases analyzed, but the former finite speed of propagation would be more desirable to represent the dynamics of the actual physical system. 89
ScienceDirect
IFAC PapersOnLine 52-7 (2019) 63–67
Hyperbolicity Hyperbolicity Hyperbolicity Hyperbolicity
of of of of
the the the the
heat heat heat heat
equation equation equation equation∗∗
Guilherme Guilherme Ozorio Ozorio Cassol Cassol and and Stevan Stevan Dubljevic Dubljevic ∗ Guilherme Ozorio Cassol and Stevan Dubljevic ∗∗ Guilherme Ozorio Cassol and Stevan Dubljevic ∗ ∗ Department of Chemical and Materials Engineering ∗ Department of Chemical and Materials Engineering of Chemical and Materials Engineering University of Edmonton ∗ Department University of Alberta, Alberta, Edmonton Department of2V4 Chemical [email protected]) Materials Engineering University of Alberta, Edmonton Alberta, Canada T6G (e-mail: Alberta, Canada University T6G 2V4 (e-mail: [email protected]) of Alberta, Edmonton Alberta, Canada T6G 2V4 (e-mail: [email protected]) Alberta, Canada T6G 2V4 (e-mail: [email protected]) Abstract: In In this this manuscript, manuscript, aa comparison comparison between between the the parabolic parabolic and and hyperbolic hyperbolic partial partial Abstract: Abstract: In this manuscript, a comparison between the parabolic and hyperbolic partial differential equations for heat diffusion is studied. First, numerical results and an eigenvalue differential equations for heat diffusion is studied. First,the numerical results an eigenvalue Abstract: In thistwo manuscript, a comparison between parabolic and and hyperbolic partial differential for heatofdiffusion is are studied. First, numerical results and an eigenvalue analysis for for equations these types equations shown, which help to understand understand the difference analysis these twofor types ofdiffusion equations are shown, which help to the difference differential equations heat is studied. First, numerical results and an eigenvalue analysis for these two types of equations are shown, which help to understand the difference in the system dynamics. Then, both equations are also considered in a Stefan problem for the in the system dynamics. Then, equations are alsowhich considered a Stefan problem for the analysis foranthese two in types of both equations are shown, help toin understand the difference in the system dynamics. Then, both equations are also considered in a Stefan problem for the melting of ice block a one-dimensional setting. The results show that the hyperbolic partial melting of an ice block in aThen, one-dimensional setting. The results showinthat the hyperbolic partial in the system dynamics. bothspeed equations are also considered a Stefan problem for the melting of an ice blockshows in a one-dimensional The results show the hyperbolic partial differential equation shows a finite finite ofsetting. propagation of heat heat andthat can represent the system system differential equation a speed of propagation of and can represent the melting of an ice block in a one-dimensional setting. The results show that the hyperbolic partial differential equation shows a finite speed of propagation of heat and can represent the system as properly properly as as the the parabolic parabolic equation. equation. as differential shows aequation. finite speed of propagation of heat and can represent the system as properly equation as the parabolic as properly the parabolic equation. © 2019, IFACas(International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Distributed-parameter Distributed-parameter systems; systems; Moving Moving boundary boundary conditions; conditions; Heat Heat flows; flows; Linear Linear Keywords: Keywords: Distributed-parameter systems; Moving boundary conditions; Heat flows; Linear systems; Eigenvalue Eigenvalue problems; systems; problems; Keywords: Distributed-parameter systems; Moving boundary conditions; Heat flows; Linear systems; Eigenvalue problems; systems; Eigenvalue problems; 1. INTRODUCTION INTRODUCTION This manuscript manuscript considers considers both both hyperbolic hyperbolic and and parabolic parabolic 1. This 1. INTRODUCTION This manuscriptThe considers bothand hyperbolic and heat equations. equations. The derivation and analysis of of theparabolic numeriheat derivation analysis the numeri1. INTRODUCTION This manuscript considers hyperbolic and heat equations. The derivation and analysis theparabolic numerical results results for both both cases areboth considered in the theof next section. cal for cases are considered in next section. heat equations. Thecases derivation and analysis ofnext the numeriTransport-reaction distributed parameter system (DPS) cal results for both are considered in the section. After that, the analysis of the eigenvalues associated with Transport-reaction distributed parameter system (DPS) cal After that,for theboth analysis of the eigenvalues associated with results cases of are considered in the nextare section. Transport-reaction distributed parameter manufacturing system (DPS) After models present in chemical, petrochemical, that, the analysis the eigenvalues associated with each equation is shown. Finally, both equations conmodels present in chemical, petrochemical, manufacturing each equation is shown.ofFinally, both equations are conTransport-reaction distributed parameter manufacturing system (DPS) that, the eigenvalues associated models present in chemical, petrochemical, and process process industry take the the mathematical form given given by After each equation isanalysis shown. Finally, both equations are with considered in aathe moving boundary problem, which numerical sidered in moving boundary problem, which numerical and industry take mathematical form by models present in chemical, petrochemical, manufacturing equation is shown. Finally,problem, both equations are conand process industry take the mathematical formThe given by each hyperbolic partial differential equations (PDEs). The main sidered in ashown moving boundary which numerical results are shown in the the subsequent section. hyperbolic partial differential equations (PDEs). main results are in subsequent section. and process industry take the mathematical form given by sidered in a moving boundary problem, which numerical hyperbolic partial equations The main results are shown in the subsequent section. conservation laws differential are embedded embedded in aa (PDEs). modelling variety conservation laws are in modelling variety hyperbolic partial differential equations (PDEs). The main results are shown in the subsequent section. conservation lawshyperbolicity are embedded modelling variety provided by by the hyperbolicity of in thea transport transport systems 2. THE THE HEAT HEAT EQUATION EQUATION provided the of the systems conservation lawshyperbolicity are embedded in a transport modellingassystems variety 2. provided by the of the which is physically relevant and desired property action 2. THE HEAT EQUATION which is physically relevant and desired property as action provided byisthe hyperbolicity of the transport THE HEAT EQUATION which is physically relevant and desired property assystems action at distance distance precluded and physically meaningful finite Considering a 2. at is precluded and physically meaningful finite fluid at at rest rest with with constant constant density density and and nenewhich is physically relevant and desired property as action Considering a fluid at distance is precluded and physically meaningful finite speed of phenomena propagation is ensured. On the other Considering a fluid at rest with constant density and nespeed of phenomena propagation is ensured. On the other glecting non-linear terms in gradients and time-derivatives, at distance is precluded and physically meaningful finite glecting non-linear in gradients and time-derivatives, speed of propagation is ensured. On the hand, thephenomena hyperbolicity mathematically ensures the other well- Considering a fluidterms at with constant density and nehand, the hyperbolicity mathematically ensures the wellglecting non-linear terms in gradients and time-derivatives, the energy energy balance forrest the system leads to the following following speed of phenomena propagation is ensured. On the other the balance for the system leads to the hand, the hyperbolicity mathematically ensures the wellposedness of local Cauchy problems (Fischer and Marsden glecting non-linear terms in gradients and time-derivatives, the energy balance for the system leads to the following posedness of local Cauchy problems (Fischer and Marsden parabolic equation: hand, the of hyperbolicity ensures well- the parabolic equation: posedness local Cauchymathematically problems (Fischer and the Marsden (1972)). energyequation: balance for the system leads to the following parabolic (1972)). posedness of local Cauchy problems (Fischer and Marsden (1972)). parabolic equation: ∂T Diffusive transport transport of of heat heat across across macroscopic macroscopic length length ∂T (1972)). = α∆T α∆T (1) Diffusive = (1) ∂T ∂t = Diffusive transport of heat across macroscopic length scales is well described by Fourier’s law. However, any ∂t α∆T (1) ∂T scales is well described by Fourier’s law. However, any Diffusive transport heat across body macroscopic length ∂t = α∆T (1) scales well described by material Fourier’s law. is However, any initial is disturbance inof the the material is propagated initial disturbance in body propagated The heat heat flux flux can be given ∂t scales is well described by Fourier’s law. However, any The can be given by by the the Fourier’s Fourier’s law law and and initial disturbance in the material body is propagated instantly due due to to the the parabolic parabolic nature nature of of the the partial partial differdifferinstantly The heat flux given by Fourier’s considering onlycan one be dimension it is isthe defined as: law and initial disturbance the material propagated considering only one dimension it defined as: instantly due toobtained the in parabolic nature ofJordan theispartial differential equation equation obtained (Christov andbody Jordan (2005)). To The heat flux given by the Fourier’s law and ential (Christov and (2005)). To considering onlycan one be dimension it is defined as: instantly due to the parabolic nature of the partial differT is defined as: (2) q= = −k∂ −k∂ζitT ential equation obtained (Christov Jordanmodified (2005)).the To considering only one dimension eliminate this unphysical unphysical feature, and Cattaneo modified the (2) q ζ eliminate this feature, Cattaneo ential equation obtained (Christov and Jordan (2005)). To q = −k∂ζ T (2) eliminate this unphysical feature, Cattaneo modified the Fourier law to take into account the thermal inertia, which Fourier lawthis to take into account theCattaneo thermal inertia, which = −k∂ζ conductivity. T (2) eliminate unphysical feature, modified the where where k k is is the the material materialqthermal thermal The energy energy Fourier law to take into account the thermal inertia, which avoids the phenomenon of infinite propagation (Cattaneo conductivity. The avoids the phenomenon of infinite propagation (Cattaneo Fourier law to take into account the thermal inertia, which where k is the material thermal conductivity. The energy balance can be rewritten as: avoids phenomenon infinite (Cattaneo balance (1958);the Vernotte (1958); of Jou et al. al. propagation (2001)). can be material rewrittenthermal as: is the conductivity. The energy (1958); Vernotte (1958); Jou et (2001)). avoids phenomenon infinite (Cattaneo where balancek can be rewritten as: (1958);the Vernotte (1958); of Jou et al. propagation (2001)). ∂ T + ∂ q = 00 → ∂ T (3) t ζ tT balance can be rewritten as: ∂ T + ∂ q = → ∂ T = = α∂ α∂ζζ (3) t ζ t ζζ T Moreover, a Stefan problem is also considered, which is (1958); Vernotte (1958); Jou et Moreover, a Stefan problem is al. also(2001)). considered, which is ∂t T + ∂ζ q = 0 → ∂t T = α∂ζζ T (3) Moreover, a Stefan problem is also considered, which is a specific type of boundary value problem for a partial ∂ T + ∂ q = 0 → ∂ T = α∂ T (3) k k t= k , ζρ is the material t ζζ a specific type of boundary value problem for a partial k with α = density, c is the Moreover, aequation, Stefan problem is focusing alsoproblem considered, which is with α = ρc adifferential specific type of boundary value for a partial differential generally in the heat disk = C k , ρ is the material density, c is the ρc C equation, generally focusing in the heat dis- material with α =heat = Ck , ρ is material density, is the capacity pertheunit unit of mass mass and C Cc is is the adifferential specific in type of boundary value problem for a moving partial k capacity ρc equation, generally focusing in the dis- material per of and the tribution, phase changing medium. Since theheat with α =heat = C , ρ is the material density, c ρc capacity tribution, inequation, aa phase changing medium. Since the moving material heat capacity per unit of mass and C is the heat per unit of volume. For an initial differential generally focusing in the heat dismaterial heat capacity per unit of volume. For an initial tribution, in a phase changing medium. Since the moving interface is unknown a priori, the solution needs to take heat capacity capacity per=unit unit ofvolume. mass and C is the interface isinunknown a priori, medium. the solution needs to take material material heat per of For an initial value problem (T (ζ, t = 0) T (ζ)) in an one dimensional tribution, a phase changing Since the moving 0 problem (T (ζ, t = per 0) =unit T0 (ζ)) in an oneFor dimensional interface is unknown a priori, theof needs to take value into account account for the the determination determination ofsolution the moving moving boundary material heat capacity of volume. an initial into for the boundary value problem (T (ζ, t = 0) = in an dimensional space, solution of this PDE is given as: interface is unknown a of priori, to take infinite 0 (ζ)) into account forexample the determination ofsolution the moving position. One example thesethe problems isneeds theboundary diffusion infinite space, the the solution of T PDE is one given as: value problem (T (ζ, t =∞ 0) = Tthis in one dimensional an 0 (ζ)) position. One of these problems is the diffusion infinite space, the solution of this PDE is given 2as: into account for the determination of the moving boundary ∞ 2 1 (z − ζ) position. One example of these problems is the diffusion of heat in the melting of ice: as the melting occurs the infinite space, the solution of this PDE is given as: 1 (z − ζ) of heat inOne the example melting of of these ice: asproblems the melting occurs the ∞ T (z, 0)exp − T (ζ, t) = 2 position. is the diffusion T (ζ, t) = (4παt) dz (4) (4) 1 3/2 − ζ)2 dz −∞ of heat in (interface the melting of ice:the as solid the melting occurs the boundary (interface between the solid and liquid liquid phase) ∞ T (z, 0)exp − (z 4αt 3/2 boundary between and phase) 4αt T (ζ, t) = (4παt) dz (4) T (z, 0)exp − 1 (z − ζ) of heat in the melting of ice: as the melting occurs the −∞ 3/2 boundary (interface between the solid and liquid phase) will be changing position. Although the natural occurrence 4αt (4παt) T (ζ, t) = dz (4) T (z, 0)exp − will be changing position. Although the natural occurrence 3/2 −∞ boundary (interface between the solid and liquid phase) As long as (4παt) 4αtfrom 0, this will be changing position. Although the natural occurrence −∞ the initial condition is different of a Stefan problem is mostly associated with the melting of a be Stefan problem is mostly associated with the melting As long as the initial condition is different from 0, this will changing position. Although the occurrence long as the initial condition different fromwhich 0, this of problem is mostly associated with the melting As solution predicts instant heatispropagation, propagation, which is anda Stefan solidification problems, there arenatural some Stefan-like solution predicts aa instant heat is and solidification problems, there are some Stefan-like As longthe as the initial condition ispropagation, different fromwhich 0, this of a Stefan problem is mostly associated with the melting solution predicts a instant heat is and solidification problems, there are some Stefan-like called Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous called the Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous solution predicts a instant heat propagation, which is and solidification problems, there are some Stefan-like called the Heat Conduction Paradox. Cattaneo (1958) problems related, for instance, to the fluid flow in porous wrote a paper in which he addressed the question of the media or even shock waves in gas dynamics (Rubenstein wrote athe paper in Conduction which he addressed theCattaneo question (1958) of the media or even shock waves in gas dynamics (Rubenstein called Heat Paradox. problems related, for instance, to the fluid flow in porous wrote a paper which he addressed theFourier’s questionlaw of and the media of He (1971)).or even shock waves in gas dynamics (Rubenstein paradox paradox of heat heatinconduction. conduction. He modified modified (1971)). a paper which he addressed theFourier’s questionlaw of and the media paradox of heatinconduction. He modified Fourier’s law and (1971)).or even shock waves in gas dynamics (Rubenstein wrote paradox of heat conduction. He modified Fourier’s law and (1971)). 2405-8963 © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.
Copyright © 2019 85 Copyright 2019 IFAC IFAC 85 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 2019 IFAC 85 10.1016/j.ifacol.2019.07.011 Copyright © 2019 IFAC 85
IFAC TFMST 2019 64 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
The parabolic equation predicts an instant heat propagation (the initial slopes show that for t > 0, T (ζ, t) > 0 ∀ ζ), while the initial response for the 2nd order hyperbolic consider that there is some time-lag for the initial propagation (there is a finite propagation velocity in the system). These results can be interesting for chemical engineering process (for instance in tubular reactors Dochain (2016); Fogler (2016)), as generally the models of reaction-diffusion nature are represented as parabolic equations. Therefore, it would be interesting to study these models represented as a hyperbolic equation instead, as it would not characterize a instantaneous transport shown in parabolic partial differential equations. 3. EIGENVALUE ANALYSIS In this section, a comparison between the eigenvalues of the parabolic and second order hyperbolic is made to understand how the system dynamics are different for the two scenarios. The eigenvalue problem is defined as Aφ = λφ, where A is the operator considered. For the ∂2 parabolic scenario, A = ∂ζ 2 , with the boundary conditions considered ∂ζ x(ζ = 0) = ∂ζ x(ζ = 1) = 0. The solution for the eigenvalue problem for the parabolic equations is well known and given by: λn = −αn2 π 2 , n = 1, ..., ∞ (8) This result indicates that the eigenvalues grow unbounded. The partial differential equation shown in Eq. (7) can be written as: 0 1 T T ∂ = (9) α ∂ 2 −1 ∂T ∂t ∂T 2 τ ∂ζ τ ∂t ∂t And an operator A can be defined as:
Fig. 1. Comparison between the results from the Heat Equations: (Upper) Parabolic; (Bottom) 2nd order Hyperbolic
based his argument on the elementary kinetic theory of gases and argued that there is a time-lag between the start of the particles at their point of departure and the time of passage through the middle layer. Therefore, if the temperature changes in time, it is the heat flux at a certain time depends on the temperature gradient at an earlier time. This assumption leads to the following definition of a modified heat flux: (5) q = −k(1 − τ ∂t )∂ζ T If τ is small:
(1 − τ ∂t )−1 ≈ (1 + τ ∂t )
T ∂x A = α ∂ 2 −1 → = Ax, x = ∂T ∂t τ ∂ζ 2 τ ∂t 0
1
(10)
The eigenvalue problem is defined as Aφ = λφ, which gives the following system of equations: φ2 = λφ1 α ∂ 2 φ1 λφ1 φ2 α ∂ 2 φ1 = λ2 φ1 → − 2 = λφ − τ ∂ζ τ 2 τ ∂ζ 2 τ
(11)
Solving this second order ordinary differential equation with the boundary conditions ∂ζ φ1 (ζ = 0) = ∂ζ φ1 (ζ = 1) = 0 results in the following condition: τ λ2n + λn + αn2 π 2 = 0, n = 1, ..., ∞ (12) and for each value of n: −1 ± 1 − 4τ (αn2 π 2 ) λn = (13) 2τ
(6)
And the energy balance can be rewritten as a second order hyperbolic PDE: (7) ∂t T + ∂ζ q = 0 → τ ∂tt T + ∂t T = α∂ζζ T If τ = 0, the original parabolic PDE is obtained. Considering the heat conduction equations shown in (3) and (7) with ∂ζ T (ζ = 0) = ∂ζ T (ζ = 1) = 0, the results shown in Figure 1 are obtained for the same initial condition 2 (T0 (ζ) = 1 − (ζ−0.5) 0.252 , for 0.25 ≤ ζ ≤ 0.75; 0, otherwise) and time interval.
And this results in a very different set of eigenvalues when compared to the parabolic equation. This system does not have a infinite set of real fast eigenvalues as shown before (for the parabolic equation, λn → −∞ as n → ∞), in fact, the whole set of eigenvalues needs to be in the region 86
IFAC TFMST 2019 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
Liquid TI (0 ≤ ζ ≤s, t)
1 where −1 τ < Re(λn ) < 0. Furthermore, when n > 4τ απ , the eigenvalues start to present imaginary parts, which is different from the eigenvalues obtained for the parabolic equation that have just real parts. 4. HEAT EQUATION EQUATION AND THE STEFAN PROBLEM To analyze the difference between the results obtained by a hyperbolic and a parabolic partial differential equation, a heat diffusion problem with phase change is considered. One of the simplest mathematical model of the change of phase phenomenon is called Stefan problem. If a phase change of a material occurs at a given point latent heat is released or absorbed and the temperature of the material at that point remains constant. It is assumed that the temperature of phase change is constant and that the temperature across the material is continuous. For simplicity, an one dimensional case is considered in a material with constant cross sectional area. If the phase change happens at a point ζ = s(t), then: (14) TI (ζ = s(t) − δ, t) = TII (ζ = s(t) + δ, t) = TC where δ is a small spatial variation such that δ → 0, TC is the phase change temperature and the subscripts I and II are used to represent two different phases. Therefore, this condition implies that the temperature close to the phase transition point is continuous and is the same in both phases. The Stefan condition is obtained by applying a energy balance between two instants t0 and t1 in the volume of a material that is going through phase change. If the phase change point moved from s(t0 ) to s(t1 ) and that s(t1 ) > s(t0 ), the energy involved in the phase change is given by the volume of the material that went through phase change between this two instants and the latent heat required for the transition: (15) Q = Cl A[s(t1 ) − s(t0 )] where Q is the energy involved in the phase transition, Cl is the latent heat (per unit of volume) and A is the constant cross sectional area. This amount of energy should be equivalent to the heat provided by heat diffusion in both sides of the phase change interface in this time interval: t1 t1 (qI − qII )dAdt = A (qI − qII )dt (16) Q= t0
A
t0
Cl [s(t1 ) − s(t0 )] =
t0
TII (s ≤ ζ ≤1, t) = 0
ζ =1 T =0
Fig. 2. Representation of the system with phase transition ature of T0 (ζ) = 0 is considered throughout the system. As the phase change temperature considered is the same as the temperature at the boundary ζ = 1 and the initial temperature in the solid phase, there is no heat flux that flows through this phase (as ∂ζ TII (s≤ζ ≤1, t) = 0 and no heat sinks or sources are considered in 0<ζ <1). This simplifies the model to a so called one phase problem, as only one phase (in this case the liquid phase) has a temperature that is time and space dependent that needs to be found (TI (0≤ζ ≤s, t)). Therefore, the Stefan condition for the moving boundary for the original Fourier’s law can be simplified to: dt s = qI = −k ∂ζ T |ζ = s(t) → dt s = −β ∂ζ T |ζ=s(t)
(19)
where, for simplicity of the notation, T is now used to represent TI and β is defined as β = Ckl , with k as the material thermal conductivity and Cl the latent heat per unit of volume, defined previously. If the modified heat flux shown in Eq. (5) is used, the Stefan condition obtained is: dt s = qI = −k(1 − τ ∂t ) ∂ζ T |ζ=s(t) → τ dtt s + dt s = −β ∂ζ T |ζ=s(t)
(20)
This results in a second order ordinary differential equation for the modified heat flux, while the original Fourier’s law results in a first order ordinary differential equation. These equations need to be solved along with the corresponding hyperbolic and parabolic diffusion equations (Eq. (7) and Eq. (3), respectively).
t0
t1
Solid
ζ =0 ζ =s T =1 T =0
where qI and qII represent the heat flux from both phases and the fact that the cross sectional area does not change was used to simplify the integral. Then, t1 (qI − qII )dt → Q = Cl A[s(t1 ) − s(t0 )] = A
65
As these partial differential equations have a moving boundary, the following change of variables is used to change from a time-varying domain to a fixed domain: ζ → ∈ [0, 1] (21) = s(t)
(17)
(qI − qII )dt
Dividing the equation by t1 −t0 and setting t1 → t0 , finally gives the Stefan condition: (18) Cl dt s(t) = qI (t, ζ = s − δ) − qII (t, ζ = s + δ)
And the derivatives can be rewritten for this new spatial variable as: 1 ∂ ∂T ( , t) = ∂ T ( , t) (22) ∂ζ T (ζ, t) = ∂ζ ∂ s
The system considered here is shown in Figure 2, which represents melting ice with fixed temperatures at ζ = 0 and ζ = 1 ( T (ζ = 0) = 1 and T (ζ = 1) = 0).
∂ ∂ζζ T (ζ, t) = ∂
The transition temperature (melting temperature) of the solid phase is considered to be TC = 0. A initial temper87
∂ ∂T ( , t) ∂ζ ∂
1 ∂ = 2 ∂ T ( , t) ∂ζ s
(23)
IFAC TFMST 2019 66 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67 Louvain-la-Neuve, Belgium, July 3-5, 2019
∂t T (ζ, t) = ∂t T (, t) +
∂ ∂T (, t) = ∂t ∂
∂t T (, t) − dt s∂ T (, t) s
(24) 1 0.8
2 ∂tt T (, t) − dt s∂ ∂t T (, t)+ s 2 ∂ T (, t) 2 (dt s)2 − dtt s s s
Temperature
∂tt T (ζ, t) = ∂ ∂T (, t) ∂ ∂T (, t) ∂ 2 ∂tt T (, t) + 2 + = ∂ ∂t ∂t ∂ ∂t2
0.6 0.4 0.2 0 0.6 0 0.4
0.2 0.4 0.2
0.6 0.8
Time
(25)
Finally, the parabolic partial differential equation can be written as: 1 (26) α 2 ∂ T (, t) = ∂t T (, t) − dt s∂ T (, t) s s And for the hyperbolic: 1 α 2 ∂ T (, t) = s 2 τ ∂tt T (, t) − dt s∂ ∂t T (, t)+ s (27) 2 2 ∂ T (, t) 2 (dt s) − dtt s + s s ∂t T (, t) − dt s∂ T (, t) s
0
1
Fig. 3. Spatial profile obtained for the Stefan problem using a parabolic heat equation (s(t = 0) = 0.5)
1
Temperature
0.8
0.6
0.4
0.2
0 0.6 0.4
The Stefan conditions also need to be properly modified. The first order is changed to: 1 dt s = −β ∂ T (, t) (28) s While the second order is given by: 1 τ dtt s + dt s = −β ∂ T (, t) (29) s The boundary conditions also need to be redefined: u(ζ = 0, t) = 1 → u( = 0, t) = 1 (30) u(ζ = s(t), t) = 0 → u( = 1, t) = 0
0.2
Time
0
1
0.8
0.6
0.4
0.2
0
Fig. 4. Spatial profile obtained for the Stefan problem using a hyperbolic heat equation interface position at s(t = 0) = 0.5 was used to show the difference on the initial diffusion between the two types of partial differential equation. The results for α = β = 1 and τ = 10−3 are shown in Figure 3 and Figure 4 for the parabolic and hyperbolic diffusion equation, respectively. Both results show that initially there is a time for which the boundary does not move. This happens because the initial condition in the liquid is the same as the transition temperature and it takes some time for the heat to flow through the liquid phase. As the diffusion starts to happen and the heat reaches the interface between the two phases, the boundary starts to move (the solid phase starts to melt) and this happens until s(t) = 1. When this condition is reached, the boundary stops to move and the system is represented just by diffusion in a fixed-domain. After some time, the system reaches the expected steady-state.
For the second order hyperbolic equation, boundary conditions are also necessary for the first order time derivatives, which can be obtained directly from the boundary conditions: u( = 0, t) = 1 → ∂t u( = 0, t) = 0 (31) u( = 1, t) = 0 → ∂t u( = 1, t) = 0 To guarantee that s(t) ≤ 1, a condition is imposed to stop the growth of the moving boundary. If s(t) > 1, the first and second derivative of s(t) are set to zero. This makes Eq. (26) and Eq. (27) convert back to the original diffusion problem (Eq. (3) and Eq. (7), respectively), as expected if there is no more change in the domain.
The main difference between the results for the parabolic and the second order hyperbolic can be seen at the initial diffusion. As shown before, in the analysis of the results of the heat equation (Figure 1), the hyperbolic equation shows a finite velocity of propagation, while the parabolic equation represents a situation of instant diffusion. These results can be seen again in Figure 3 and Figure 4. For the parabolic equation, if t > 0 all points between 0 ≤ ζ ≤ s(t) present a temperature different from zero. For the
5. RESULTS This section presents the results for the hyperbolic and parabolic partial differential equations (and the respective second order and first order ordinary differential equations) with the conditions previously chosen. A initial 88
IFAC TFMST 2019 Louvain-la-Neuve, Belgium, July 3-5, 2019 Guilherme Ozorio Cassol et al. / IFAC PapersOnLine 52-7 (2019) 63–67
1
1
0.9
0.9 0.8 Par( =1, =1)
0.7
-3
Cat( =1, =1, =10 )
s (Moving boundary)
s (Moving boundary)
0.8
67
0.6 0.5 0.4 0.3
0.7 0.6 0.5
Par( =1, =1)
0.4
Cat( =1, =1, =10 -3 ) Par( =1, =10)
0.3
Cat( =1, =10, =10 -3 ) Par( =0.1, =10)
0.2
Cat( =0.1, =10, =10 -3 ) Par( =0.1, =1)
0.2 0.1
Cat( =0.1, =1, =10 -3 )
0.1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0
Time
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time
Fig. 5. Moving boundary dynamics for the Stefan problem using a parabolic and hyperbolic heat equation (s(t = 0) = 0.5)
Fig. 6. Moving boundary dynamics for the Stefan problem using a parabolic and hyperbolic heat equation for different values of α and β (s(t = 0) = 0)
hyperbolic, the initial slope after t > 0 shows the finite speed of propagation. Even though the initial diffusion is different, both hyperbolic and parabolic present a similar growth of the boundary, as it is shown in Figure 5 for s(t = 0) = 0.5. The difference between the partial differential equations is barely noticeable and the hyperbolic equation seems to be also adequate to model the system.
ACKNOWLEDGEMENTS The support for this work is provided by CAPES 88881.128514/2016-01 (Brazil) and support for Guilherme Ozorio Cassol is gratefully acknowledged. REFERENCES
Figure 6 shows the results for s(t = 0) = 0 using different values of α and β. As the value of α increases, faster is the growth of the moving boundary, as the heat gets to the interface of the two phases in a faster pace. This can be seen if the red line is compared to the black line, and the blue line with the green line. A faster growth of the boundary is also seem if β. The higher the value of this parameter, the smaller is the heat necessary for phase transition, which makes the boundary moves faster. This can be seen if the red line is compared to the green line and the black line with the blue line. Also, a decrease in one of these parameters is not compensated if the other is increased, as it can be conclude if the green and blue lines are compared. For all cases it is possible to see that the difference between the hyperbolic and parabolic equations is barely noticeable. Once again, the hyperbolic is adequate to represent the system.
Cattaneo, C. (1958). On a form of heat equation which eliminates the paradox of instantaneous propagation. C. R. Acad. Sci. Paris, 431–433. Christov, C.I. and Jordan, P.M. (2005). Heat conduction paradox involving second-sound propagation in moving media. Physical review letters, 94 15, 154–301. Dochain, D. (2016). Analysis of the multiplicity of equilibrium profiles in tubular reactor models. IFACPapersOnLine, 49(18), 903–908. Fischer, A.E. and Marsden, J.E. (1972). The Einstein evolution equations as a first-order quasi-linear symmetric hyperbolic system, I. Communications in Mathematical Physics, 28(1), 1–38. Fogler, H.S. (2016). Elements of chemical reaction engineering. Prentice Hall, 5 edition edition. Jou, D., Casas-Vazquez, J., and Lebon, G. (2001). Extended irreversible thermodynamics (3rd Edition). Springer. Rubenstein, L. (1971). The Stefan problem, by L.I. Rubenstein. Translations of mathematical monographs. American Mathematical Society. Vernotte, P. (1958). Les paradoxes de la theorie continue de l’equation de la chaleur. C. R. Acad. Sci. Paris, 246, 3154–3155.
6. CONCLUSION This manuscript compared the results obtained for the heat diffusion considering hyperbolic and parabolic equations. First, the simple case of one-dimensional heat diffusion was considered and the eigenvalue analysis was performed to access the difference between these two equations, which shows the finite speed of propagation of the hyperbolic equation. Finally, the Stefan problem was considered and the results show that the interface dynamic for the hyperbolic equation are similar to the parabolic in the cases analyzed, but the former finite speed of propagation would be more desirable to represent the dynamics of the actual physical system. 89