Technische Universität MünchenFakultät für MathematikFinite Element Solutionsof Heat Conduction Problemsin Complicated 3D GeometriesUsing the Multigrid MethodDiplomarbeitBastian PentenriederAufgabensteller: Prof. Dr. Christoph ZengerBetreuer:Prof. Dr. Sergey SlavyanovAbgabetermin:29. Juli 2005

Hiermit erkläre ich, dass ich die Diplomarbeit selbständig angefertigt und nur die angegebenen Quellen und Hilfsmittel verwendet habe.München, den 29. Juli 2005.Bastian Pentenrieder

Contents1 Introductory remarks1.1 Aims and motivation of this thesis . . . . . . . . . . . . . . . . . . . . . .1.2 Overview of the individual chapters . . . . . . . . . . . . . . . . . . . . .1.3 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 Modelling2.1 The different modes of heat transfer . . .2.2 Derivation of Fourier’s partial differential2.3 Boundary conditions . . . . . . . . . . .2.4 Formulation of the problem . . . . . . .2.5 A possible generalization . . . . . . . . 04 Results4.1 Graphical user interfaces for data I/O . . . . . . . . . . . .4.2 Calculated heat fluxes and effective thermal conductivities4.3 Visualization and interpretation . . . . . . . . . . . . . . .4.4 Findings from the generalized problem . . . . . . . . . . .6363667282. . . . .equation. . . . . . . . . . . . .3 Solution3.1 Finite element method . . . . . . . . . . . . .3.1.1 The basics . . . . . . . . . . . . . . . .3.1.2 Application to the posed problem . . .3.1.3 Application to the generalized problem3.2 Multigrid method . . . . . . . . . . . . . . . .3.2.1 Motivation . . . . . . . . . . . . . . . .3.2.2 Communication between the grids . . .3.2.3 Hierarchical generating system . . . . .3.2.4 Sketch of a multigrid cycle in peano3d3.2.5 Computation of the diagonal elements3.3 Physical interpretation of the algorithm . . . .3.4 Implementation . . . . . . . . . . . . . . . . .3.4.1 Algorithm . . . . . . . . . . . . . . . .3.4.2 Geometry description . . . . . . . . . .3.4.3 Postprocessing . . . . . . . . . . . . .5

Contents5 Summary and outlook855.1 Achievements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 855.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85A The element stiffness matrix for the tensor generalization689

1 Introductory remarks1.1 Aims and motivation of this thesisPartial differential equations (PDE s) play an important role in a wide range of disciplines. They emerge as the governing equations of problems arising in such different fieldsof study as biology, chemistry, physics and engineering—but also economy and finance.In order to solve them and thus gain a deeper insight into the corresponding issue, oneis in need of powerful algorithms and data structures.A popular and widely used approach to the solution of partial differential equations isthe finite element method (FEM), which, as often in numerical mathematics, reducesthe initial problem to the task of solving a system of linear equations. For this purpose, especially when dealing with a large number of unknowns (e.g. 106 ), classicaldirect solvers turn out to be inappropriate, and more modern iterative schemes like themultigrid method come into play.But not only a sophisticated algorithm is essential. If we fail to arrange the unknownsin the memory in an intelligent way, we will encounter a typical bottleneck related tocontemporary computer architecture. It is caused by the fact that the speed of the CPUhas become faster and faster over the last years, whereas the memory could not keep pacewith that development. The following metaphor might give us an idea of the problem:We imagine a man working on an assembly line. He represents the processor. The piecesdelivered by the line can be thought of as boxes (data packages), each one comprising acertain amount of units (values) that come from the depths of the factory (the memory).Since the man works extremely fast in comparison to the velocity of the assembly line, wehave to ensure that only those units actually needed for the current job (computation)are put into the box. Otherwise, the worker will come across unnecessary parts insteadof required ones, and he will have to wait for the next box to arrive doing nothing inthe meantime (the processor performs idle cycles). Roughly speaking, this means thatan efficient way has to be found how to store things inside the factory, take them out ofthe storage onto the assembly line and afterwards back again.Of course, this conclusion sounds to some extent like a commonplace. However, the truthis that the problem of organizing the data in adequate memory structures in order toavoid the described bottleneck is not trivial. Within the scope of his doctoral thesis,Markus Pögl developed and implemented a method that traverses the computationaldomain on a space-filling curve of Peano type, while it moves the data values througha system of 27 stacks1 ([Pögl 04]). Thereby, it achieves perfect efficiency in terms of the1A stack is a simple structure within the memory which permits only two kinds of data access:pop (to take data from the stack) and push (to put data onto the stack).7

1 Introductory remarksso-called cache hit rate: Whenever the algorithm enters a new cell along the curve, thecorresponding data can be found on top of a priori known stacks. After performing thecomputations, the data is pushed onto (in general) different stacks. When proceeding tothe next cell, the situation there will be the same, i.e. we never have to pop any actuallynot required data from the stacks in order to get access to the elements beneath (fordetails see [Pögl 04]).Since the main objective of the doctoral thesis consisted in demonstrating the cacheefficiency of the new concept, only the Poisson equation as a model problem was implemented in the original code: 4u(x) f (x)u(x) 0 x Ω (0, 1)3 x ΩThe aim of my thesis was to provide the necessary adjustments and extensions to overcome this restriction to a very small class of mathematical problems while maintainingthe cache-efficiency. In the end, the algorithm should be able to solve PDE s of the type div(T (x) u(x)) f (x)with T (x) R3,3 positive definite.The major challenge was the adaptation of the existing multigrid method to coefficientsthat depend on the spatial coordinates, especially in terms of keeping the excellentconvergence properties of the initial code from [Pögl 04].Suitable real-life problems should be modelled and simulated by the resulting implementation, thus illustrating the power of the new concept with respect to practical issues.Interplay with the research work of Prof. SlavyanovIt was in the year 2003 that I got to know Prof. Slavyanov during a student seminarin Saint Petersburg. Together with Prof. Zenger, he supervised the course MathematicalModelling and Numerical Simulation of Problems in Physics. At that time, I myselfparticipated in the program and contributed a talk about the principles of the finiteelement method.One year later, Prof. Slavyanov told us that he was interested in the German HeatInsulation Ordinance and the succeeding EnergieEinsparVerordnung (EnEV). His aimwas to establish an awareness of the necessity of energy-saving among Russian architectsand construction engineers. Connections to the brick-manufacturing company Knauffinally gave rise to the idea of investigating the thermal characteristics of ceramic blocks2with complicated geometry, as they are used in present-day buildings. Choosing thisproblem to be the first test for the algorithm, made particular sense since there alreadyexisted a paper with reference results which was published by V. Grikurov, R. Trepkovand Prof. Slavyanov ([Slavyanov et al. 04]).Therefore, I decided to spend the main period of my diploma thesis on the spot—at theFaculty of Physics of the Saint Petersburg State University.28Throughout this thesis, ceramic block refers to a brick featuring gaps that are filled with air. Whentalking about the brick as a whole, we say ceramic block, in contrast to its two components brickand air.

1.2 Overview of the individual chapters1.2 Overview of the individual chaptersFirst of all, in chapter 2, a brief introduction to heat transfer is given. The differentmechanisms of transportation are presented in some detail to help us understand theboundary conditions of Fourier’s PDE. Besides of that, you will find a derivation for thelatter one, the exact formulation of the problem we are aiming to solve, and also a smalldiscussion of alternative modelling approaches and a possible generalization.After setting up the problem, chapter 3 then deals with its solution. It explains the basicsof the finite element/multigrid method and shows how these techniques can be used forour simulation of heat conduction within ceramic blocks. By means of the knowledgefrom chapter 2, we will be able to recognize that the mathematical algorithm (at leastto some extent) imitates the physical processes inside the material. The last subchapteris dedicated to details of the implementation.Chapter 4 presents the calculated results in numbers and pictures. Both of these arediscussed under various aspects. Particular attention will be paid to the promising convergence properties of the tensor generalization.Finally, chapter 5 gives an outlook on possibilities how to continue the research on thetopic, summarizes the achievements and points out how they can be used in relatedwork.1.3 AcknowledgementsA lot of people contributed to the success of this project. In the first instance, I amdeeply indebted to Prof. Dr. Christoph Zenger and Prof. Dr. Sergey Slavyanov for theinspiring and patient supervision of my thesis and their great help in organizing the stayin Saint Petersburg.Furthermore, I would like to thank Markus Pögl, Frank Günther and Andreas Krahnkefor all the explanations during the initial phase of my work. Without them, the processof understanding the starting-point implementation of the algorithm would have takenmuch longer. The useful hints of Dr. Miriam Mehl concerning the multigrid methoddeserve a special mention. Wolfgang Herder and Ellen Maier were so kind to proof-readthe final version of this thesis.The pleasant atmosphere among the fellow graduands and scientists at the chair of Prof.Zenger and the amiability of Prof. Slavyanov and all the Russian people at his facultywith whom I became acquainted have been a very procreative environment for the wholetime.Last–but not least–the share of my parents and my girlfriend must not be forgotten.They always encouraged me in my plans. I am very grateful that my mother and myfather made it possible for me to attend university, and I want them to know that Idon’t take for granted all the support they have given me through the years.9

1 Introductory remarks10

2 ModellingIn scientific computing, the examination of a certain problem is always a three-stageprocess (figure 2.1). Foremost, we need to think about how to model the task in thelanguage of mathematics. Then, the next step is to discretize the obtained (continuous)equations and solve them by an appropriate numerical method. Finally, we have topresent the results in an adequate way. This may be done by giving some key values orby applying visualization techniques to the computed datasets.Figure 2.1: Scientific computing—a three-stage processAccording to this scheme, we start with introducing some fundamental laws of heattransfer that will help us to translate the heat conduction problem within ceramic blocksinto mathematical equations. For profound studies on this branch of engineering, theinterested reader is recommended the definitive textbooks [Incropera/DeWitt 02] and[Baehr/Stephan 03].2.1 The different modes of heat transferBy definition, heat is the energy that flows from the higher level of temperature to thelower (without any work being performed), whenever there exists a temperature diffe-11

2 Modellingrence or gradient respectively. The equations which quantify the amount of transferredheat all fit into the patternflow transport coefficient potential gradient,where the potential gradient represents a derivative or some difference expression. As flowvalue, we encounter heat flux q̇ [W/m2 ]: energy per time and area or heat transfer rate Q̇ [W ]: energy per time (passing through a fixed reference area).The transport coefficient depends on the particular mode of transfer . . .ConductionHeat conduction is the diffusive transport of thermal energy. In liquids and gases, it iscaused by the interaction of moving atoms and molecules, in solids by lattice oscillationsand in electroconductive material additionally by unbound electrons.For instance, consider a swimming pool that is warm on the left and cold on the righthand side (cf. figure 2.2). We suppose that the water is completely motionless. On aFigure 2.2: Brownian motion of H2 O molecules (red: warm, blue: cold)microscopic level, neighboring molecules are in mutual interaction and exchange kineticenergy—more or less chaotically. But if we change the scale of our observation to themacroscopic level, it is evident that, in total, more kinetic energy will be transferredfrom the left to the right than vice versa. The system tends towards an equilibriumstate with the same temperature everywhere. During this levelling diffusion process, atransport of thermal energy occurs.If we know the internal temperature distribution for a given moment in time, we cancalculate the heat flu