14 |
\section{Steady-state Heat Refraction} |
\section{Steady-state Heat Refraction} |
15 |
\label{STEADY-STATE HEAT REFRACTION} |
\label{STEADY-STATE HEAT REFRACTION} |
16 |
|
|
17 |
In this chapter we will learn how to handle more complex geometries. |
In this chapter we show how to handle more complex geometries. |
18 |
|
|
19 |
Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex models. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modeling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in \esc to define properties like material constants and source locations. |
Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in \esc to define properties like material constants and source locations. |
20 |
|
|
21 |
We will introduce in two stages. First we will look at a very simple geometry. In fact, we will solve |
We proceed in this chapter by first looking at a very simple geometry. In fact, we solve |
22 |
the steady-state heat equation over a rectangular domain. This is not a very interesting problem but in a second |
the steady-state heat equation over a rectangular domain. This is not a very interesting problem but then we extend our script by introducing an internal, curved interface. |
|
step we will extend our script by introducing an internal, curved interface. |
|
23 |
|
|
24 |
\begin{figure}[ht] |
\begin{figure}[ht] |
25 |
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} |
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} |
27 |
\label{fig:pycad rec} |
\label{fig:pycad rec} |
28 |
\end{figure} |
\end{figure} |
29 |
|
|
30 |
\section{Example 4: Creating the model with \pycad} |
\section{Example 4: Creating the domain with \pycad} |
31 |
\sslist{example04a.py} |
\sslist{example04a.py} |
32 |
|
|
33 |
We will modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: firstly we will look the steady state |
We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look the steady state |
34 |
case with slightly modified boundary conditions and secondly we will use a more flexible tool |
case with slightly modified boundary conditions and use a more flexible tool |
35 |
to generate to generate the geometry. Lets look at the geometry first. |
to generate to generate the geometry. Lets look at the geometry first. |
36 |
|
|
37 |
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and |
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and |
38 |
insulation to both sides we want to calculate the steady-state temperature distribution. |
insulation to both sides we want to calculate the steady-state temperature distribution. |
39 |
|
|
40 |
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; |
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; |
41 |
the ones we will use are: |
the ones we use are: |
42 |
\begin{python} |
\begin{python} |
43 |
from esys.pycad import * |
from esys.pycad import * |
44 |
Point() #Create a point in space. |
Point() #Create a point in space. |
87 |
from esys.pycad.gmsh import Design |
from esys.pycad.gmsh import Design |
88 |
d=Design(dim=2, element_size=200*m) |
d=Design(dim=2, element_size=200*m) |
89 |
\end{python} |
\end{python} |
90 |
The argument \verb|dim| defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verb|element_size| defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be choose with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m |
The argument \verb|dim| defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verb|element_size| defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be chosen with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m |
91 |
and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So will end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily. |
and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So we end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily. |
92 |
We can easily add our domain \verb|rec| to the \verb|Design|; |
We can easily add our domain \verb|rec| to the \verb|Design|; |
93 |
\begin{python} |
\begin{python} |
94 |
d.addItem(rec) |
d.addItem(rec) |
95 |
\end{python} |
\end{python} |
96 |
We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this |
We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this |
97 |
but \pycad offers a more convenient technique called tagging. With this techniques items in the domain are |
but \pycad offers a more convenient technique called tagging. With this technique items in the domain are |
98 |
named using the \verb|PropertySet| class. We can then later use this name to set values secifically for |
named using the \verb|PropertySet| class. We can then later use this name to set values secifically for |
99 |
those sample points located on the named items. Here we name the bottom face of the |
those sample points located on the named items. Here we name the bottom face of the |
100 |
domain where we will set the heat influx: |
domain where we will set the heat influx: |
116 |
# write domain to an text file |
# write domain to an text file |
117 |
domain.write("example04.fly") |
domain.write("example04.fly") |
118 |
\end{python} |
\end{python} |
119 |
and then for recovery in another script; |
and then for reading in another script; |
120 |
\begin{python} |
\begin{python} |
121 |
# read domain from text file |
# read domain from text file |
122 |
from esys.finley import ReadMesh |
from esys.finley import ReadMesh |
131 |
|
|
132 |
Before we discuss how we solve the PDE for the |
Before we discuss how we solve the PDE for the |
133 |
problem it is useful to present two additional options of the \verb|Design|. |
problem it is useful to present two additional options of the \verb|Design|. |
134 |
They allow the user accessing the script which is used to by \gmsh to generate the mesh as well as |
They allow the user accessing the script which is used by \gmsh to generate the mesh as well as |
135 |
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: |
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: |
136 |
\begin{python} |
\begin{python} |
137 |
d.setScriptFileName("example04.geo") |
d.setScriptFileName("example04.geo") |
139 |
\end{python} |
\end{python} |
140 |
Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and |
Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and |
141 |
the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. |
the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. |
142 |
Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesher can be visualise with \gmsh from the command line using |
Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using |
143 |
\begin{verbatim} |
\begin{verbatim} |
144 |
gmsh example04.geo # show geometry |
gmsh example04.geo # show geometry |
145 |
gmsh example04.msh # show mesh |
gmsh example04.msh # show mesh |
155 |
\section{The Steady-state Heat Equation} |
\section{The Steady-state Heat Equation} |
156 |
\sslist{example04b.py, cblib} |
\sslist{example04b.py, cblib} |
157 |
|
|
158 |
Steady state is reached in the temperature is not changing in time. So to calculate the |
Steady state is reached in the temperature when it is not changing in time. So to calculate the |
159 |
the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero; |
the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero; |
160 |
\begin{equation}\label{eqn:Tform nabla steady} |
\begin{equation}\label{eqn:Tform nabla steady} |
161 |
-\nabla \cdot \kappa \nabla T = q\hackscore H |
-\nabla \cdot \kappa \nabla T = q\hackscore H |
162 |
\end{equation} |
\end{equation} |
163 |
This PDE is easier then the PDE in \refEq{eqn:hdgenf2} we needed to solve in each time step. We can just drop |
This PDE is easier to solve, than the PDE in \refEq{eqn:hdgenf2} in each time step. We can just drop |
164 |
the \verb|D| term; |
the \verb|D| term; |
165 |
\begin{python} |
\begin{python} |
166 |
mypde=LinearPDE(domain) |
mypde=LinearPDE(domain) |
172 |
x=Solution(domain).getX() |
x=Solution(domain).getX() |
173 |
mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop) |
mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop) |
174 |
\end{python} |
\end{python} |
175 |
Notice that we have use the \verb|sup| function to calculate the maximum $y$ coordinator of the relevant sample points. |
Notice that we use the \verb|sup| function to calculate the maximum of $y$ coordinates of the relevant sample points. |
176 |
|
|
177 |
In all cases so far we have assumed that the domain is insulated which translates |
In all cases so far we have assumed that the domain is insulated which translates |
178 |
into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling |
into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling |
185 |
where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero |
where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero |
186 |
for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face. |
for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face. |
187 |
Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. |
Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. |
188 |
One could use the masking techniques we have already extensively used to define $q\hackscore{S}$ in \esc |
We could define $q\hackscore{S}$ by using the masking techniques demonstrated earlier. The tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise |
|
but the tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise |
|
189 |
constant functions such as $q\hackscore{S}$. You need to remember now that we |
constant functions such as $q\hackscore{S}$. You need to remember now that we |
190 |
marked the bottom face with the name \verb|linebottom| when we defined the domain. |
marked the bottom face with the name \verb|linebottom| when we defined the domain. |
191 |
We can use this now to create $q\hackscore{S}$; |
We can use this now to create $q\hackscore{S}$; |
233 |
Our heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependent upon the material properties and shape. |
Our heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependent upon the material properties and shape. |
234 |
|
|
235 |
|
|
236 |
We will modify the example of the previous section by subdividing the block into two parts. The curve |
We modify the example of the previous section by subdividing the block into two parts. The curve |
237 |
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points |
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points |
238 |
used to define the curve may be measurement but for simplicity we assume here that there coordinates are |
used to define the curve may be measurement but for simplicity we assume here that there coordinates are |
239 |
known in analytic form. |
known in analytic form. |
242 |
|
|
243 |
It is now possible to start defining our domain and boundaries. |
It is now possible to start defining our domain and boundaries. |
244 |
|
|
245 |
The curve defining our clinal structure is located in approximately the middle of the domain and has a sinusoidal shape. We define the curve by generating points at discrete intervals; $51$ in this case, and then create a smooth curve through the points using the \verb Spline() function. |
The curve defining our clinal structure is located approximately in the middle of the domain and has a sinusoidal shape. We define the curve by generating points at discrete intervals; $51$ in this case, and then create a smooth curve through the points using the \verb Spline() function. |
246 |
\begin{python} |
\begin{python} |
247 |
# Material Boundary |
# Material Boundary |
248 |
x=[ Point(i*dsp\ |
x=[ Point(i*dsp\ |
275 |
#clockwise check |
#clockwise check |
276 |
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) |
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) |
277 |
\end{python} |
\end{python} |
278 |
The last few steps in creating the model take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc. |
The last few steps in creating the model are: take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc. |
279 |
To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements that need to be applied to the domain \verb element_size . It then becomes a simple task of adding the sub-domains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the sub-domain properties in the solution script. This is done using the \verb PropertySet() function. The geometry and mesh are then saved so the \esc domain can be created. |
To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements that need to be applied to the domain \verb element_size . It then becomes a simple task of adding the sub-domains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the sub-domain properties in the solution script. This is done using the \verb PropertySet() function. The geometry and mesh are then saved so the \esc domain can be created. |
280 |
\begin{python} |
\begin{python} |
281 |
# Create a Design which can make the mesh |
# Create a Design which can make the mesh |
313 |
\sslist{example05b.py and cblib.py} |
\sslist{example05b.py and cblib.py} |
314 |
|
|
315 |
We want to investigate the profile of the data of the last example. |
We want to investigate the profile of the data of the last example. |
316 |
We are particularly interested in the deepth profile of the heat flux which is |
We are particularly interested in the depth profile of the heat flux which is |
317 |
the second component of $-\kappa \nabla T$. We will extend the script developed in the |
the second component of $-\kappa \nabla T$. We extend the script developed in the |
318 |
previous section to show how for instance vertical profile can be plotted. |
previous section to show how for instance vertical profile can be plotted. |
319 |
|
|
320 |
The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and |
The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and |
322 |
the flux is the product of the piecewise constant function $\kappa$ and |
the flux is the product of the piecewise constant function $\kappa$ and |
323 |
the gradient of the temperature $T$ which has a kink across the rock interface. |
the gradient of the temperature $T$ which has a kink across the rock interface. |
324 |
Before plotting this function we need to smooth it using the |
Before plotting this function we need to smooth it using the |
325 |
\verb|Projector()| factory; |
\verb|Projector()| class; |
326 |
\begin{python} |
\begin{python} |
327 |
from esys.escript.pdetools import Projector |
from esys.escript.pdetools import Projector |
328 |
proj=Projector(domain) |
proj=Projector(domain) |
335 |
\begin{python} |
\begin{python} |
336 |
xiq,yiq,ziq = toRegGrid(qu[1],200,200) |
xiq,yiq,ziq = toRegGrid(qu[1],200,200) |
337 |
\end{python} |
\end{python} |
338 |
using the \verb|toRegGrid| function from the cookbook library which we are already using for the contour plot. |
using the \verb|toRegGrid| function from the cookbook library which we are using for the contour plot. |
339 |
At return \verb|ziq[j,i]| is the value of vertical heat flux at point |
At return \verb|ziq[j,i]| is the value of vertical heat flux at point |
340 |
(\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by |
(\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by |
341 |
plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script |
plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script |
358 |
|
|
359 |
\section{Arrow plots in \mpl} |
\section{Arrow plots in \mpl} |
360 |
\sslist{example05c.py and cblib.py} |
\sslist{example05c.py and cblib.py} |
361 |
We would like to visualize the distribution of the flux $-\kappa \nabla T$ over the domain |
We would like to visualise the distribution of the flux $-\kappa \nabla T$ over the domain |
362 |
and produce a plot like shown in \reffig{fig:hr001qumodel}. |
and produce a plot like shown in \reffig{fig:hr001qumodel}. |
363 |
The plot puts together three components. A contour plot of the temperature, |
The plot puts together three components. A contour plot of the temperature, |
364 |
a colored representation of the two sub-domains where colour represents the thermal conductivity |
a colored representation of the two sub-domains where colour represents the thermal conductivity |