Discrete Simulation of Dynamical Boundary Value ... - Semantic Scholar
1
Discrete Simulation of Dynamical Boundary Value Problems R. Rabenstein a a Lehrstuhl fur Nachrichtentechnik, Universitat Erlangen{Nurnberg, Cauerstr. 7, D-91058 Erlangen, Germany The application of methods from multidimensional systems theory and digital signal processing opens new ways for the ecient simulation of time and space dependent problems like wave propagation or heat and mass transfer. The basic idea is the application of suitably chosen functional transformations for the time and space coordinates. Initial and boundary conditions are considered in a systematic way. This leads to the structure of a discrete system which closely models the behaviour of the given continuous system and which is suitable for computer implementation. Numerical results show considerable savings in computer time over existing numerical methods.
1. INTRODUCTION Problems which depend on continuous variables like time and space are generally modelled by di erential equations. If the quantities in this model are considered as input and output signals, then such an idealized description is also called a continuous system. Purely time dependent problems are described by ordinary di erential equations (ODE), leading to a onedimensional (or lumped parameter) system. Time and space dependent problems like wave propagation or heat and mass transfer are represented by partial di erential equations (PDE), leading to multidimensional (or distributed parameter) systems. Here, initial and boundary values as well as possible excitation functions are the input signals and the desired solution is the output signal. Unfortunately, continuous systems cannot be implemented directly on a digital computer, since no equivalent for the operations di erentiation or integration is available. Instead, it is necessary to derive from the continuous system a discrete one, which contains only shift or delay operations. Frequently used approaches for the discretization of onedimensional systems are transfer function and numerical integration methods, while multidimensional systems are usually treated with nite di erence and nite element methods. This contribution presents a method for the derivation of a discrete model from a continuous multidimensional system, which is based on the use of functional transformations for the independent variables (time and space). It is an extension of transfer functionbased-methods for onedimensional systems. The latter ones will be reviewed shortly in the next section.
2
2. ONEDIMENSIONAL SYSTEMS As an example for a onedimensional system with time as the independent variable, we consider the simple electrical network shown in g. 1. Its mathematical model is the ODE (1), which relates the input and output voltages v(t) and y(t). a is the inverse time constant. Additionally, the voltage y(0) for t = 0 is given as an initial condition. y_ (t) + a y (t) = a v (t); y (0) = yi : (1) This model can be replaced by a system description with transfer functions by application of the Laplace transformation 1 Z Lfy(t)g = Y (s) = y(t)e?
st dt;
Lfy_(t)g = sY (s) ? y(0)
(2)
:
0
Its di erentiation theorem simpli es the ODE (1) in two ways: At rst, it replaces the time derivative in (1) by a multiplication with the frequency variable s. Further, it includes the initial value yi as an additive term. Thus, the ODE with an extra condition (1) is converted into one algebraic equation with no side constraints (s + a)Y (s) = aV (s) + yi : (3) Solving (3) for Y (s) gives the transfer function description shown in gure 1 (top right) with V (s) and yi as inputs and Y (s) as output. V (s) v (t)
y (t)
Laplace transformation
s+a Y (s)
yi
continuous system
a
1 s+a
? discrete system vd (k )
b T
c
Vd (z )
inverse yd (k ) Z-transformation yi
b
z ?c
Yd (z )
z
z ?c
Figure 1. Onedimensional continuous and discrete systems This transfer function description allows an easy transition to discrete systems. One frequently used approach assumes a certain form of the input signal v(t), typically encountered in applications (e.g. impulse train, staircase function, polygon). The corresponding
3 output signal is sampled with step size T and transformed with the Z-transformation to obtain the corresponding transfer function of the discrete system. Fig. 1 shows the result for the so called step invariant transformation with c = exp(?aT ), b = 1 ? c (bottom right), which is custom made for staircase function input. Inverse Z-transformation gives a di erence equation, which can be readily cast into a computer program yd (k ) = c yd (k ? 1) + b vd (k ? 1);
yd (0) = yi :
(4) Fig. 1 shows a block diagram of the discrete system. It consists only of the operations addition, multiplication, and delay by T . Note that this simple system gives the exact output of the continuous system in the sense of yd(k) = y(kT ) for any staircase function v (t) as input signal and any step size T . The application of these methods to onedimensional systems is well known and has been proven to be numerically e ective. (See [3,7] for more detailed descriptions.) It is a characteristic point of this approach, that an appropriate functional transformation takes care of the inital condition in an analytically exact way. We will address now the question, how this method can be extended to multidimensional systems.
3. MULTIDIMENSIONAL SYSTEMS After having considered onedimensional systems, we now look for a similar transformation for the space variables of a multidimensional system. If we were able to include the boundary conditions of the PDE description in a similar way into an algebraic equation as we did for the initial condition, then the above method could be applied also to time and space dependent systems. The e ect would even be more striking, since the boundary conditions have a paramount and lasting in uence on the behaviour of multidimensional systems.
3.1. Continuous system
We will develop this idea by considering the following PDE y_ (x; t) + Lfy (x; t)g = v (x; t); x 2 V ; y (x; 0) = yi (x); x 2 V ; (5) fb fy (x; t)g = (x; t); x 2 S : with the vector x of space variables in the domain V , the operator L of spatial derivatives, and the boundary operator fb de ned on the surface S . fb may describe boundary conditions of the rst, second, or third kind. The excitation function v(x; t), the initital value yi(x), and the boundary values (x; t) are assumed to be known. The operator L determines the type of problem described by (5). For an important class of parabolic problems (heat conduction, mass di usion, current distribution), L takes the form 1 div ((x)grad y(x; t)) : Lfy (x; t)g = ? (6) c(x) In the case of heat conduction, y is the temperature, c the heat capacity per unit volume, and the thermal conductivity.
4 Next, we search for suitable functional transformations, which transform (5) into an algebraic equation which contains the excitation, the initial, and the boundary values. According to our experience with onedimensional systems, the Laplace transformation appears to be a good candidate for the time variable. Indeed, application of (2) to (5) turns the initial-boundary value problem into a boundary value problem with the initial value as an additive term: sY (x; s) + LfY (x; s)g = V (x; s) + yi (x); x 2 V (7) fb fY (x; s)g = (x; s); x2S : We could treat the boundary condition in the same way, if we had a transformation T for the space variables x such that the transform of LfY g could be expressed by the transform of Y and the known boundary values . This requirement can be written in analogy to (1) as
T fY (x; s)g = Y ( ; s); T fLfY (x; s)gg = 2Y ( ; s) ?
ZZ G (x; ) (x; s) dS b
S
(8)
The factor G b (x; ) depends on the exact form of the transformation T , which will now be determined. We show the procedure for the operator L in the form of (6). For simplicity, the temporal frequency variable s is dropped for the moment. We start with the approach T fY (x)g = Y ( ) = c(x) Y (x)K (x; ) dV; (9)
ZZZ V
where the material parameter c is a weighting function and K (x; ) is the kernel of the transformation. The latter shall be chosen such that T fLfY (x)gg takes the desired form of requirement (8). Application of the transform (9) to L in the form of (6) gives
ZZ Z ZZ Z T fLfY (x)gg = c(x) LfY (x)gK (x; ) dV = ? div ((x) grad Y (x)) K (x; ) dV (10) V
V
Using Green's integral theorem and the fact that L is a self adjoint operator, we can split the volume integral over c LfY gK into a volume integral over c LfK gY plus a surface integral. The result does indeed take the desired form of requirement (8), provided that K (x; ) ful lls the conditions ? 2K (x; ) + L fK (x; )g = 0; (11) fb fK (x; )g = 0 : The details of the derivation are omitted for the sake of brevity (see [2,4]). Thus, the transformation kernel K (x; ) can be determined from a boundary value problem (11) with the same structure as the original problem (7). However, conditions (11) di er from (7) in three important aspects: They do not contain the temporal frequency variable s. They describe a homogeneous problem.
5
They belong to the class of Sturm-Liouville problems with the useful properties: { Nontrivial solutions exist only for discrete values ; = 1; 2; 3 : : : of the
frequency variable (Eigenvalues). { The corresponding solutions K (x; ) are mutually orthogonal (Eigenfunctions). Using these results, not only the the transformation kernel can be determined, but also the inverse transformation can be constructed. It is given by an in nite series, due to the discrete nature of the spectrum and the orthogonality of K (x; ) 1 T ?1fY ( )g = Y (x) = 1 Y ( ) K (x; ) : (12)
X
=1
N
is the L2-Norm of the Eigenfunctions. The transform pair (9,12) is also called generalized Fourier transformation or Sturm-Liouville transformation [1]. Its application turns the boundary value problem (9) into an algebraic equation, which can be easily solved for Y ( ) and which contains exactly the given functions on the right hand side of (7) in terms of their transforms N
sY ( ; s) + 2Y ( ; s) = V ( ; s) + y0( ) +
ZZ G (x; ) (x; s) dS b
S
:
(13)
Thus, we have arrived at a transfer function description of the initial-boundary value problem (5) which parallels the corresponding relation for onedimensional systems (see (3) and g. 1). The boundary values show up as an additional term on the right hand side.
3.2. Discrete system
The frequency domain description (13) suggests to apply the same techniques for the time discretisation as we did previously for onedimensional systems. The result is a discrete-time continuous-space system shown in g. 2. The coecients b( ) and c( ) are calculated in the same way as before with a from (3) replaced with 2 from (13). The b( ) y( ; k ) vd(x; k ) yd (x; k ) T T T ?1
c( )
Figure 2. Multidimensional discrete-time continuous-space system simplest way for the discretization of the space variable x is to evaluate (12) at the discrete
6 points of interest. Special measures to avoid convergence problems may be necessary. Fast Fourier transform (FFT) techniques can be applied in the case of equidistant eigenvalues. An alternative is the application of discrete convolution techniques [5,6]. The result is a fully discrete system in form of a multidimensional digital lter. Its inputs are sampled versions of the excitation, initial and boundary values and its output closely resembles the output of the corresponding continuous systems at the sampling points in time and space.
4. EXAMPLE brick
insulation
1
temperature in K
temperature in K
brick
2
0.8 0.6 0.4 0.2
0
temperature
0 0.1 0.2 0.3 space in m
0.4 0
100
200
300
400
500
1.5 1 0.5 0 0 600
0.1
time in s
0.2 0.3 space in m
0.4 0
100
200
300
400
500
time in s
Figure 3. Heat transfer through a multi layer slab As an example for the application, consider the heat ow through a multi layer slab with convection boundary conditions (third kind). Fig. 3 shows the description of the problem as well as the simulation result for step like and periodic variation of the temperature at the left side.
REFERENCES 1. R.V. Churchill. Operational Mathematics. McGraw-Hill, New York, 2. edition, 1958. 2. R.M. Cotta. Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton, 1993. 3. Z. Kowalczuk. Discrete approximation of continuous-time systems: a survey. IEE Proceedings-G, 140(4):264{278, Aug. 1993. 4. M.D. Mikhailov and M.N. Ozisik. Uni ed Analysis and Solutions of Heat and Mass Di usion. John Wiley & Sons, 1984. 5. R. Rabenstein. Simulation of linear continuous systems with distributed parameters. Simulation Practice and Theory, 1:93{107, 1993. 6. R. Rabenstein. Discrete simulation of linear multidimensional continuous systems. Multidimensional Systems and Signal Processing, 5:7{40, 1994. 7. H.W. Schuler. A signalprocessing approach to simulation. FREQUENZ, 35(7):174{184, 1981.
600
Discrete Simulation of Dynamical Boundary Value Problems R. Rabenstein a a Lehrstuhl fur Nachrichtentechnik, Universitat Erlangen{Nurnberg, Cauerstr. 7, D-91058 Erlangen, Germany The application of methods from multidimensional systems theory and digital signal processing opens new ways for the ecient simulation of time and space dependent problems like wave propagation or heat and mass transfer. The basic idea is the application of suitably chosen functional transformations for the time and space coordinates. Initial and boundary conditions are considered in a systematic way. This leads to the structure of a discrete system which closely models the behaviour of the given continuous system and which is suitable for computer implementation. Numerical results show considerable savings in computer time over existing numerical methods.
1. INTRODUCTION Problems which depend on continuous variables like time and space are generally modelled by di erential equations. If the quantities in this model are considered as input and output signals, then such an idealized description is also called a continuous system. Purely time dependent problems are described by ordinary di erential equations (ODE), leading to a onedimensional (or lumped parameter) system. Time and space dependent problems like wave propagation or heat and mass transfer are represented by partial di erential equations (PDE), leading to multidimensional (or distributed parameter) systems. Here, initial and boundary values as well as possible excitation functions are the input signals and the desired solution is the output signal. Unfortunately, continuous systems cannot be implemented directly on a digital computer, since no equivalent for the operations di erentiation or integration is available. Instead, it is necessary to derive from the continuous system a discrete one, which contains only shift or delay operations. Frequently used approaches for the discretization of onedimensional systems are transfer function and numerical integration methods, while multidimensional systems are usually treated with nite di erence and nite element methods. This contribution presents a method for the derivation of a discrete model from a continuous multidimensional system, which is based on the use of functional transformations for the independent variables (time and space). It is an extension of transfer functionbased-methods for onedimensional systems. The latter ones will be reviewed shortly in the next section.
2
2. ONEDIMENSIONAL SYSTEMS As an example for a onedimensional system with time as the independent variable, we consider the simple electrical network shown in g. 1. Its mathematical model is the ODE (1), which relates the input and output voltages v(t) and y(t). a is the inverse time constant. Additionally, the voltage y(0) for t = 0 is given as an initial condition. y_ (t) + a y (t) = a v (t); y (0) = yi : (1) This model can be replaced by a system description with transfer functions by application of the Laplace transformation 1 Z Lfy(t)g = Y (s) = y(t)e?
st dt;
Lfy_(t)g = sY (s) ? y(0)
(2)
:
0
Its di erentiation theorem simpli es the ODE (1) in two ways: At rst, it replaces the time derivative in (1) by a multiplication with the frequency variable s. Further, it includes the initial value yi as an additive term. Thus, the ODE with an extra condition (1) is converted into one algebraic equation with no side constraints (s + a)Y (s) = aV (s) + yi : (3) Solving (3) for Y (s) gives the transfer function description shown in gure 1 (top right) with V (s) and yi as inputs and Y (s) as output. V (s) v (t)
y (t)
Laplace transformation
s+a Y (s)
yi
continuous system
a
1 s+a
? discrete system vd (k )
b T
c
Vd (z )
inverse yd (k ) Z-transformation yi
b
z ?c
Yd (z )
z
z ?c
Figure 1. Onedimensional continuous and discrete systems This transfer function description allows an easy transition to discrete systems. One frequently used approach assumes a certain form of the input signal v(t), typically encountered in applications (e.g. impulse train, staircase function, polygon). The corresponding
3 output signal is sampled with step size T and transformed with the Z-transformation to obtain the corresponding transfer function of the discrete system. Fig. 1 shows the result for the so called step invariant transformation with c = exp(?aT ), b = 1 ? c (bottom right), which is custom made for staircase function input. Inverse Z-transformation gives a di erence equation, which can be readily cast into a computer program yd (k ) = c yd (k ? 1) + b vd (k ? 1);
yd (0) = yi :
(4) Fig. 1 shows a block diagram of the discrete system. It consists only of the operations addition, multiplication, and delay by T . Note that this simple system gives the exact output of the continuous system in the sense of yd(k) = y(kT ) for any staircase function v (t) as input signal and any step size T . The application of these methods to onedimensional systems is well known and has been proven to be numerically e ective. (See [3,7] for more detailed descriptions.) It is a characteristic point of this approach, that an appropriate functional transformation takes care of the inital condition in an analytically exact way. We will address now the question, how this method can be extended to multidimensional systems.
3. MULTIDIMENSIONAL SYSTEMS After having considered onedimensional systems, we now look for a similar transformation for the space variables of a multidimensional system. If we were able to include the boundary conditions of the PDE description in a similar way into an algebraic equation as we did for the initial condition, then the above method could be applied also to time and space dependent systems. The e ect would even be more striking, since the boundary conditions have a paramount and lasting in uence on the behaviour of multidimensional systems.
3.1. Continuous system
We will develop this idea by considering the following PDE y_ (x; t) + Lfy (x; t)g = v (x; t); x 2 V ; y (x; 0) = yi (x); x 2 V ; (5) fb fy (x; t)g = (x; t); x 2 S : with the vector x of space variables in the domain V , the operator L of spatial derivatives, and the boundary operator fb de ned on the surface S . fb may describe boundary conditions of the rst, second, or third kind. The excitation function v(x; t), the initital value yi(x), and the boundary values (x; t) are assumed to be known. The operator L determines the type of problem described by (5). For an important class of parabolic problems (heat conduction, mass di usion, current distribution), L takes the form 1 div ((x)grad y(x; t)) : Lfy (x; t)g = ? (6) c(x) In the case of heat conduction, y is the temperature, c the heat capacity per unit volume, and the thermal conductivity.
4 Next, we search for suitable functional transformations, which transform (5) into an algebraic equation which contains the excitation, the initial, and the boundary values. According to our experience with onedimensional systems, the Laplace transformation appears to be a good candidate for the time variable. Indeed, application of (2) to (5) turns the initial-boundary value problem into a boundary value problem with the initial value as an additive term: sY (x; s) + LfY (x; s)g = V (x; s) + yi (x); x 2 V (7) fb fY (x; s)g = (x; s); x2S : We could treat the boundary condition in the same way, if we had a transformation T for the space variables x such that the transform of LfY g could be expressed by the transform of Y and the known boundary values . This requirement can be written in analogy to (1) as
T fY (x; s)g = Y ( ; s); T fLfY (x; s)gg = 2Y ( ; s) ?
ZZ G (x; ) (x; s) dS b
S
(8)
The factor G b (x; ) depends on the exact form of the transformation T , which will now be determined. We show the procedure for the operator L in the form of (6). For simplicity, the temporal frequency variable s is dropped for the moment. We start with the approach T fY (x)g = Y ( ) = c(x) Y (x)K (x; ) dV; (9)
ZZZ V
where the material parameter c is a weighting function and K (x; ) is the kernel of the transformation. The latter shall be chosen such that T fLfY (x)gg takes the desired form of requirement (8). Application of the transform (9) to L in the form of (6) gives
ZZ Z ZZ Z T fLfY (x)gg = c(x) LfY (x)gK (x; ) dV = ? div ((x) grad Y (x)) K (x; ) dV (10) V
V
Using Green's integral theorem and the fact that L is a self adjoint operator, we can split the volume integral over c LfY gK into a volume integral over c LfK gY plus a surface integral. The result does indeed take the desired form of requirement (8), provided that K (x; ) ful lls the conditions ? 2K (x; ) + L fK (x; )g = 0; (11) fb fK (x; )g = 0 : The details of the derivation are omitted for the sake of brevity (see [2,4]). Thus, the transformation kernel K (x; ) can be determined from a boundary value problem (11) with the same structure as the original problem (7). However, conditions (11) di er from (7) in three important aspects: They do not contain the temporal frequency variable s. They describe a homogeneous problem.
5
They belong to the class of Sturm-Liouville problems with the useful properties: { Nontrivial solutions exist only for discrete values ; = 1; 2; 3 : : : of the
frequency variable (Eigenvalues). { The corresponding solutions K (x; ) are mutually orthogonal (Eigenfunctions). Using these results, not only the the transformation kernel can be determined, but also the inverse transformation can be constructed. It is given by an in nite series, due to the discrete nature of the spectrum and the orthogonality of K (x; ) 1 T ?1fY ( )g = Y (x) = 1 Y ( ) K (x; ) : (12)
X
=1
N
is the L2-Norm of the Eigenfunctions. The transform pair (9,12) is also called generalized Fourier transformation or Sturm-Liouville transformation [1]. Its application turns the boundary value problem (9) into an algebraic equation, which can be easily solved for Y ( ) and which contains exactly the given functions on the right hand side of (7) in terms of their transforms N
sY ( ; s) + 2Y ( ; s) = V ( ; s) + y0( ) +
ZZ G (x; ) (x; s) dS b
S
:
(13)
Thus, we have arrived at a transfer function description of the initial-boundary value problem (5) which parallels the corresponding relation for onedimensional systems (see (3) and g. 1). The boundary values show up as an additional term on the right hand side.
3.2. Discrete system
The frequency domain description (13) suggests to apply the same techniques for the time discretisation as we did previously for onedimensional systems. The result is a discrete-time continuous-space system shown in g. 2. The coecients b( ) and c( ) are calculated in the same way as before with a from (3) replaced with 2 from (13). The b( ) y( ; k ) vd(x; k ) yd (x; k ) T T T ?1
c( )
Figure 2. Multidimensional discrete-time continuous-space system simplest way for the discretization of the space variable x is to evaluate (12) at the discrete
6 points of interest. Special measures to avoid convergence problems may be necessary. Fast Fourier transform (FFT) techniques can be applied in the case of equidistant eigenvalues. An alternative is the application of discrete convolution techniques [5,6]. The result is a fully discrete system in form of a multidimensional digital lter. Its inputs are sampled versions of the excitation, initial and boundary values and its output closely resembles the output of the corresponding continuous systems at the sampling points in time and space.
4. EXAMPLE brick
insulation
1
temperature in K
temperature in K
brick
2
0.8 0.6 0.4 0.2
0
temperature
0 0.1 0.2 0.3 space in m
0.4 0
100
200
300
400
500
1.5 1 0.5 0 0 600
0.1
time in s
0.2 0.3 space in m
0.4 0
100
200
300
400
500
time in s
Figure 3. Heat transfer through a multi layer slab As an example for the application, consider the heat ow through a multi layer slab with convection boundary conditions (third kind). Fig. 3 shows the description of the problem as well as the simulation result for step like and periodic variation of the temperature at the left side.
REFERENCES 1. R.V. Churchill. Operational Mathematics. McGraw-Hill, New York, 2. edition, 1958. 2. R.M. Cotta. Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton, 1993. 3. Z. Kowalczuk. Discrete approximation of continuous-time systems: a survey. IEE Proceedings-G, 140(4):264{278, Aug. 1993. 4. M.D. Mikhailov and M.N. Ozisik. Uni ed Analysis and Solutions of Heat and Mass Di usion. John Wiley & Sons, 1984. 5. R. Rabenstein. Simulation of linear continuous systems with distributed parameters. Simulation Practice and Theory, 1:93{107, 1993. 6. R. Rabenstein. Discrete simulation of linear multidimensional continuous systems. Multidimensional Systems and Signal Processing, 5:7{40, 1994. 7. H.W. Schuler. A signalprocessing approach to simulation. FREQUENZ, 35(7):174{184, 1981.
600
