# Mass Transport

## Introduction

This tutorial gives an introduction to modeling mass transport of diluted species. Equations and boundary conditions that are relevant for performing mass transport analysis are derived and explained.

Mass transport is a discipline of chemical engineering that is concerned with the movement of chemical species. The two mechanisms of mass transport are mass diffusion and mass convection. The driving force behind a mass diffusion is the difference in a species concentration at different locations. Mass convection, on the other hand, only occurs when species are transported in a moving fluid medium. The combination of these two mechanisms leads to changes in the species concentration field over time and is modeled with a mass balance equation.

The modeling process results in a partial differential equation (PDE) that can be solved with NDSolve. Furthermore, in this tutorial different types of mass transport boundary conditions are introduced. For the boundary conditions given the modeling of various real world chemical species interactions is shown.

This notebook provides an introduction overview of how to model mass transport. The focus is to present simple yet illustrative examples. Extended examples that show industrial applications of mass transport modeling can be found in the Mass Transfer Model Collection.

Many of the animations of the simulation results shown in this notebook are generated with a call to Rasterize. This is to reduce the disk space this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high quality graphics remove or comment out the call to Rasterize.

To get high fidelity visualizations comment out the rasterization process.

The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.

## Boundary Conditions in Mass Transport

The most common boundary conditions in mass transport modeling can be modeled with DirichletCondition, NeumannValue and PeriodicBoundaryCondition, and can be categorized in the following four types:

• Dirichlet type boundary conditions. This type of boundary condition specifies the species concentration at the boundary, and can be modeled with DirichletCondition.
• Neumann type boundary conditions. This type of boundary condition specifies the mass flux at the boundary, and can be modeled with NeumannValue.
• Robin type boundary conditions. This type of boundary condition specifies the relation between the species concentration and its normal derivatives at the boundary, and can modeled with a NeumannValue since Robin type boundary conditions are technically generalized Neumann boundary conditions.
• Periodic boundary conditions. This type of boundary condition specifies the species concentration at one part of the boundary to be the same at another part, and can be modeled with PeriodicBoundaryCondition.

Under these four types, the following boundary conditions are introduced:

• Dirichlet Type
• Robin Type
• Periodic Type

The following section describes several physical boundary conditions common in modeling mass transport and how they can be represented with the use of DirichletCondition, NeumannValue and PeriodicBoundaryCondition. For this purpose the boundary condition currently discussed is always at the left hand side of the example simulation domain. In some examples additional boundary conditions are at the right hand side to better demonstrate the behaviour of the boundary condition at the left hand side.

Note that the setup of Neumann type boundary conditions will be slightly different between the conservative model and non-conservative model. The detail of this difference is presented in the following section.

### Neumann Values for Conservative and Non-conservative Models

In a previous section the derivation and setup of the conservative and non-conservative mass transport model was presented. The difference of the model formulation has an impact on how Neumann type boundary conditions are to be set up and what they mean.

With the conservative formulation the NeumannValue specifies the boundary values for . For NeumannValue[val,XΓb] this means:

With the non-conservative formulation the NeumannValue specifies the boundary values for . For NeumannValue[val,XΓb] this means:

Note that for the conservative formulation the term: , which is the mass transported by a flow of the medium across the boundary, is also included in the Neumann value.

In some of the boundary conditions presented the form of the Neumann value needs to changed. This means that, for example a Neumann value in the conservative form like eqn (1) will have to be transformed into

Transformations of this sort pose a problem as usually the normal in a Neumann formulation like the left side of eqn (2) is never actually computed but replaced with the value of the right hand side of eqn (3). Once a normal appears on the right hand side of the equation that value of the normal will have to be computed and this is done automatically. If the boundary unit normal is simple it can be specified with the model parameter "BoundaryUnitNormal". This may increase computational efficiency both in speed and memory usage.

Before the mass transport boundary conditions for both conservative model and non-conservative model are presented some example model parameters are set up.

### Model Parameter Setup

The following model parameters are made use of in the examples of mass transport boundary conditions. These parameters define the simulation domain and the simulation end time .

Set up model parameters for the domain and the simulation end time .

In some examples a smoothed step function will be used to prescribe a time profile for a transient parameter such as for example the mass flux or the concentration on the boundary. The smoothed step function is defined as follows:

Here the minimum value and the maximum value of the function can reach are denoted by and . The location of the step is controlled by and the smoothed steepness is controlled by .

Define and visualize the smoothed step function .

### Concentration Boundary Condition

#### Purpose

The purpose of a concentration boundary condition is to set a specific species concentration on some part of the boundary.

#### Formulation

With a specified concentration on the boundary , the concentration boundary condition is given for both conservative and non-conservative models by:

Since the concentration boundary condition is a Dirichlet type condition and thus is not associated with the model equation, its formulation is the same for both the conservative model and non-conservative model.

A concentration boundary for the dependent variable modeled with DirichletCondition.

#### Derivation

A given species concentration prescribed on a boundary is called a concentration boundary condition. This prescribed concentration can either be a constant or time-dependent value, and is set with a DirichletCondition in the mass transport PDE model.

#### Example

To model, for example, the supply of a given species on a boundary, a transient species concentration can be set up at the left end. Note that a Neumann zero condition is implicitly applied at the right end as an outflow boundary.

Here a smoothed step function is used to described the profile of the species concentration on the boundary from to . The parameters and are arbitrarily chosen to simulate the supply of species from the boundary.

Specify the species concentration on the boundary using a smoothed step function.
Set up the concentration boundary condition with DirichletCondition.

Next, set up the mass transport model with a flow field , an initial concentration field and a diffusion coefficient .

Set up the transient mass transport PDE.
Solve the PDE with NDSolveValue.

#### Visualization

Make an animation of the solution using Plot and ListAnimate.

The simulation begins with an undisturbed domain where . As the concentration increases at the left boundary, the excess species are then transported to the right and brings up the overall species concentration within the domain. The speed of the mass transport depends on the species diffusivity and the flow field.

### Outflow Boundary Condition

#### Purpose

If the mass transport occurs in a system where the flow velocity , then an outflow boundary condition is used to model an outlet where the species can leave the domain with the fluid flow.

#### Formulation

When modeling mass transport in a fluid medium, the outflow boundary condition at the outlet is given by:

An outflow boundary for the conservative mass transport model using NeumannValue.
An outflow boundary for the non-conservative mass transport model using NeumannValue.

Note that for non-conservative model an outflow boundary condition is essentially a Neumann zero condition, which means it will be implicitly applied if no boundary condition is specified at a given boundary.

#### Derivation

When modeling mass transport with a fluid flow, the diffusion mass flux at the flow outlet is set to zero to ensure that the mass transferred at the outlet boundary is by mass convection only and the mass diffusion is ignored: .

With the non-conservative formulation, the NeumannValue[val,XΓb] specifies the boundary values (4): . So for it is required that be . The outflow boundary condition for the non-conservative mass transport model is given by:

With the conservative formulation the NeumannValue[val,XΓb] specifies the boundary values (5): . For it is required that:

The outflow boundary condition can only be applied on fully-developed flows. That is, at the flow outlet the velocity profile is unchanging in the flow direction.

In a case where there is recirculation through the outlet boundary, which often happens for turbulent flow, the reentering flow will affect the concentration field of the flow inside the domain and break the zero diffusion flux assumption. In this situation the outflow boundary condition is no longer applicable.

#### Example

In the following example an outflow boundary condition will be set at the left end to model the outgoing mass flux from the domain. The setup for both conservative model and non-conservative model is presented.

To highlight the effect of an outflow boundary, assume the species diffusivity to be zero. Therefore the mass transport is by convection only with the fluid flow .

Set up of model variables and parameters.

A concentration boundary condition is used to model the supply of a species into the domain from the right side.

Set up of a concentration condition.

Next, the boundary condition setup for both conservative model and non-conservative model is presented. Note, that besides the added conservative model form, both setups that are to follow, are the same.

##### Conservative Mass Transport Model with an Outflow Boundary Condition
Set up the conservative form.
Set up the outflow boundary condition for conservative mass transport model with a left-going flow of and a normal at the outlet.
Set up the conservative mass transport PDE.
Solve the PDE with NDSolveValue.

A message is generated to indicate that the PDE model is convection dominated. This is expected since there is no diffusion mass flux in this case, but is stable for this example. More information about this issue can be found in the Finite Element Usage Tips tutorial.

##### Non-Conservative Mass Transport Model with an Outflow Boundary Condition
Set up the default non-conservative form.
Set up the outflow boundary condition for non-conservative mass transport model.

If no boundary condition is specified on any part of the boundary then by default a Neumann zero boundary condition is implicitly used. This implies that an outflow boundary is the default boundary condition used for non-conservative mass transport model.

Set up the transient non-conservative mass transport PDE.
Solve the PDE with NDSolveValue.

#### Visualization of an Outflow Boundary Condition

Make an animation of the solution using Plot and ListAnimate.

With an outflow boundary condition applied at the left boundary, the species was transported out of the domain by the fluid flow without reflection. Note that the results of conservative model and non-conservative model are consistent.

### Mass Flux Boundary Condition

#### Purpose

The purpose of a mass flux boundary condition is to model the amount of a given species entering or leaving the domain through some part of the boundary.

#### Formulation

With a prescribed mass flux on the boundary , the mass flux boundary condition is given by:

A mass flux boundary for the conservative mass transport model using NeumannValue.
A mass flux boundary for the non-conservative mass transport model using NeumannValue.

#### Derivation

A boundary where the mass flux normal to the boundary is specified and not equal to zero is called a mass flux boundary:

By convention a negative sign is added in front of to indicates that the mass flux is specified opposite to the outward normal . Therefore, a positive value of denotes the inward mass flux where a given species enters the domain, and a negative denotes an outward flux.

As described in the previous section Transport Equation Derivation, the mass flux consists of the diffusion flux and the convection flux :

Recall that for the conservative formulation, the NeumannValue[val,XΓb] specifies the boundary values (6): . Therefore, by inserting (7) into (8), the mass flux boundary condition for the conservative model is given by:

In this case is and this results in NeumannValue[(t,X),XΓb].

For the non-conservative formulation the NeumannValue[val,XΓb] specifies the boundary values (9): . So by casting the boundary convection flux in (10) to the right hand side, the mass flux boundary condition for the non-conservative model is given by:

The unit of a mass flux depends on the dimension of the boundary. In 1D (), 2D () and 3D domain (), has a unit of , and , respectively.

#### Example

In the following example a transient mass flux is applied at the left boundary to model the supply of a species into the domain with no actual chemical reaction involved. The right boundary is set up as an outflow boundary condition to model the outgoing mass flux at the right end.

The profile of the mass flux is defined as:

Set up the transient mass flux.
Set up of model variables and parameters.

Next, the boundary condition setup for both conservative model and non-conservative model is presented. Note, that besides the added conservative model form, both setups that are to follow, are the same.

##### Conservative Mass Transport Model with a Mass Flux Boundary Condition
Set up the conservative form.
Set up the mass flux boundary condition for conservative mass transport model.

Specify a flow field and an outflow boundary condition to model the outgoing mass flux at the right end.

Set up an outflow boundary condition for conservative mass transport model.

The initial concentration field and the diffusion coefficient are given by and , respectively.

Set up the time-dependent conservative mass transport PDE.
Solve the PDE with NDSolveValue.
##### Non-Conservative Mass Transport Model with a Mass Flux Boundary Condition
Set up the default non-conservative form.
Set up the mass flux boundary condition for non-conservative mass transport model with a flow velocity and a normal at the outlet.
Set up the transient non-conservative mass transport PDE.

Note that for non-conservative model an outflow boundary condition is a Neumann zero condition, and is implicitly applied at the right end of the domain.

Solve the PDE with NDSolveValue.

#### Visualization of a Mass Flux Boundary Condition

Make an animation of the solution using Plot and ListAnimate.

With a mass flux applied at the left boundary, the species are gradually transported into the domain. After the mass flux is turned off at time , the remaining species continue to be transferred out of the domain with the fluid flow and the overall concentration field is brought down. Note that the results of the conservative model and the non-conservative model are consistent.

### Impermeable Boundary Condition

The impermeable boundary condition is a special case of the mass flux boundary condition for the case where the flux across the boundary is .

#### Purpose

An impermeable boundary condition is to model a boundary where a species can not pass through and there is no mass flux across it.

#### Formulation

An impermeable boundary condition is given by:

An impermeable boundary for the conservative mass transport model using NeumannValue.

Note that for a conservative model an impermeable boundary condition is essentially a Neumann zero condition, which means it will be implicitly applied if no boundary condition is specified at a given boundary.

An impermeable boundary for the non-conservative mass transport model uses NeumannValue.

#### Derivation

An impermeable boundary condition denotes a boundary where there is no mass flux across it:

Inserting (11) into the mass flux boundary condition (12), yields the formulation of an impermeable boundary condition for both conservative and non-conservative models.

Impermeable boundary condition for the conservative mass transport model.
Impermeable boundary condition for the non-conservative mass transport model.

#### Example

In the following examples a Gaussian distribution is used to describe the initial concentration field within the domain, and an impermeable boundary is positioned on both ends to prohibit mass flux across boundaries.

Set up a Gaussian distribution as the initial concentration field .

To highlight the effect of an impermeable boundary, assume that the flow field , which means that the mass transport is by diffusion only. The diffusion coefficient of the species is given by .

Set up of model variables and parameters.

Next, the boundary condition setup for both conservative model and non-conservative model is presented. Note, that besides the added conservative model form, both setups that are to follow, are the same.

##### Conservative Mass Transport Model with an Impermeable Boundary Condition
Set up the conservative form.
Set up the impermeable boundary condition for conservative mass transport model.

If no boundary condition is specified on any part of the boundary then by default a Neumann zero boundary condition is implicitly used. This implies that an impermeable boundary is the default boundary condition when used for a conservative mass transport model.

Set up the time-dependent conservative mass transport PDE.
Solve the PDE with NDSolveValue.
##### Non-Conservative Mass Transport Model with an Impermeable Boundary Condition
Set up the default non-conservative form.
Set up the impermeable boundary condition for non-conservative mass transport model.
Set up the transient non-conservative mass transport PDE.
Solve the PDE with NDSolveValue.

#### Visualization of an Impermeable Boundary Condition

Make an animation of the solution using Plot and ListAnimate.
Compute the net concentration at and for both models.

It is seen that the initial concentration field was gradually smoothed out by the internal diffusion, however, the net concentration within the domain remains unchanged. This is because there is no mass flux across the impermeable boundary on both ends. Note that the results of conservative model and non-conservative model are consistent.

### Surrounding Flux Boundary Condition

A surrounding flux boundary condition is a flux condition applicable in the special case when the flow velocity field is zero . The surrounding flux boundary condition then extends a mass flux boundary condition for the case where the flux across the boundary depends on an exterior concentration and a mass transfer coefficient by: .

#### Purpose

A surrounding flux boundary condition is used to model mass transfer between the modeled system and the surrounding environment via diffusion across a boundary where the flow velocity is zero .

#### Formulation

Given the profile of an exterior concentration and a mass transfer coefficient on the boundary , the surrounding flux boundary condition is given by:

Since there is no fluid flow involved in the model, the formulation of the surrounding flux boundary condition is the same for both the conservative model and non-conservative model.

A surrounding flux boundary for the dependent variable modeled with NeumannValue.

#### Derivation

When there is no fluid flow at the boundary, diffusion becomes the only mechanism for species to pass through the boundary. The rate of the diffusion flux across a boundary depends on the concentration gradient and physical properties of the transported species and properties of the medium, such as its phase and density.

Unfortunately, the relationship between the diffusion flux and these physical properties is not easily determined. To work around that the mass transfer coefficient is defined to lump these factors together. The diffusion mass flux can then be expressed as:

where is the exterior concentration in the surroundings of the modeled domain. The mass transfer coefficient can be determined experimentally [13, 14]. The mass transfer coefficient typically ranges from for gas phases to for liquid phases.

Inserting (15) into the mass flux boundary condition (16) and set the flow field , then the surrounding flux boundary condition can be written as:

#### Example

Consider a 1D example where carbon dioxide, , diffuses out of the domain from the left and right boundaries. The concentration outside the domain is comparably dilute and is assumed to be . The diffusion mass flux through the boundaries is modeled by a surrounding flux boundary condition with a given mass transfer coefficient .

In this case a Gaussian distribution is used to describe the initial concentration field , and the diffusion coefficient of is given by .

Set up of model variables and parameters.
Specify an exterior concentration and the mass transfer coefficient on the boundary.
Set up the transient mass transport PDE.
Solve the PDE with NDSolveValue.

To better understand the effects of the surrounding flux boundary, compare the result of an analytical solution without any boundary condition. That is, the molecules keep diffusing out as if it had an infinite extent of the domain.

Solve the PDE analytically with no restriction on the boundary.

#### Visualization of a Surrounding Flux Boundary Condition

Make an animation of the solution using Plot and ListAnimate.

With surrounding flux conditions applied on both ends, the molecules were transported out of the domain at a faster rate than the one without it, which means the mass diffusion at boundaries is more efficient than the internal diffusion in this case.

Note that the diffusion flux to the surrounding environment varies from case to case, and is controlled by the values of the exterior concentration and the mass transfer coefficient .

Compare the net concentration at and .

With the surrounding flux condition the net concentration has been reduced from roughly to in three hours.

### Symmetry Boundary Condition

#### Purpose

A symmetry boundary condition is used to reduce the extent of the computational domain to a symmetric sub-domain of the full physical model geometry by effectively reflecting across a linear boundary of the simulation domain where the flow velocity normal to the boundary is zero. This allows for a faster solution process with a lower memory requirement.

#### Formulation

The symmetry boundary condition is given by:

The formulation of the symmetry boundary condition is the same for both the conservative model and non-conservative model. Note that there is no boundary convection flux term in (17) since the normal flow velocity on the symmetry boundary will remain at zero at all time.

A symmetry boundary for the dependent variable modeled with NeumannValue.

Note that a symmetry boundary condition is essentially a Neumann zero condition, which means it will be implicitly applied if no boundary condition is specified at a given boundary.

#### Derivation

To maintain the symmetry across a linear boundary, there should be no mass flux and no flow normal to the boundary at any time:

By inserting (18) into the mass flux boundary condition (19), the surrounding flux boundary condition can be written as:

#### Example

Consider the case of solving the concentration field of a 1D system from to . If the concentration pattern is expected to have a mirror symmetry along , the simulation domain can effectively be constructed with only the left half of the system. Then a symmetry boundary condition should be applied at .

Due to the symmetry, the concentration gradient and the flow velocity on the symmetry boundary will remain at zero at all time. Strictly speaking, a symmetry boundary is essentially an impermeable boundary. By setting the flow velocity to zero in the impermeable boundary condition (20), the symmetry boundary condition can be written as:

### Periodic Boundary Condition

#### Purpose

The purpose of a periodic boundary condition is to map the species concentration from one boundary to another in order to exploit the periodicity of the domain.

#### Formulation

Given a function that maps the species concentration from the periodic boundary to the targeted boundary , the periodic boundary condition can be written as:

The formulation of a periodic boundary is the same for both the conservative model and non-conservative model.

A periodic boundary condition for the dependent variable modeled with PeriodicBoundaryCondition.

#### Derivation

A periodic boundary condition is applied to compute mass transfer in spatially periodic domains. Given a targeted boundary , the species concentration on a periodic boundary can be mapped to the concentration on the targeted surface by a prescribed function . The boundary condition is set by the PeriodicBoundaryCondition in the mass transport PDE model.

#### Example

As an example species transport within a circular tube is presented. With the usage of the periodic boundary condition it is possible to perform the simulation in a 1D domain.

To highlight the effect of a periodic boundary, assume the species diffusivity to be zero. Therefore the mass transport is by convection with the fluid flow only.

Set up model variables and the species diffusivity and the flow velocity .

The circular tube is converted into a 1D model of the length equals to the perimeter of the tube. The periodic boundary is set at the left end so that when the species pass through the right side of the domain , it re-appears at the left side with the same magnitude.

Set up the PeriodicBoundaryCondition at the left end with the mapping function .

Next, a Gaussian distribution is used as the initial concentration field .

Set up the initial concentration field .
Set up the transient mass transport PDE.
Solve the PDE with NDSolveValue.

A message is generated to indicate that the PDE model is convection dominated. This is as expected since there is no diffusion mass flux in this case. More information about this issue can be found in the Finite Element Usage Tips tutorial.

#### Visualization of a Periodic Boundary Condition

Make an animation of the solution using Plot and ListAnimate.

It is seen that the species was transported to the right by the fluid flow within a spatially periodic domain. As the species pass through the right boundary it re-appears at the left side because of the periodicity of the domain set up with the periodic boundary condition. Since the species diffusivity is set as zero and the mass transported is by convection only, the pattern and the magnitude of the species concentration remains constant at all times.

## Appendix

### Special Cases of the Mass Balance Equation

As presented in the section Mass Balance Equation Derivation, the conservative and non-conservative forms of the mass balance equation are given by:

The following sections present special cases of the mass balance equation expressed in cylindrical coordinates and spherical coordinates.

#### Stationary Case

If the concentration field is in a steady state, the transient terms in the mass balance equation (21) and (22) vanish and the mass balance equation simplifies to:

#### Mass Transport by Diffusion Only

If the mass transport occurs in a solid medium, then, because a solid can not have an internal velocity field by definition, the convection terms in the mass balance equation (23) and (24) are set to zero, yielding:

In this case the mass transport is by diffusion only.

#### Mass Transport by Convection Only

When modeling mass transport by fluid flow, the diffusion term may be ignored if the mass transport is dominated by the mass convection, leading to:

In this case the mass transport is by convection only.

In a situation where the convection components and become large, the mass transport model will result in convection dominated PDEs. The solution of this type of PDEs may become unstable and extra stabilization methods may be needed. More information about this issue can be found in the Finite Element Usage Tips tutorial.

#### Mass Balance Equation in Cylindrical Coordinates

When modeling a mass transport problem, sometimes it is not convenient to describe the model in Cartesian coordinates. The heat equation may also be expressed using a cylindrical or spherical coordinate system.

A graphics showing cylindrical coordinates:

In the cylindrical coordinate system , and denote the radial, azimuthal and vertical directions, respectively. In terms of the Cartesian coordinates , the cylindrical coordinates are defined by:

or

By inserting the coordinates relation (25) into equation (26) and (27), the mass balance equation can be expressed in cylindrical coordinates as:

If the mass transport within a model is rotationally symmetric about the axis, the resulting concentration field will be invariant in the direction. Equation (28) then simplifies to:

In that case, a 3D mass transport problem can be modeled in a 2D domain by making use of this symmetric property. This type of model is known as an axisymmetric model.

#### Mass Balance in Spherical Coordinates

A graphics showing spherical coordinates:

In the spherical coordinate system , and denote the radial, azimuthal and polar directions, respectively. In terms of the Cartesian coordinates: , the spherical coordinates are defined by:

or

By inserting the coordinates relation (29) into equation (30) and (31), the mass balance equation can be expressed in spherical coordinates as:

## References

1.  Fick, A. Über Diffusion. Annalen der Physik (in German), 94 (1): 5986 (1855).

2.  Wolf, E. E. and Alfani, F. Catalysts Deactivation by Coking, Catalysis Reviews: vol 24 329-371 (1982).

3.  Calvo, E. G., Arranz, M. A. and Leton, P. Effects of Impurities in the Kinetics of Calcite Decomposition, Thermochimica Acta. 170: 711 (1990).

4.  Fu, J. C., Hagemeir, C. and Moyer, D. L., An unified mathematical model for diffusion from drug polymer composite tablets, Journal of Biomedical Materials Research: vol. 10, no. 5, pp. 743758 (1976).

5.  Parra-Guevara, D. and Skiba, Y. N. Elements of the mathematical modeling in the control of pollutants emissions, Ecological Modelling: vol. 167, no. 3, pp. 263275 (2003).

6.  Goldstein, R.J. and Cho, H. H. A review of mass transfer measurements using naphthalene sublimation, Experimental Thermal and Fluid Science: vol. 10, no. 4, pp. 416434 (1995).

7.  Cañizares, P., García-Gómez, J., Fernández de Marcos, I., Rodrigo, M. A. and Lobato, J. Measurement of Mass-Transfer Coefficients by an Electrochemical Technique, Journal of Chemical Education: vol. 83, no. 8 (2006).

8.  LeVeque, R. J. Numerical Methods for Conservation Laws, Birkhduser (1992).

9.  Whitman, G. Walter. The two film theory of gas absorption, International Journal of Heat and Mass Transfer: vol. 5, no. 5, pp. 429-433 (1962).

10.  Prausnitz, J. M., Lichtenthaler R. N. and de Azevedo, E. G. Molecular Thermodynamics of Fluid Phase Equilibria 3rd Ed., Prentice Hall PTR, New Jersey (1999).