PDEModelica

(1)

Linköping Studies in Science and Technology Dissertation No. 1016

A High-Level Language for Modeling with

Partial Differential Equations

by

Levon Saldamli

Department of Computer and Information Science Linköpings universitet

(2)

(3)

“. . . to boldly go where

no man has gone before.”

(4)

(5)

Abstract

This thesis describes work on a new high-level mathematical modeling language and framework called PDEModelica for modeling with partial differential equa-tions. It is an extension to the current Modelica modeling language for object-oriented, equation-based modeling based on differential and algebraic equations. The language extensions and the framework presented in this thesis are consistent with the concepts of Modelica while adding support for partial differential equa-tions and space-distributed variables called fields.

The specification of a partial differential equation problem consists of three parts: 1) the description of the definition domain, i.e., the geometric region where the equations are defined, 2) the initial and boundary conditions, and 3) the ac-tual equations. The known and unknown distributed variables in the equation are represented by field variables in PDEModelica. Domains are defined by a geo-metric description of their boundaries. Equations may use the Modelica derivative operator extended with support for partial derivatives, or vector differential opera-tors such as divergence and gradient, which can be defined for general curvilinear coordinates based on coordinate system definitions.

The PDEModelica system also allows the partial differential equation models to be defined using a coefficient-based approach, where PDE models from a library are instantiated with different parameter values. Such a library contains both con-tinuous and discrete representations of the PDE model. The user can instantiate the continuous parts and define the parameters, and the discrete parts containing the equations are automatically instantiated and used to solve the PDE problem numerically.

Compared to most earlier work in the area of mathematical modeling languages supporting PDEs, this work provides a modern object-oriented component-based approach to modeling with PDEs, including general support for hierarchical mod-eling, and for general, complex geometries. It is possible to separate the geometry definition from the model definition, which allows geometries to be defined sep-arately, collected into libraries, and reused in new models. It is also possible to separate the analytical continuous model description from the chosen discretization and numerical solution methods. This allows the model description to be reused, independent of different numerical solution approaches.

(6)

work affords a clearer abstraction and defines a new type of variable. Arrays of such field variables can be defined in the same way as arrays of regular, scalar variables. The PDEModelica language supports a clear, mathematical syntax that can be used both for equations referring to fields and explicit domain specifica-tions, used for example to specify boundary conditions. Hierarchical modeling and decomposition is integrated with a general connection concept, which allows connections between ODE/DAE and PDE based models.

The implementation of a Modelica library needed for PDEModelica and a pro-totype implementation of field variables are also described in the thesis. The PDEModelica library contains internal and external solver implementations, and uses external software for mesh generation, requisite for numerical solution of the PDEs. Finally, some examples modeled with PDEModelica and solved using these implementations are presented.

(7)

Acknowledgments

First of all, I would like to thank my supervisor Peter Fritzson, for the inspiration and the motivation I needed to finish this thesis, and for commenting both my work and my language. I would also like to thank my co-supervisor Bernhard Bach-mann, for help with parts of the implementation in this work, and for his valuable comments on the contents of this thesis. Thanks also to Hansjürg Wiesmann for help with the implementation and for interesting discussions.

Other colleagues at PELAB I would like to thank, who have contributed to the OpenModelica platform, where I could experiment with language extensions, are Peter Aronsson, Kaj Nyström, Adrian Pop, Peter Bunus, Håkan Lundvall, and David Broman. Thanks to the rest of the PELAB members for interesting dis-cussions about both work and other more or less relevant subjects, especially to Jens Gustavsson, Andreas Borg, and John Wilander for our very interesting work together on research methodology and computer science. Special thanks to An-drzej Bednarski for being both a colleague and a friend. Many thanks to Bodil Mattsson-Kihlström for her support in many ways, and to Lillemor Wallgren for her administrative help.

Finally, I would like to thank my parents for enhancing my motivation, and my friends, for their company and support that helped me finish this work.

This work has been supported by the Swedish Foundation for Strategic Research in the ECSEL graduate school and the VISIMOD project, the European Commis-sion in the RealSim project, and by Vinnova in the VISP project.

(8)

(9)

Glossary

class A user-defined, reusable component type in Modelica, containing variables,

instances of other classes, algorithms and equations.

restricted class A specialized Modelica class, with some restrictions on the

contents it may have, for example occurence of equations.

extends Term used for subclassing of Modelica classes, for defining inheritance

hierarchies of classes.

modification In Modelica, a mechanism to change values of parameters, types

of declared elements etc., while inheriting or instantiating a class.

connector A specialized Modelica class to represent interfaces between

compo-nents, that can be used with the connect operator to automatically generate equations.

domain The definition domain of the independent variables. For time-dependent

functions, the domain is the time interval where the variable is defined. For two or three-dimensional space-dependent functions, the domain is the two and three-dimensional geometric regions where the function is defined, re-spectively.

boundary The bounding limit of the definition domain.

field A function of independent variables, for example a space-dependent, two

dimensional function f(x,y).

field value Value of a field at a point in the definition domain.

partial derivative Derivative of an expression or function with respect to one of

the independent variables, e.g., one of the space coordinates.

PDE Partial differential equation, i.e., an equation containing partial derivatives

of functions.

ODE Ordinary differential equation. containing functions of one independent

variable and derivatives of functions and expressions with respect to that independent variable.

(14)

DAE Differential and algebraic equations. An algebraic system of equations

con-taining ordinary derivatives.

boundary condition Conditions to determine a unique solution to a PDE,

de-fined on the boundary of the definition domain where the PDE is solved. The PDE is usually valid only in the interior of the domain, while the boundary conditions give the values of the unknown or its derivatives on the boundary of the domain.

initial condition In time-dependent equations, the starting values of the unknown

variables, needed to determine a unique solution to the time-dependent prob-lem.

Dirichlet A type of boundary condition, specifying the value of the unknown

vari-able on the boundary.

Newman A type of boundary condition, specifying the value of the derivative of

the unknown variable on the boundary.

Robin A type of boundary condition, consisting of an equation containing both

the unknown variable and its derivative. Also called mixed boundary condi-tion.

elliptic Usually steady-state PDEs containing no time derivatives, e.g., the

Pois-son equation.

parabolic Usually time-dependent PDEs containing first-order time derivatives,

e.g., the time-dependent heat transfer equation.

hyperbolic Usually time-dependent PDEs containing second-order time

deriva-tives, e.g., the wave equation.

stationary (steady-state) The time-independent equation, corresponding to the

time-dependent equation with the time elapsed enough long so that the solu-tion becomes time-independent, e.g., the transients have vanished.

finite difference An approximation of the derivative using the derivative

defini-tion and the values of the unknown at discrete grid points and the distance between the grid points.

central difference A finite difference method using discrete values from both

previous and next discrete point.

forward difference A finite difference method using discrete values from the

(15)

Contents

backward difference A finite difference method using discrete values from the

previous discrete point.

finite difference method A method to solve PDE problems using finite

differ-ences to approximate partial derivatives and solve the resulting equation sys-tems.

finite element method (FEM) A method to solve PDE problems by

subdivid-ing the domain in for example triangles and expresssubdivid-ing the unknown in terms of known basis functions and solving the resulting equation systems.

method of lines A method to solve time-dependent PDE problems by replacing

the time- or space-dependent parts of the equations with discrete formula-tions, resulting in a only time- or only space-dependent problem. In this work, spatial discretization is used, generating DAEs from time-dependent PDEs.

(16)

(17)

Chapter 1 Introduction

“Its five year mission: To explore strange new worlds, to seek out new life and new civilizations, to boldly go where no man has gone before.” Captain Kirk, Starship Enterprise

In mathematical modeling of physical systems, the behavior of a system is typi-cally modeled using differential equations or a system of differential and algebraic equations, i.e. algebraic equations containing derivatives.

For large scale models, some properties of the physical system that is modeled are often approximated; for example when controlling the temperature of a fluid in a container, the temperature can be assumed to be the same everywhere inside the container. This way, the space dependency is eliminated, resulting in models con-taining only ordinary differential equations (ODEs), or differential and algebraic equations (DAEs), i.e. equations where all derivatives are with respect to time only. This greatly simplifies the solution of the equations and the simulation.

However, sometimes a more detailed model with explicit space dependency is needed for the whole or some part of the system. Such space-dependent models are used to study the behavior in more detail, and contain partial derivatives, i.e., derivatives with respect to space variables like x, y and z in a Cartesian coordinate system. Differential equations containing partial derivatives are called partial dif-ferential equations (PDEs), and are used in many fields to model physical behavior, for example structural mechanics, computational fluid dynamics, and electrostat-ics.

There are many tools for modeling and simulation with ordinary differential equations as well as with partial differential equations. Most are however spe-cialized for certain application domains or special kinds of models.

(18)

An effort to define a general, application domain and tool-independent model-ing language is made by Modelica Association, with the Modelica language as a result. Modelica®_{[1, 2] is a standard language designed for component-based} modeling, where previously defined models can be reused as components in new models. Object oriented constructs, equation-based and declarative modeling sup-port, as well as a connection concept make this possible. Currently, Modelica only supports models with ordinary differential equations, algebraic equations, discrete equations, hybrid (mixed time-discrete and time-continuous) equations but not par-tial differenpar-tial equations. Several Modelica implementations and simulation tools exist [3, 4].

Thanks to the progress in computer performance, the use of detailed models con-taining partial differential equations is becoming increasingly common. The aim of this work is to define a language based on Modelica, in order to handle partial differential equations in models in addition to differential and algebraic equations. This language is called PDEModelica, and this thesis discusses the initial language extensions needed for the desired support.

1.1 PDE-based Model Example

Consider a simplified model of heat distribution in a room where only the distri-bution in thexandydirections is studied and the temperature in thezdirection is assumed to be constant. A heater is installed on one of the walls and there is a window on another wall. The domain of this problem, i.e., the region (including geometry) on which it is defined, can be seen in Figure 1.1. A heater is represented by assigning constant temperature boundary condition to the middle part of the upper wall. The window is modeled by non-zero heat flow through the left wall. The heat flow is proportional to the temperature difference between the inside tem-perature and the outside temtem-perature. The rest of the walls are insulated, i.e., no heat flow occurs through these walls. The right wall is curved in order to illustrate modeling with complex geometries.

Heat transfer by conduction in two dimensions can be modeled using Poisson’s equation:

−(∂2T(x, y) ∂x2 +

∂2T(x, y)

∂y2 ) = g(x, y) on Ω The possible boundary conditions are the following:

• Dirichlet1boundary condition for the heated wall: T(x, y) = h₁(x, y) on ∂Ω

(19)

1.1 PDE-based Model Example Heater

Window

left

bottom

top3 top2 top1

right

x y

h

w

Figure 1.1: Example of time-dependent heat transfer by conduction.

• Neumann2_{boundary condition for insulated walls:} ∂T(x, y)

∂n = h2(x, y) on ∂Ω

• Robin3(also called mixed) boundary condition for the window which is not perfectly insulated:

∂T(x, y)

∂n + qT (x, y) = h3(x, y) on ∂Ω

If heat distribution over time is studied, a time derivative is added to the Laplace equation which gives the heat diffusion equation:

ρC∂T(x, y, t) ∂t − k( ∂2T(x, y, t) ∂x2 + ∂2T(x, y, t) ∂y2 ) = g(x, y, t) on Ω with ρ, C and k being material constants, and a given initial condition:

T(x, y, 0) = T₀(x, y) on Ω

This problem can be modeled in PDEModelica with the geometry and the model description shown in Figure 1.2. The domain geometry is defined by describing

2_{The value of the outward normal derivative of the unknown variable is known on the boundary.}

The outward normal derivative is the space derivative in the outward normal direction. See also Section 2.2.4.

3_{A linear combination of the unknown variable and its outward normal derivative is known on the}

(20)

its boundary in a specific direction in order to find out on what side of the bound-ary the actual domain is. The boundbound-ary is built up of several sections defined as named components in the domain description. Hence, when defining the PDE model, boundary conditions can be assigned to each section using the name of each boundary component. The definition of the boundary conditions and the PDE are not shown in this example; they can be defined as described in Chapter 4. The solution to a similar example can be found in Section 6.2.1.

1.1.1 Connection to ODE/DAE models

A system with active temperature control can be modeled using a time-dependent model as shown in the previous section. A temperature sensor can be approximated by reading the computed temperature at some point in the domain, and the temper-ature can be used as input to a heat controller. The model overview is illustrated in Figure 1.3, and the PDEModelica code is listed in Figure 1.4.

The controller is a proportional, integrating controller (PI-controller), which contains an ordinary differential equation:

Terror= Tsensor− Tgoal

Theater(t) = kpTerror+ ki

Terrordt

The resulting problem consists of a PDE part modeling temperature distribution, and an ODE part modeling the controller. This can be compared to a simplified, lumped system where the temperature is assumed to be the same at all points, i.e. where the temperature distribution in the domain is instant. Temperature dis-tribution in such a system is modeled using an ODE-based model, resulting in a complete system with only ODEs. However, if a more realistic temperature model is used, involving PDEs, where spatial variation in temperature is studied as well, coupled ODE and PDE models are needed. Spatial dependency can then be re-moved by spatial discretization as a first step, using the method of lines, resulting in an equation system of ODEs also for the spatially distributed temperature model. Thus, the final discretized system will only contain ODEs which is solved using ex-isting ODE/DAE solvers.

This example is simplified by the fact that the sensor is replaced with a simple field access to read the temperature value. The connection between the controller and the PDE model is accomplished manually by writing down the equations. In the future, the connection between PDE and ODE models can be more complex, using connectors and field reduction, see Section 4.3.4 and Section 9.3.2.

(21)

1.1 PDE-based Model Example c l a s s C u r v e d R e c t a n g u l a r " G eom etr y " e x t e n d s C a r t e s i a n 2 D ( b o u n d a r y ={ bottom , r i g h t , t o p 1 , t o p 2 , t o p 3 , l e f t } ) ; parameter P o i n t p0 ; parameter R e a l w ; parameter R e a l h ; parameter R e a l cw ; parameter R e a l ch =h ; parameter P o i n t c c =p0 + {w, h / 2 } ; Line2D b o t t o m ( p1=p0 , p2=p0 + {w , 0 } ) ; Line2D t o p 1 ( p1 =p0 + {w, h } , p2= p0 + {2*w/ 3 , h } ) ; Line2D t o p 2 ( p1 =p0 + {2*w/ 3 , h } , p2=p0 + {w/ 3 , h } ) ; Line2D t o p 3 ( p1 =p0 + {w / 3 , h } , p2=p0 + { 0 , h } ) ; Line2D l e f t ( p1 =p0 + { 0 , h } , p2= p0 ) ; B e z i e r 2 D r i g h t ( n =8 , p= f i l l ( cc , 8 ) + { { 0 . 0 ,−0.5} , {0.0 , −0.2} , { 0 . 0 , 0 . 0 } , {−0.85 , −0.85} , { −0.85 ,0.85} , { 0 . 0 , 0 . 0 } , { 0 . 0 , 0 . 2 } , { 0 . 0 , 0 . 5 } } *{{cw , 0 } , { 0 , ch } } ) ; end C u r v e d R e c t a n g u l a r ; model D i f f u s i o n " H eat d i f f u s i o n w i t h h e a t e r " . . . imp ort M o d e l i c a . S I U n i t s . T e m p e r a t u r e ; parameter T e m p e r a t u r e h e a t e r V a l u e ; parameter R e a l h_3 ; parameter R e a l q ; C u r v e d R e c t a n g u l a r omega ; D i f f u s i o n 2 D pde ( dom ain =omega ) ;

e q u a t i o n d er ( pde . u , n ) = 0 on omega . b o t t o m ; d er ( pde . u , n ) = 0 on omega . r i g h t ; d er ( pde . u , n ) = 0 on omega . t o p 1 ; pde . u = h e a t e r V a l u e on omega . t o p 2 ; d er ( pde . u , n ) = 0 on omega . t o p 3 ;

d er ( pde . u , n ) + q* pde . u = h_3 on omega . l e f t ; end D i f f u s i o n ;

Figure 1.2: The PDEModelica code to define the heat diffusion problem in Sec-tion 1.1. Predefined PDE modelDiffusion2Dand detailed parameter declarations have been left out. The syntax is explained in Chapter 4.

(22)

Heater

Glass layer

PI-Controller

Sensor

Figure 1.3: A coupled PDE and ODE model. Time dependent heat transfer model with a controller.

1.2 Research Method

The research method used in this work can be called explorative design [5]. The approach starts with a design hypothesis which is believed to be a possible solution to the problem. A prototype implementation of this design is developed. During implementation, flaws in the design are discovered and new knowledge about the problem is acquired. The implementation is modified to handle these issues. Thus, the system that is implemented differs from that in the original design. On the other hand, the implementation also differs from the one that would have been implemented from start using the new knowledge. For this reason, at some point, the prototype is discarded and a new design and prototype are developed. See Figure 1.5 for an illustration of the work flow when using this method.

Thus, important parts of the research results are the original design, the actual implemented prototype, and a proposed, improved design. In this thesis, we present the previous iterations of the design and implementation in Chapter 7. The final design of PDEModelica in this thesis work is presented in Chapter 4 and the final implementation is described in Chapter 3 and Chapter 5.

(23)

1.2 Research Method c l a s s C u r v e d R e c t a n g u l a r " G eom etr y " e x t e n d s C a r t e s i a n 2 D ( b o u n d a r y ={ bottom , r i g h t , t o p 1 , t o p 2 , t o p 3 , l e f t } ) ; parameter P o i n t p0 ; parameter R e a l w ; parameter R e a l h ; parameter R e a l cw ; parameter R e a l ch =h ; parameter P o i n t c c =p0 + {w, h / 2 } ; Line2D b o t t o m ( p1=p0 , p2=p0 + {w , 0 } ) ; Line2D t o p 1 ( p1 =p0 + {w, h } , p2= p0 + {2*w/ 3 , h } ) ; Line2D t o p 2 ( p1 =p0 + {2*w/ 3 , h } , p2=p0 + {w/ 3 , h } ) ; Line2D t o p 3 ( p1 =p0 + {w / 3 , h } , p2=p0 + { 0 , h } ) ; Line2D l e f t ( p1 =p0 + { 0 , h } , p2= p0 ) ; B e z i e r 2 D r i g h t ( n =8 , p= f i l l ( cc , 8 ) + { { 0 . 0 ,−0.5} , {0.0 , −0.2} , { 0 . 0 , 0 . 0 } , {−0.85 , −0.85} , { −0.85 ,0.85} , { 0 . 0 , 0 . 0 } , { 0 . 0 , 0 . 2 } , { 0 . 0 , 0 . 5 } } *{{cw , 0 } , { 0 , ch } } ) ; end C u r v e d R e c t a n g u l a r ; model C o n t r o l l e d D i f f u s i o n " H eat d i f f u s i o n w i t h c o n t r o l l e d h e a t e r " . . . imp ort M o d e l i c a . S I U n i t s . T e m p e r a t u r e ; parameter T e m p e r a t u r e g o a l V a l u e ; parameter R e a l h_3 ; parameter R e a l q ; / / Use P I r e g u l a t o r f r o m S t a n d a r d L i b r a r y M o d e l i c a . B l o c k s . C o n t i n u o u s . P I r e g u l a t o r ; C u r v e d R e c t a n g u l a r omega ;

D i f f u s i o n 2 D pde ( dom ain =omega ) ; T e m p e r a t u r e h e a t e r V a l u e ; T e m p e r a t u r e s e n s o r V a l u e ; e q u a t i o n d er ( pde . u , n ) = 0 on omega . b o t t o m ; d er ( pde . u , n ) = 0 on omega . r i g h t ; d er ( pde . u , n ) = 0 on omega . t o p 1 ; pde . u = h e a t e r V a l u e on omega . t o p 2 ; d er ( pde . u , n ) = 0 on omega . t o p 3 ;

d er ( pde . u , n ) + q* pde . u = h_3 on omega . l e f t ;

h e a t e r V a l u e = r e g u l a t o r . y [ 1 ] ;

r e g u l a t o r . i n P o r t . s i g n a l [ 1 ] = g o a l V a l u e − s en s o r V a l u e ; s e n s o r V a l u e = pde . u ( 2 . 6 1 , 2 . 3 3 ) ;

end C o n t r o l l e d D i f f u s i o n ;

Figure 1.4: The PDEModelica code to define the heat diffusion problem with a con-trolled heater in Section 1.1.1. Predefined PDE modelDiffusion2D

(24)

Prototype implementation Design Hyptohesis Modified Design Hyptohesis Improved Design Hyptohesis New knowledge Iterations

Figure 1.5: Explorative design.

1.3 Contributions

The result of this work is two-fold. First, a design and implementation of a Model-ica library for PDEs is presented, which illustrates the concepts needed to specify and simulate PDE-based models. Second, a design and a partial implementation of language extensions for the Modelica language is presented, which integrates the concepts from the PDE library into the Modelica language. The proposed lan-guage design for an object-oriented modeling lanlan-guage, PDEModelica, supports the following:

• Object-oriented, component-based modeling with PDEs.

• General, complex geometry description with lines, polygons, parametric curves and a combination of these.

• Component-based geometry definition, i.e., separation of geometry defini-tion from model definidefini-tion.

• Separation of analytical, continuous model description from discretization and numerical solution.

• Field concept for general declaration of spatially distributed variables inte-grated in the language.

(25)

1.4 Overview of the Thesis • Clear, mathematical syntax for equations containing fields and explicit

do-main specification.

• Hierarchical modeling and decomposition.

• Modeling of combined ODE/DAE and PDE problems.

We believe that the combination of these aspects in a single modeling language is rather new. See also the discussion and comparison to other related systems in Section 8.3.

A prototype translator and solver environment is set up as well, with a basic PDE solver interface for adding new solvers and mesh generators.

Some of the work discussed in this thesis has been published previously in the following publications:

1. L. Saldamli and P. Fritzson. Object-Oriented Modeling with Partial Differ-ential Equations. In Proceedings of the Modelica Workshop 2000, Lund, Sweden, October 2000 [6].

2. L. Saldamli and P. Fritzson. A Modelica-Based Language for Object-Orient-ed Modeling with Partial Differential Equations. In ProceObject-Orient-edings of the 4th International EUROSIM Congress, Delft, The Netherlands, June 2001 [7]. 3. L. Saldamli and P. Fritzson. Domains and Partial Differential Equations in

Modelica. In Proceedings of the 42nd SIMS Conference, Porsgrunn, Norway, October 2001 [8].

4. L. Saldamli, P. Fritzson and B. Bachmann. Extending Modelica for Partial Differential Equations. In Proceedings of the 2nd International Modelica Conference, Oberpfaffenhofen, Germany, March 2002 [9].

5. L. Saldamli and P. Fritzson. Field Type and Field Constructor in Model-ica. In Proceedings of SIMS 2004, the 45th Conference on Simulation and Modelling, Copenhagen, Denmark, September 2004 [10].

6. L. Saldamli, B. Bachmann, H Wiesmann and P. Fritzson. A Framework for Describing and Solving PDE Models in Modelica. In Proceedings of the 4th International Modelica Conference, Hamburg, Germany, March 2005 [11]. Section 4.3 contains material not previously published.

1.4 Overview of the Thesis

This thesis is organized as follows. Chapter 2 contains background information rel-evant for the thesis: a short overview of the Modelica language, a basic introduction

(26)

to partial differential equations and some existing numerical solution methods for partial differential equations.

Chapter 3 describes a PDE library consisting of Modelica packages written us-ing standard Modelica. The PDE library can be used to describe and solve PDE models on general domains using the finite element method and the method of lines discretization.

Chapter 4 presents proposed extensions to the Modelica language for definition of space-distributed models, e.g. models containing fields and partial derivatives.

Chapter 5 contains the implementation details. The PDE library, prototypes for language extensions and the interface to external solvers are described here.

Chapter 6 contains examples to demonstrate the use of the PDE library and the language extensions.

Chapter 7 discusses the previous design iterations that were implemented and tested.

Chapter 8 presents an overview of related work. Different low- and high-level tools with varying modeling language support are summarized here.

Finally, some conclusions and possible future work directions are discussed in Chapter 9.

(27)

Chapter 2 Background

“Logic is the beginning of wisdom; not the end.”

Mr. Spock

The basic concepts of Modelica are presented in the first part of this chapter. In the second part, the topic of partial differential equations (PDEs) and classifica-tion of PDEs are briefly presented, and finally an overview of different numerical solution methods is given.

2.1 Modelica

Modelica [2, 12] is a modeling language for equation-based, object-oriented mod-eling and simulation of physical systems. Using object-oriented concepts, it allows hierarchical, component-based modeling which in turn makes reuse of existing models possible. The general modeling concepts in Modelica allow it to be used in different application domains and in multi-domain modeling, for example when defining a combined electrical and mechanical model. There is a free Modelica Standard Library [13] and other free Modelica libraries with packages of existing models from different domains that can be used as components in specific models. Acausal modeling, i.e., modeling using equations - not assignment statements, al-lows single models to be used in different data flow contexts since variables are not explicitly declared as input or output. This causality information is rather derived from the context where the model is used. In his PhD thesis [14], H. Tummescheit discusses design and implementation of modeling libraries using Modelica.

(28)

Models in Modelica are built hierarchically from sub-models that are defined separately, which in turn can be defined in the same way. Hence, a model can con-tain one or more instances of other models with a different set of parameter values for each instance, and a set of connections between these components/instances. Additionally, each model can have variables of built-in types, and equations that define the relationship between these variables as well as between variables in sub-models.

2.1.1 Classes

The basic structural unit in Modelica is class. A class can contain other classes, variable declarations, equations, and algorithms. An example of a class can be seen in Figure 2.1. Other keywords can be used to denote classes, restricted classes that are special cases of classes with some restrictions. These are record, type, connector, model, block, package and function. Different restricted classes and restrictions that apply to them can be seen in Table 2.1.

model MyModel " S h o r t d e s c r i p t i o n o f t h i s m odel "

R e a l x , y ;

e q u a t i o n

x + y = 5 ; x + 2 y = 1 1 ;

end MyModel ;

Figure 2.1: Example of a Modelica class.

2.1.2 Variable Declarations

Variable declarations consist of a type and a variable name. The type can be one of the primitive types or another class name. The primitive types are Real, In-teger, String, Boolean and Enumeration. Each declaration can also have type modifiers, such as the variability type prefix constant, discrete, or pa-rameter, or a causality prefix, input or output, used in function and block classes to declare the direction of a variable.

Variables with the constant or parameter prefix keep their value constant during a simulation, i.e., they are constants. The difference between parameters and constants is that values of parameters can be changed before each simula-tion without needing to recompile the Modelica code. An excepsimula-tion is structural parameters , which require recompilation of the model because they change the model structure. A parameter used as number of components in an array declara-tion is an example of a structural parameter. Parameters also appear in tool menus

(29)

2.1 Modelica Table 2.1: Restricted classes in Modelica.

Restricted class Restrictions

record No equations are allowed.

type May only be used to extend primitive types, enumerations, or records.

connector No equations are allowed.

model Model instances may not be used in connections. block Fixed causality. Each variable must be declared

as input or output.

package May only contain classes and constants. Is allowed to import from.

function Similar to block. No equations, at most one algorithm section. May be called using positional arguments. May be recursive.

to allow the user to change them before simulation. Variables of type Real de-clared without the constant or parameter modifiers are implicitly time-dependent, i.e. functions of time. Time derivatives of variables can be represented using the der()operator. Using the time derivatives, ordinary differential equations can be expressed, and full systems of differential and algebraic equations (DAEs) can be specified.

2.1.3 Subtyping and Inheritance

Modelica models can be defined using inheritance. The extends clause is used together with a class name to inherit that class. Multiple inheritance is allowed us-ing several extends clauses, one for each inherited class. Inheritance is equiva-lent to inserting the contents of the inherited class at the place where the extends clause resides.

Type equivalence in Modelica is defined as follows: Two types T and U are equivalent if and only if

• they denote the same primitive type (see Section 2.1.2), or

• T and U are classes containing the same elements (according to their names) and the elements types are equivalent.

Subtypes in Modelica are defined independently from the inheritance mecha-nism. A class C is a subtype of a class S if

(30)

• both of the following statements hold:

– every public element of S also exists in C (according to its name). – the type of each such element in C is a subtype of the type of the

corre-sponding element in S.

If a class C is a subtype of a class S then S is called the supertype of C. Subtypes and supertypes do not necessarily need to be in an inheritance hierarchy, but a class that is inherited from, a base class, is a supertype of a class that inherits from it, a derived class. An example is shown in Figure 2.2.

p a r t i a l model O n e P o r t P i n p , n ; V o l t a g e v " V o l t a g e d r o p " ; e q u a t i o n v = p . v − n . v ; p . i + n . i = 0 ; end O n e P o r t ; model R e s i s t o r ; e x t e n d s O n e P o r t ; parameter R e a l R ( u n i t = "Ohm" ) " R e s i s t a n c e " ; e q u a t i o n v = p . i * R ; end R e s i s t o r ; model T e m p R e s i s t o r " T e m p e r a t u r e d e p e n d e n t r e s i s t o r " e x t e n d s O n e P o r t ; parameter R e a l R ( u n i t = "Ohm" ) " R e s i s t a n c e a t r e f e r e n c e t e m p e r a t u r e " ;

parameter R e a l RT ( u n i t = "Ohm / degC " ) = 0

" T e m p e r a t u r e d e p e n d e n t r e s i s t a n c e " ;

parameter R e a l T r e f ( u n i t = " degC " ) = 20 " R e f e r e n c e temp . " ;

R e a l Temp = 20 " A c t u a l t e m p e r a t u r e " ;

e q u a t i o n

v = p . i * (R + RT * ( Temp − T r e f ) ) ;

end T e m p R e s i s t o r ;

Figure 2.2: Subtyping in Modelica. TempResistor does not inherit from Re-sistor, but it is a subtype of Resistor. They both inherit from OnePort.

The TempResistor class cannot extend the Resistor class because it has a different equation in the equation section, but it is still a subtype of Resistor because it contains all the elements in Resistor, and additional elements. Both classes inherit from the OnePort1_{class, though.}

(31)

2.1 Modelica

2.1.4 Modifications

A modification in Modelica is a short-hand Modelica syntax for expressing class specialization. For example, the modificationcompa(pa=1)in ModelB in Figure 2.3 expresses that the equationpa=0in ModelA should be replaced by the equation

pa=1to create a temporary, specialized class derived from ModelA and used for creation of the compa instance.

Since models are built up hierarchically, modifications can be overridden, in which case the topmost (outermost) modification is applied. A modification can also be hierarchical, in order to modify a lower level parameter directly. Each modification consists of a component reference and an expression that is evaluated in the context where the declaration resides. Example of modifications can be seen in Figure 2.3. Declarations preceded by the final keyword are final elements and cannot be modified using modifications. The final keyword can also precede modifications, in which case it prevents further modifications of that variable in outer modifications. model ModelA parameter R e a l pa = 0 ; / / d e f a u l t v a l u e end ModelA ; model ModelB ModelA compa ( pa = 1 ) ; / / m o d i f i c a t i o n parameter R e a l pb = 0 ; end ModelB ; model MyModel

parameter R e a l mypa =3 , mypb = 4 ;

/ / h i e r a r c h i c a l m o d i f i c a t i o n

ModelB compb ( pb =mypb , compa ( pa =mypa ) ) ;

end MyModel ;

Figure 2.3: Modifications in Modelica. The resulting value of the variable compb-.compa.pais mypa, i.e., 3.

Modifications can also be applied directly when extending a class, in the ex-tendsclause. For example:

model MyNewModel

e x t e n d s MyModel ( mypb = 5 ) ;

When instantiating components of MyNewModel the default value of mypb will be 5.

(32)

2.1.5 Equations and Algorithms

The keywords equation and algorithm denote equation and algorithm sec-tions, respectively. There can be several such sections of each kind in a class. During compilation, all equation sections are merged into a single equations sec-tion.

Equation sections can contain standard equation clauses, connect equations, con-ditional equations, for equations and when equations. Standard equations consist of two expressions and the equality operator denoted by =. The assignment oper-ator := is not allowed in equation sections, whereas the equality operoper-ator = is not allowed in algorithm sections. Besides the common operators like arithmetic operators and function calls, Modelica has if-expressions, which can be used in situations like:

e q u a t i o n

x = i f _{( z > 0 ) 2* z e l s e (−2* z ) ;}

Connect equations consist of the connect() operator with two arguments that are references to connectors. Connectors are instances of connector restricted classes and usually contain two kinds of variables, flow (e.g. current) and non-flow (e.g. potential) variables, the former being denoted by the flow keyword. Connect equations are translated into standard equations during compilation, where for each connection the non-flow variables of the involved connectors are set equal, and the sum of all the flow variables is set to be equal to zero. Figure 2.4 illustrates the use of connectors and connections, and the resulting equations.

Conditional equations can be written using the if-elseif-else construct with equations in the body. for equations are useful for repetitive equation struc-tures involving arrays, whereas when equations and when statements are used for discrete event simulation.

Algorithm section can contain assignments, conditional statements, for-loops, while-loops and when-statements. The difference from equation sections is that algorithm sections contain assignments instead of equations, and the contents of the algorithms are executed sequentially. An assignment consists of a variable reference as the left-hand side, the assignment operator := and an expression as the right-hand side.

2.1.6 Replaceable Elements

Elements in classes in Modelica can be declared as replaceable, in which case they can be replaced by another type of element that is type compatible when the containing class is extended or instantiated. Generic classes with type parameters is supported in Modelica through the replaceable mechanism.

The type restriction of a redeclaration is that the new declared type must be a subtype of the original declared type, a constraining type.

(33)

2.1 Modelica

B

A

p

C

(a) Visual view of a connector and connection example.

c o n n e c t o r P i n " E l e c t r i c a l p i n " f l o w C u r r e n t i ; V o l t a g e v ; end P i n ; model OnePin P i n p ; end OnePin ; model T e s t OnePin a , b , c ; e q u a t i o n c o n n e c t ( a . p , b . p ) ; c o n n e c t ( b . p , c . p ) ; end T e s t ;

(b) Textual view of a connector and connection example.

e q u a t i o n

a . p . i + b . p . i + c . p . i = 0 ; a . p . v = b . p . v ;

b . p . v = c . p . v ;

(c) Resulting equations

Figure 2.4: Example of connectors and connections and the resulting simple equa-tions after compilation.

(34)

2.2 Differential Equations

Differential equations are commonly used in mathematical models of physical phe-nomena that are distributed in time and/or space. Examples include heat transfer, fluid flow, structural mechanics, wave propagation, etc. Such mathematical models relate certain variables, called dependent variables, and their derivatives to other variables, called independent variables. For example the time variable or space coordinates occur commonly as independent variables.

2.2.1 Fields

A field is a mapping from a domain to scalar or vector values. The geometrical region, i.e., the domain where the field is defined, is called the fields definition domain. A field can be seen as a continuous variable distributed in spatial coordi-nates. Example of a scalar field and a vector field can be seen in Figure 2.5. The scalar field is a mappingR2 → R while the vector field is a mapping R2 → R2, i.e. for each tuple(x, y), the fields give a value u and (u, v), respectively.

omega x y u omega x y v u

Figure 2.5: Two examples of fields. A scalar field defined on two-dimensional do-main is depicted on the left. A vector field defined on the same dodo-main can be seen on the right.

2.2.2 Ordinary Differential Equations and Differential and

Algebraic Equations

In some cases a variable that is distributed in space and/or time, e.g. a temperature field, can be approximated by a scalar variable. For example, when the temperature of a fluid in a container is studied, the temperature could be assumed to be the same throughout the entire container at each specific time instant, and one can study the temperature change as a function of time based on the incoming and outgoing fluid flow. Since the simplified model depends only on one variable, time, derivatives that occur are time derivatives, and the equations of the model

(35)

2.2 Differential Equations are ordinary differential equations (ODEs) or differential and algebraic equations (DAEs).

2.2.3 Partial Differential Equations

In more detailed models when the temperature distribution over the container is studied, with different temperature values at different positions inside the container, the temperature is a function of the independent variables. representing the space coordinates inside the container, and also of the time variable if time-dependent be-havior is studied. Derivatives that occur in such models can be partial derivatives , i.e. derivatives with respect to one of the space variables or the time variable. Dif-ferential equations involving partial derivatives are thus called partial difDif-ferential equations (PDEs). If the equations contain the time variable or any time-dependent variable, the problem is called time-dependent, otherwise it is called stationary.

The order of a partial differential equation is defined as the highest differentia-tion order that occurs in that equadifferentia-tion. The most commonly used PDEs in models of many practical systems are second-order PDEs, containing second order deriva-tives. A partial derivative can be stated in different ways. The first-order partial derivative of a dependent variable u with respect to the independent variable x is represented using the following three variants of mathematical notation:

∂u

∂x = ∂xu= ux A second order derivative can be written accordingly as

∂2u

∂x2 = ∂xxu= uxx

The second differentiation can be done with respect to a different variable y than the first differentiation, e.g. ∂xyu= uxy, which is also a second order derivative, a

mixed partial derivative.

A general, second order PDE with an unknown dependent variable u and the independent variables x and y can be stated as

f(u, ux, uxx, uy, uyy, uxy, uyx, x, y) = 0 (x, y) ∈ Ω (2.1) The definition domain, i.e., the geometric region on which the PDE is defined, denotedΩ, is an open2 _{and in most cases bounded region in space as defined by} the independent variables, see Figure 2.6.

2_{An open region is a region where any point in the region can be enclosed by a ball that is completely}

(36)

∂Ω n

Ω

Figure 2.6: The definition domainΩ, its boundary ∂Ω and the outward normal vec-tor n.

2.2.3.1 Vector notation

The differential operators gradient and divergence are commonly used to express PDEs using vector notation. These operators are dimension and coordinate name independent. The gradient operator is defined as follows:

∇u = grad u(x1, . . . , xn) =

∂(u) ∂x₁ . . . ∂(u) ∂xn

The divergence operator is defined as follows: ∇ · u = div u(x1, . . . , xn) = ∂u1

∂x₁ + . . . + ∂un

∂xn

2.2.4 Boundary Conditions

Boundary conditions in a PDE problem specify the behaviour of the model on the boundary of the domainΩ, denoted ∂Ω, and usually contain derivatives of one order less than the PDE. A general boundary condition corresponding to Equa-tion (2.1) can be expressed as

g(u, ux, uy, x, y) = 0 (x, y) ∈ ∂Ω (2.2)

Boundary conditions usually contain a directional derivative instead of the partial derivative with respect to one variable. A common directional derivative is the outward normal derivative, which is obtained by the scalar product of the outward normal on the domain boundary and a vector of the partial derivatives. The outward normal derivative is denoted ∂n and is commonly used in boundary conditions to specify e.g. heat flux through the boundary.

(37)

2.2 Differential Equations Three kinds of boundary conditions that often occur have been given special names:

Dirichlet: u= h₁ on ∂Ω N eumann: _∂n∂ u= h₂ on ∂Ω Robin: _∂n∂ u+ qu = h₃ on ∂Ω

In cases where the right-hand side is zero, the corresponding condition is called homogeneous, otherwise non-homogeneous.

2.2.5 Initial Conditions

Initial conditions are needed in order to find a unique solution to a time-dependent differential equation. Due to the zero derivatives of constants, infinite number of solutions can satisfy a differential equation. Giving the value of the dependent vari-able at some point in time, usually time value where simulation starts, the constant can be determined and a unique solution can be found.

2.2.6 Classification

PDEs are classified based on linearity properties. An equation expressed using the linear operator3_{L() in the form}

L(αu + βv) = f is linear if

L(αu + βv) = αL(u) + βL(v)

If a PDE with the dependent variable u and the independent variables x and y can be expressed in the following form:

auxx+ 2buxy+ cuyy= d (2.3)

where the coefficients a, b and c are constants, the PDE is linear. If the coefficients are functions of the independent variables only, x and y in this example, the PDE is semi-linear, and if the coefficients also depend on u or its first-order derivatives the PDE is quasi-linear. Equations that cannot be expressed linearly in the second order derivatives as in Equation (2.1) are nonlinear. In this text, semi-linear PDEs will be regarded as linear, while quasi-linear PDEs will be classified together with the nonlinear.

(38)

Equation (2.3) is classified further into three categories depending on the value of b2− ac:

Elliptic: b2− ac < 0 P arabolic: b2− ac = 0 Hyperbolic: b2− ac > 0

Elliptic equations typically occur in stationary heat conduction models, whereas parabolic equations arise in time-dependent models. The wave equation used for modeling propagation of waves such as sound waves in gas or electro-magnetic waves is an example of a hyperbolic equation.

2.3 Numerical Solution Methods

A very small number of PDE problems can be solved symbolically, and even fewer of these are solvable with practically useful boundary conditions. For this rea-son, much research has been done on numerical solution of PDEs and different methods have been developed. Unlike the case for solution of ordinary differential equations, there is no unified method for solution of general PDE problems, how-ever, and an appropriate numerical solver must be selected for the specific kind PDE being solved.

The solution method depends on whether a stationary or a time-dependent PDE is solved. Only space discretization is needed for stationary problems, whereas time-dependent problems need both space and time discretization. The methods described below can be used for both these cases. Discretization of only time or space variables can be done as well. For example, an equation can be discretized with respect to space coordinates, leaving a system of algebraic equations contain-ing only time derivatives (the method of lines). Time-dependency can be removed similarly, replacing time derivatives with a finite difference, leaving a stationary PDE problem which can be solved iteratively at different time steps.

Space discretization is performed in different ways, depending on the numerical solution method. The finite difference method requires a grid, which is usually regular, i.e., with equidistant discrete points. This also requires that the definition domain is rectangular. The finite element method can be used on more general domains, because the domain is discretized using a triangular mesh.

2.3.1 Finite Difference Methods

In order to find the unknown function, the domain, i.e. the geometry on which the problem is defined is discretized into a set of grid points and the values of the

(39)

2.3 Numerical Solution Methods function at these points are calculated. When the domain is discretized, the partial derivatives can be approximated by equations using the values at the grid points. Using Taylor’s theorem, different approximations can be derived [15]. The first-order derivative of u with respect to x using a grid with the distance h between the points can be approximated by for instance:

∂u ∂x =

u(x + h) − u(x − h)

2h + O(h2)

Here,O(h2) represents the approximation error. This equation is called a central-difference approximation. Other approximations can be derived, for example

∂u ∂x

u(x + h) − u(x) h which is called a forward-difference formula, and

∂u ∂x

u(x) − u(x − h) h

which is called a backward-difference formula. Higher-order derivatives can be ap-proximated similarly. The derivatives in the PDE are then replaced by the approx-imations and repeated for each grid point, generating an equation system with the values of u at the grid points as unknowns. The boundary conditions are applied at the grid points neighboring the boundary, depending on the kind of boundary con-ditions. With Dirichlet boundary conditions, the values on the boundary are known and can be used directly. With Neumann and Robin conditions, the derivative is ap-proximated in the same way as in the PDE and equations involving the grid points outside the domain are generated, that can be used to express the unknown values on the boundary.

When solving time-dependent problems, the time derivatives can be approxi-mated with finite differences as well. Note however that especially when using finite difference methods, the discretization steps must be chosen according to cer-tain conditions in order to maincer-tain numerical stability during simulation.

2.3.2 Finite Element Methods

A common way of solving a PDE numerically is to approximate the dependent variable u with a functionˆu belonging to a finite dimensional space with known basis functions βk(x), also called trial functions. Using a vector space of

dimen-sion N , the approximate function is defined as

ˆu(x) =

N

k=1

(40)

Each set of coefficients c₁...cN is the coordinates of a function in this space. The

problem is then reduced to finding the coefficients for a function û that satisfy the equation as well as possible, i.e. minimizing the error u− û. Because u is unknown, the error cannot be measured directly. Instead, the partial differential equation is used to measure the error introduced by using the approximation û instead of u. Different measurements lead to different solution methods, such as the least squares method or the weighted residual method. Since there are a finite number of basis functions and finite number of unknown coefficients, it is possible to calculate the coefficients numerically, given the PDE and known values at some points, i.e. the boundary conditions.

The Galerkin method is a basic method for calculating the unknown coefficients. For a partial differential equation of the form

L(u(x)) = g(x)

whereL() is a linear operator containing u and its derivatives, a residual error R of an approximationˆu is defined by

R(ˆu(x)) = L(ˆu(x)) − g(x)

The residual error measures how well the approximation satisfies the partial differ-ential equation. If it is identical to zero, i.e.R(ˆu(x)) ≡ 0 for all x in the domain, then the approximation is the exact solution. A criterion for measuring how close a function is to being zero is to look at its orthogonality to a chosen set of functions. Orthogonality of two functions u(x) and v(x) is defined as

Ω

u(x)v(x) dx = 0

The residual of the exact solutionR(u(x)) is orthogonal to all functions according to this definition. The residual of the approximation is checked against a set of test functions from an appropriate test space. It is sufficient to check the basis vectors of the test space to assure orthogonality to all functions in the test space.

The Galerkin method uses the same space for both trial and test functions, so the basis functions βkare used as test functions as well:

ΩR(ˆu(x))βk(x) dx = 0 k = 1..N

Substituting ˆu with its definition in Equation (2.4) and calculating the integrals gives an equation system of size N that can be solved to finally obtain the coeffi-cients to the trial functions and consequently the approximation to the solution.

One problem with this method is to find basis functions that can be easily com-bined to approximate the unknown function. Basis functions like polynomials that

(41)

2.3 Numerical Solution Methods has global support, e.g. where each basis function contribute to the entire domain, have the problem that, in order to increase the accuracy of the approximation, the degree of the polynomials must be increased, leading to very high degree of poly-nomials.

A more flexible method, called the finite element method, is to divide the domain into M subdomains, called partitions or elements, and use basis functions with only local contribution. Each basis function ϕkis zero in the entire domain except for a

small number of elements. The hat function is an example of a basis function with local support, defined as:

ϕi(x) =      0, x /∈ [x_i−1, x_i+1], x−xi−1 xi−xi−1, x∈ [xi−1, xi], x−xi+1 xi−xi+1, x∈ [xi, xi+1]. (2.5)

With the finite element method, may smaller sized elements can be used for better approximation instead of polynomials of higher order. Element size and shape can be varied as well, using adaptive triangulation/subdivision to create smaller ele-ments for the parts of the domain where the error is too large, for example where the solution has steep changes. Subdivision of the domain in two or three dimen-sions can be performed in different ways, triangulation being a common method.

2.3.3 Method of Lines

The method of lines is a solution method that converts a PDE into a set of ODEs involving only time-dependent functions and time-derivatives. The conversion is done by discretizing the PDE in space, leaving a number of unknowns and their time derivatives. For the space discretization, any of the methods described earlier can be used. For example, if the finite difference method is used, the space dis-cretization leads to one unknown and its time derivative at each grid point on the domain, i.e. a set of ODEs.

One advantage of the method of lines is that advanced numerical solution meth-ods exist for solving general ODEs that do not yet exist for PDEs. There are, for instance, solvers with automatic step adjustment to find a solution with required accuracy. Another advantage is that coupled systems containing both ODE and PDE based models become easier to solve because the space discretization of the PDEs results in ODEs that can be solved together with the already existing ODEs. One disadvantage, though, is that the space discretization is independent of the error controlled step adjustment that is done in the ODE solution process. Thus, even though the ODE solver solves the given ODEs with a desired accuracy, error introduced during the space discretization can be much larger if space discretization is performed without caution.

(42)

Also, the relation between the step sizes in the space and time domain must fulfill certain criteria [16] in order to get a numerically stable solution. However, this is not checked by a pure ODE solver.

(43)

Chapter 3 A Modelica Library for

PDE-based Modeling

“You can’t argue with a machine, Jim!”

Dr. McCoy

This chapter describes a set of packages and models in Modelica to describe and solve PDE-based models.

3.1 Introduction

An object-oriented approach has been taken in the design of the PDE library. As described earlier, all PDE models have certain common components. These are found in this library on an abstract level, without any information about specific models or solution methods. Models and classes that contain specific informa-tion for different soluinforma-tion methods are separated into sub-packages. Here, we first present the continuous parts of the library and show an example of how these can be used to define a PDE-based model. Then, we describe the discrete parts, specif-ically the finite element package which is the main solver package used in the library.

Due to restrictions in the Modelica language, the implementation of the PDE library requires a specific design. Therefore, the design is presented in two steps. First, we describe the general design using classes as they are known in object-oriented design terminology. Then, we describe the packages and models in the library in more detail. Mapping from the general design to a Modelica-specific

(44)

design using packages, called The Package Design, is described separately in order to clarify the concept.

3.1.1 Separation into Continuous and Discrete Parts

The PDE library and the existing models are divided into a continuous and a dis-crete part, in order to separate the mathematical model itself from discretization details. However, only the discrete part of the models contain actual model equa-tions, because of the current lack of Modelica language constructs for expressing space-distributed equations, such as field variables, partial derivatives etc. Such language constructs are described in Chapter 4. Thus, the continuous part contains a number of PDE-based models that can be reused with different parameters (coef-ficients), with the actual equations given in the corresponding discrete models. This approach is similar to the coefficient-based approach, described in Section 7.2, ex-cept that here the discrete models can be fully described in the Modelica language itself instead of in external libraries, allowing Modelica library developers to define PDE-model libraries using the framework.

3.2 Continuous Part Design

The class design of the continuous part of the PDE library is described using object-oriented terminology and UML diagrams. Currently, only two-dimensional do-mains are supported. For discussion of generalization to one and three dimensional models, see Section 3.10. Some UML components used in the following descrip-tion is listed in Appendix A.

The general classes in the PDE library areField, Domain, andBoundary. Some help classes are also included, namelyBoundaryCondition,Point and

BPoint. These classes are independent of the solution methods and specific PDE models. The UML diagram in Figure 3.1 shows an overview of these classes and their relationship.

A field represents a variable that is distributed over a spatial domain. TheField

class is parameterized on the field value type. Thus, the field values can be of any type, including vector-valued types. A field has an operation calledvaluewhich returns the field value at a given point in the domain. Hence, a field can be seen as a mapping from a spatial domain to field values. Each field has exactly oneDomain

component.

TheDomainclass stores spatial domain information associated with each field. A domain is defined as the interior of a spatial region delimited by its boundary. Hence, a spatial domain is implicitly defined by its boundary component.

Each boundary object can be a parametric curve or a number of concatenated boundary parts each of which are parametric curves. TheBoundaryclass

(45)

repre-3.2 Continuous Part Design

Boundary

+shape(in u): BPoint

Domain BoundaryCondition BPoint +x: Integer +y: Integer +bcIndex: Integer Field +value(Point): FieldType FieldType: Point +x: Integer +y: Integer ConstConstField +constval: FieldType +value(Point): FieldType ConstField +value(Point): FieldType

Figure 3.1: Basic classes of the PDE library. See Section 3.4 for description of these packages including their definition in Modelica syntax.

(46)

sents the boundary of a domain using parametric curves or surfaces. A boundary has an operation calledshapewhich defines the shape by mapping parameters to spatial coordinates. In the two dimensional case, the boundary is a curve depend-ing on one parameter, e.g. theshape function is a mapping from R to R2. In order to define a finite domain, a boundary must be closed, e.g. in the two di-mensional case the parameter range[0, 1[ must define a closed boundary, where shape(0) == shape(1).

Other classes areBoundaryCondition,PointandBPoint. BoundaryCon-ditionrepresents one of the included boundary condition types and also con-tains parameters needed by that type of boundary condition. ThePointclass sim-ply represents a coordinate. TheBPointclass represents a point with additional boundary condition information, used for the discretized boundary points.

3.2.1 Standard Boundaries

There are a number of boundary classes defined in the PDE library, in the package

Boundaries. These are shown in Figure 3.2.

Line

+p1: Point +p2: Point +shape(in u): BPoint

Arc

+c: Point +r: Real +a_start: Real +a_end: Real +shape(in u): BPoint

Circle

+shape(in u): BPoint « a_start=0, a_end=2*PI »

Rectangle

+p: Point +w: Real +h: Real

Bezier

+p: array [2..n] of Point +shape(in u): BPoint

Composite

+shape(in u): BPoint parts

1..n Boundary

bottom, right, top, left 4

Figure 3.2: Standard boundary classes in theBoundariespackage.