Spatial Modeling in iThink and STELLA

STELLA and iThink provide capabilities to model spatial problems.  Version 9.1.2 expands these capabilities to allow more intuitive equation formulations.

Spatial modeling is concerned with modeling behavior over space as well as time.  Common applications include modeling flows in lakes and streams or modeling development in an urban landscape.  To demonstrate the new capabilities in version 9.1.2, the one-dimensional diffusion problem will be modeled in STELLA.

Consider the diffusion of heat through a one meter metal bar when a constant temperature is applied to both ends.  This problem is typically introduced as the following partial differential equation for temperature function u(x, t) with the given boundary conditions:

∂u/∂t – k ∂²u/∂x² = 0
u(x, 0) = 0, 0 < x < 1
u(0, t) = u₀, t ≥ 0
u(1, t) = u₀, t ≥ 0

Here, k is the diffusion coefficient (based on thermal conductivity, density, and heat capacity) and u₀ is the constant temperature applied to both ends of the bar.  The finite difference solution to this problem (with the above initial conditions) is:

u(x, t + dt) = u(x, t) + ∂u/∂t
∂u/∂t = k (u(x + dx, t) – 2u(x, t) + u(x – dx, t))/dx²

It is also possible to derive a closed-form solution.  Unless you are an astrophysicist, this solution probably does little to help you understand what is happening.  It is far more intuitive to develop a physical model of the underlying mechanics than to wrestle with the mathematics.

If you consider the metal bar as a single object that is being acted on by its environment, it is possible to consider some of the changes to that bar over time.  If the environment is warmer than the bar, there will be an inflow of heat to the bar.  If the bar is warmer than its surrounding environment, there will be an outflow of heat from the bar.  We probably all know from experience that a hot bar will radiate heat outward in all directions, so in a three-dimensional problem, heat would flow outward in all directions.

If you divide the bar into a set of consecutive finite elements (a finite element is just a fancy word for a small fixed-size section of a larger object – finite because its size is bounded), it becomes possible to imagine a model of what happens to each of these elements individually.  Connecting them in tandem should then yield the behavior of the entire bar.  This is the essence of spatial modeling.

Using the ideas above, we must have some amount of heat entering each element and some amount leaving each element.  For this problem, however, which way does the heat radiate?  The problem is intentionally limited to one dimension, x, so the only choice for the direction of radiation (from an element’s point of view) is left or right.  Thus, we can have heat radiating into the element from either the left or right side, and heat radiating out of the element from either the left or right side.  This is shown in the diagram below and is our model for a finite element of the metal bar.  Note that for clarity u has been renamed T for temperature.

If we chain a number of these elements together horizontally, we will have a model of the entire metal bar.  Using three elements leads to the following model of the metal bar.

This model is obviously greatly simplified.  Note that heat applied left is the heat arriving from the left side of T1.  The heat leaving from the left side is zero, so is not even shown.  An analogous simplification has been made to T3.  This simplified model allows us to consider the nature of each of the four flows in our modeled element, as all of them relate to the neighbors of that element.  In this example, T2 has one neighbor to the left, T1, and one neighbor to the right, T3.  Thus, any heat gained from the left must come from T1 and any heat gained from right must come from T3.  In addition, any excess heat T2 has will radiate out to T1 on the left and T3 on the right.  This simple structure allows us to arrive at the equations for the four flows in our finite element model, shown in terms of both this example and the general finite element model:

in left = C·T1 = C·T(x – dx, t)
out left = C·T2 = C·T(x, t)
out right = C·T2 = C·T(x, t)
in right = C·T3 = C·T(x + dx, t)

C is a diffusion constant describing the rate of temperature loss from a warm body, i.e., the fraction radiated by a finite element.  In the one-dimensional case, C = k/dx².  Note that summing these four flows together properly (i.e., adding the inflows and subtracting the outflows) yields the exact same ∂u/∂t equation derived mathematically and shown above.

But are three elements, two of which are held at the applied temperature, sufficient to properly model a one meter metal bar?  Not really.  The number of elements used defines dx.  Specifically, dx = bar length/number of elements, in this case 1 m/3 = 1/3 m.  Like dt, dx must be small enough to properly catch the dynamics of the system.  We need far more elements to properly model this problem.  Since we have broken the metal bar down into identical finite elements, it becomes natural to build the metal bar up using an array of finite elements rather than showing each individual element.  This model looks like this:

Version 9.1.2 allows the equations above to be naturally expressed as Apply to All equations as follows:

in left = C*T[X – 1]
out left = C*T[X]
out right = C*T[X]
in right = C*T[X + 1]

where X is name of the array dimension.  Note the value one is used instead of dx.  This is because each array element represents a finite element of width dx.  The subscript for the array counts elements rather than distance.  To get the position of the element along the bar it is necessary to multiply the index of the element by dx.  Also note that t is not represented by a subscript in the array.  This is because all models already simulate across time.  Therefore, all entities in a model are functions of time.  Running this model gives the expected behavior:

A final note:  Version 9.1.2 allows any constant, not just one, to be added to, or subtracted from, a dimension name in an equation.  The software evaluates T[X + c] to zero when X + c results in a subscript beyond the bounds of the array.  Thus, when X = 1, the first element, T[X – 1] will be forced to zero.  Likewise, when X = DIMSIZE(), the last element, T[X + 1] will be forced to zero.  As you will discover when you explore it, this particular model (1d-diffusion) forces these end conditions to the constant temperature being applied (100 in this case).

Editor’s Note:  This is part 1 of a 3-part series on spatial modeling in iThink and STELLA.  Part 2 is available here.  Part 3 will be posted the week of April 12th, 2009.