Eindhoven University of Technology BACHELOR The hotspot problem van Roosmalen, AH

37 Pages·2017·2.3 MB·English

Eindhoven University of Technology
BACHELOR
The hotspot problem
van Roosmalen, A.H. (Teun)
Award date:
2017
Link to publication
Disclaimer
This document contains a student thesis (bachelor's or master's), as authored by a student at Eindhoven University of Technology. Student
theses are made available in the TU/e repository upon obtaining the required degree. The grade received is not published on the document
as presented in the repository. The required complexity or quality of research of student theses may vary by program, and the required
minimum study period may vary in duration.
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.
• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
Eindhoven University of Technology
The Hotspot Problem
Author:
A.H. van Roosmalen
[email protected]
0900962 Supervisor:
dr.ir. J.H.M. ten Thije
Boonkkamp
[email protected]
June 19, 2017
The Hotspot Problem
Abstract
In this paper I will work with the model of an exothermic chemical reaction in some vessel. This
consists of a partial differential equation which can not be solved exact. Therefore a numerical
integration method needs to be used to solve this problem numerically. In this report I will use
the implicit Euler and exponential Euler method. I will also derive and use an improved version
ofthelatter. Eventhoughthebackgroundofthisproblemliesinthefieldofchemistry, thispaper
will focus mainly on the mathematical approach of such a model. The purpose of this study is to
gain a deeper insight of the use of numerical mathematics to solve partial differential equations. I
also want to research the use of exponential time-integration methods.
/department of mathematics and computer science 1
The Hotspot Problem
Contents
1 Introduction 3
2 The Hotspot Model 5
2.1 Construction of the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 The Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Central-Difference Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3 Conventional Integration Methods 12
3.1 Explicit Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Implicit Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4 Exponential Integration Methods 19
4.1 Exponential Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Improved Exponential Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.3 Analysis of Exponential Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5 Conclusion 27
A MATLAB Code 28
/department of mathematics and computer science 2
The Hotspot Problem
1 Introduction
Energy,heatandconsequentlyalsotemperatureareimportantpartsofoureverydaylife. Animals
have an instinct for keeping warm and ever since the invention of fire, people have been using
its warmth for different purposes. The first person to create a reliable thermometer was Gabriel
DanielFahrenheitin1724[8]. Thisgavepeoplethepossibilitytoassignquantitiestotemperature
and therefore apply operations to them. This gave the option to apply mathematics to temper-
ature. Heat diffusion turned out to be of great interest and in 1807, Fourier first stated his heat
equation[8]. It became one of the most famous and most used partial differential equations. In its
autonomousversionitiseasytounderstandandwidelyusedtoteachaboutdifferentialequations.
One of the things influenced by temperature is the rate at which chemical reactions occur.
Of almost all reactions, the reaction rate greatly depends on the temperature of the affected
substances. Many companies in the modern times use chemical reactions to manufacture their
products. It is of major concern for them to know how the temperature progresses through time.
One of the reasons for those in the chemical industry to study this is safety. Unforeseen increases
or decreases in temperature are unwanted. Rapidly increasing temperature can even cause explo-
sions. Being able to predict this behaviour gives those companies involved the opportunity to act
accordingly.
In this report I will take a look specifically at an exothermic chemical reaction. This is an
reaction that releases energy in the form of light or heat, next to the reaction products. In this
case the reaction the released energy is heat. This reaction is indeed influenced by the tempera-
ture, but also acts on the temperature itself. The process associated with these kinds of reactions
is called thermal runaway or thermal explosion. However fascinating this occurrence is, it is also
very dangerous. In 1947 two ships in Texas City exploded when the ammonium nitrate on board
heated enough to start such a thermal explosion. It resulted in at least 568 deaths and hundreds
of millions of dollars material damage[9].
Figure 1: Abandoned street, showing extensive damages after the Texas City disaster
[12].
Iwilltakealookatsuchanequationinsidea1by1domain. Twoofthesideswillbeisolated,
so no heat gets lost or gained through interaction with the outside world. The other two sides
willhaveaconstanttemperature, eitherbyusingcoolingorheating. Thisgivesinterestingresults
/department of mathematics and computer science 3
The Hotspot Problem
that I will inspect through the use of different time integration methods. These methods can only
be used if I discretise space in this problem first. This will be done by using central-difference
schemes. I will then first try to use the most basic integration methods: the implicit and explicit
Eulermethods. Next,twoslightlymoreadvancedexponentialintegrationschemeswillbeapplied.
This report will focus mainly on the mathematical part of this situation and therefore I will use a
non-dimensionalized equation further on. This slightly disconnects the equation from reality, but
keeps things clearer with less constants.
The paper is divided into sections as follows. Section 2 will contain a construction of the
modelandanexplanationofthecentral-differencescheme. Thiswillalsobeappliedtothepartial
differential equation to create a system of linear ordinary differential equations. In section 3 I will
use the implicit and explicit Euler methods on the obtained system. I will then elaborate on any
potentialproblemthatoccurs. Section4willbesomewhatthesameexceptthatIuseexponential
integration methods in this part. These methods will first be derived and then examined. Section
5 will be the conclusion of this report.
InthisreportIwillusenormalLatinorGreeklettersforscalars,boldlowercaselettersforvectors,
and bold uppercase letters for matrices. For example a and α are scalars, x is a vector and A is
a matrix. Whenever a norm is used, the standard Euclidean norm or two-norm is meant unless
specified otherwise.
/department of mathematics and computer science 4
The Hotspot Problem
2 The Hotspot Model
2.1 Construction of the Model
The model I will work with in this report describes the temperature of a reacting mixture in a
vessel. This vessel is two dimensional, given by x and y coordinates in the Euclidean plane. The
model consists of two parts. The first part is a self-accelerating reaction. We can assume the
vessel to contain a perfectly stirred mixture. This makes the temperature independent on the
locationinthevessel. Thetemperaturewillbegivenbythefunctionu=u(t),wheretistime. As
brieflymentionedbeforeintheintroduction,therateofreactiondependsonthetemperature. The
reaction itself is exothermic. This means that energy is released, in this case in the form of heat.
Combining these properties of this reaction leads to a situation where an increase of temperature
influences the reaction rate and thereby increases the temperature. The reaction can lead to a
thermal runaway. This property is described by the differential equation
du R
= (1+α−u)eδ(1−1/u) =:s(u). (2.1.1)
dt αδ
This equation is also mentioned by Hundsdorfer and Verwer in section V.3.1 [6].
The parameters R, δ and α will be specified later on in this report. The factor (1+α−u) can
be seen as the amount of fuel in the vessel. If the temperature rises, the reaction will use more
fuel, so less of it will be left. This means that α is a parameter connected to the initial amount of
fuelinthereactor. Theexponentinequation(2.1.1)describestheinfluenceofthetemperatureon
the reaction rate. A larger δ means more influence of the temperature. Finally R is a parameter
of the reaction rate.
Another important part of this model is the diffusion of heat. For this, the vessel no longer
contains a perfectly stirred mixture. I will describe the situation without a reaction at all. The
temperature is now dependent on the location, so it is given by u = u(x,t) with x ∈ Rn. The
equation to model this diffusion is the heat equation, which reads as follows
∂
u=d∆u, (2.1.2)
∂t
where ∆ is the Laplace operator and d > 0 is the diffusion coefficient. This depends on the
material that the heat spreads through. We assume the vessel to contain a homogeneous content,
so that this diffusion coefficient is constant. A thorough derivation of the equation can be found
in ’The Heat Equation’ by D. Widder [11]. In our two-dimensional case, the temperature is given
by u=u(x,y,t). The heat equation boils down to the following:
∂u (cid:18)∂2u ∂2u(cid:19)
=d + . (2.1.3)
∂t ∂x2 ∂y2
To simplify the way of writing things, I will from now on denote the derivatives with a subscript.
Now that the different effects on the temperature are modelled and described with equations,
they can be combined. Since they both alter the rate of change of temperature through time, the
right hand sides of (2.1.1) and (2.1.3) have to be added to create this equation:
u =d(u +u )+s(u), (2.1.4)
t xx yy
where s(u) is defined in equation (2.1.1).
Now that the equation describing the behaviour of the temperature is specified, the vessel needs
to be modelled. In this report a container of non-dimensional length and width of 1 will be used.
The third dimension, height, will not be taken into account, so basically we have the unit square
/department of mathematics and computer science 5
The Hotspot Problem
[0,1]2 ⊂ R2. The actual vessel could be very high. This would mean that there is no variation
of temperature in the z-direction, except for near the top and bottom. A two-dimensional vessel
models a cross section of such a three-dimensional container. For the model it helps a great deal
todescribethisproblemintwoinsteadofthreedimensions. Anadvantageisvisuallyrepresenting
the results. With a two-dimensional area of interest, a third dimension can be used to represent
thetemperature. Thisallowsthecreationofthreedimensionalsurfaceplots. Themainadvantage,
however, is simplicity. It is easier to calculate in only two dimensions instead of three. Especially
formorecomplicatedproblems,theCPUtimeandmemoryneededforsimulationsaresignificantly
lower in two dimensions.
On the boundaries of this container additional conditions are needed. For two neighbouring
sides the condition that the temperature is always 1 is added. This could be due to cooling or
heating elements at those sides. The sides where either the x- or y-coordinate is 1 are chosen for
this. At the other two sides the heat flux in the outward direction is 0. This is the result of a
perfectisolation. Aremarkhastobemadethatthisdoesnotmeanthatthetemperatureatthose
boundaries does not change. It can still change in time. These boundary conditions are written
like this
u (0,y,t)=0, u(1,y,t)=1, 0<y <1,t>0, (2.1.5a)
x
u (x,0,t)=0, u(x,1,t)=1, 0<x<1,t>0. (2.1.5b)
y
y axis
u=1
1
u =0 u=1
x
0 u =0 1 x axis
y
Figure 2: A sketch of the domain with the boundary conditions.
The only thing left now for this model is the initial condition. This will be u(x,y,0) = 1
everywhere in the square. Now, combining everything leads to the following initial and boundary
value problem (IBVP):
u =d(u +u )+s(u), 0<x,y <1,t>0, (2.1.6a)
t xx yy
u(x,y,0)=1, 0<x,y <1, (2.1.6b)
u (0,y,t)=0, u(1,y,t)=1, 0<y <1,t>0, (2.1.6c)
x
u (x,0,t)=0, u(x,1,t)=1, 0<x<1,t>0. (2.1.6d)
y
For all the parameters, the following values will be assigned: d=1,α=1,δ =20 and R=5. As
we will see later on in this report, the value of δ will play an interesting role in this problem.
/department of mathematics and computer science 6
The Hotspot Problem
2.2 The Source Term
Thesourceterminthisproblemisthemainpartthatmakesitdifficulttostudy. Withoutit,only
diffusion would be taken into account. This would result in equation (2.1.3). I could work this
out, for example with the help of Fourier transforms. However, this might not be necessary if I
just take a look at the situation.
Since the initial condition states that u=1 at t=0 for all x and y with 0≤x,y ≤1, and the
boundary conditions do not imply any change in temperature, there would be no diffusion. This
means that without the source term, u(x,y,t)=1 for all x, y and t.
Now the importance of the source term is clear. A good way to start might be to take a closer
look at this function. Substituting R and α gives us the following function
5
s(u)= (2−u)eδ(1−1/u). (2.2.1)
δ
The reason that I will not yet substitute δ, is that I want to investigate its role in the IBVP. This
role will be clear later on. We can see that the only root of this function is at u=2, determined
by the factor (2−u) in the function. It is also easy to see that it is positive for u < 2 and
negative for u>2. The maximum of the function is however greatly dependent on the exponent
in this expression. This exponent contributes to a very large derivative around this root u = 2.
To illustrate this I have plotted s for various values of δ. The results can be seen in Figure 3.
Figure 3: The funcion s on the interval 0.9<u<2.1 for various different values of δ.
The enormous differences between δ = 1 and δ = 20 are well illustrated here. The function s
lookslikeanalmoststraightlineforδ =1andevenwhenδ =10therearebarelyanyfluctuations.
The differences in temperature close to u = 2 are enormous for δ = 20, however. As we can see
/department of mathematics and computer science 7
The Hotspot Problem
from the figure, around u = 1.8 we have s(u) = 300, but at u = 2.1 the value of the function is
even smaller than -700. For now it is enough to know the behaviour of this function. Later we
will come back to this and see the effects of these great differences.
2.3 Central-Difference Scheme
To approach this problem numerically, the 1-by-1 square that is our area of interest must be
coveredbyagridofpoints. Thiswayasetofinfinitelymanypairsofcoordinatesisapproximated
byasetoffinitelymanypoints. Ahorizontallinewillbesplitinnpoints,letssay{x ,x ,...,x }
1 2 n
with 0 = x < x < ··· < x < x = 1. It is logical and simpler to evenly space the points.
1 2 n−1 n
This results in a distance between two points, called the grid size, of ∆x= 1 . Notice that the
n−1
capitaldeltaisnottheLaplaceoperatorlikebefore. Inthiscaseitisusedtodenotethedifference
in x. Likewise a vertical line can be split in m points {y ,y ,...,y } with 0 = y < y < ··· <
1 2 m 1 2
y <y =1. And the grid size for y is ∆y = 1 .
m−1 m m−1
The x- and y-coordinates are split into grids using one label each. A point on the grid in the
squarethereforeneedstwolabels. Alabeliforthex-coordinateandalabelj forthey-coordinate.
A point on the grid is given by xi,j =(xi,yj). A sketch of this grid:
y axis
(1,m) (n,m)
j .. ..
. .
(1,2) ... (n,2)
(2,2)
(1,1) (2,1) ... (n,1) x axis
i
Figure 4: An illustration of the labels for the grid.
A central-difference method can now be used to approximate the spatial derivatives u and
xx
u . For the derivation of this method I work along the lines of Vuik et al.[10]. Say we want to
yy
approximate the second derivative f(cid:48)(cid:48)(x) of a function f(x). Then the central-difference formula
for its approximation Q(∆x) is give by
α f(x−∆x)+α f(x)+α f(x+∆x)
Q(∆x)= −1 0 1 , (2.3.1)
∆x2
/department of mathematics and computer science 8