« go back

Analytical solutions for compartmental models of contaminant transport in enclosed spaces

Understanding the transport of hazardous airborne materials within buildings and other enclosed spaces is important for predicting and mitigating the impacts of deliberate terrorist releases of chemical and biological materials. Because such materials may be acutely toxic or infectious it is important to understand how concentrations may change with time to understand the hazards that different scenarios may pose. It is also relevantto the study of accidental releases of industrial materials and the impact ofenvironmental pollutants on indoor air quality.

A range of different numerical modelling approaches are regularly used to study these problems as well as experimental methods. Computational fluid dynamics (CFD) can be used for detailed studies of air and contaminant movement within enclosed spaces. However, CFD methods require highly detailed input data and have significant model creation and execution times. This can make them impractical for whole building studies in some cases.

Multizone models (CONTAM [1], COMIS [2] etc.) provide an alternative approach where the building is divided into a series of well-mixed volumes connected by paths through which air and contaminants can pass. These models have the advantage that they arequicker to execute than CFD models and typically require less input information. Thecontaminants are normally considered dilute in such approaches. Typical model size isof the order of 10-100 zones, although 1000 or more may be required in some cases.

Multizone models solve the air flow through the network (typically using a non-linear pressure solver) for a series of quasi-steady states. The contaminant dispersion resulting from the air flow solution and the contaminant initial and boundary conditions isthen calculated. Whilst these models are well developed and have been validated for arange of studies, they rely on numerical methods which can become time consuming forlarge studies (e.g. Monte Carlo analysis). In addition, little insight is gained into thesystem behaviour using a numerical approach.


Problem presented by:

Simon Parker (DSTL)

Study Group contributors:
Elif Baskaya (Karadeniz Technical University)
Stephen Cowley (University of Cambridge)
Ishak Cumhur (Karadeniz Technical University)
Patricio Farrell (University of Oxford)
Jens Gravesen (Technical University of Denmark)
Dariusz Grala (University of Warsaw)
Si^an Jenkins (University of Bath)
James Kirkpatrick (University of Oxford)
Karin Mora (University of Bath)
Margarita Schlackow (University of Oxford)
Suleyman Sengul (Karadeniz Technical University)
Xiaoting Wang (University of Cambridge)
Dave Wood (University of Warwick)
Tristram Armour (Industrial Mathematics KTN)


Related resources:

Problem brief

Initial presentation

Final presentation

Final report

Other defense projects

Other Study Group projects



Archived discussion

Reply by Tristram Armour on April 4, 2011 at 10:24pm

Here are photos of the whiteboard on Monday evening.


P1010331.JPG, 2.3 MB

P1010332.JPG, 2.4 MB


Reply by Tristram Armour on April 5, 2011 at 9:56am

The questions that Simon want answered is in the attached picture.

Calling numerical analysts... (we do not have one here)

- How do you solve for x(t) when A is not diagonalisable? (we think there is a standard solution)

- solving for when the number of rooms/zones large? (perhaps caused by these other problems we see)

Attachments: P1010334.JPG, 2.6 MB


Reply by Shengxin Zhu on April 5, 2011 at 12:42pm

[1] eigenvalue location....

The smallest/largest magnitude eigenvalue could be estimite by Gershgorin Circle/Disc Theroem. The smallest maginitude real eigenvalue could be -max{(Q_{jj}+Q_{o,j})/v_{j,j}}.; and the largerst real part of the eigenvalue is bounded by -max{Q_0j/v_{jj}}

In fact. what importat for the eigenvalue is onlythe real part of the eigenvalue, we don't care about the image part of the eigenvalue.and the real part of the eigenvalue locate in a interval , [-max{(Q_{jj}+Q_{0,j})/V_{jj}},max{Q_{0,j}/V_{j,j}}.

You can also indicate that there is very likely the the matrix A is sigular /indefinite or with eigenvalues with positive part.

[Ref: Any Numerical Linear Algebra. or Matrix Analysis (by, Horn and Johnson)

[2] one possible case could ensure that the matrix A is diagonalisable.

that is some postive number k_{i}, which satisfied that the Q_{ij}/k_{i}=Q_{ji}/k_{j}. In this case, it could be ensure that the matrix A could expressed as a digonal matrix mutiply a symmetric matix. which implies that the matrix A is selfadjiont with a diagonal postitive matrix. For Matrix Theory [Ref; Horn and Johnasn Matrix Analysis I , II, mightbe chapter 4 or 6] and Thus A is diagonalsiable. Further more with all reall eigenvalues.

[3] Nolinear term could be added and easily handeled both analysically or numerically.

For analytically, I perfer using Laplacian Transform methods to solve such problem. For example the nolinear term like Gaussion .

For Numerically, the difficult lies on not only computing the expentiall of matrix A. but also lies in the compute the Exp (A)-I, as we said before A might be sigular and it very likely that there is a case 1-1. However, the recent work of Nick Higham on matrix function and Nick Trefethen could be easily handle this porblem.

(Ref: A.Kassam, L.N. Trefethen; Fourth-order time-stepping for stiff PDEs, SIAM J Sci, Comput 2005 (20)4 1214-1233.) In the time direction the Stiff PDEs reduces to sitff ODEs. So the paper suggest stable and fast way to compute the ODE system like

u'=Lu +Nu (Linear term + nolinear term.)

Reply by Shengxin Zhu on April 5, 2011 at 2:47pm

[1]Thanks for David during the lunch break time, I found the The matrix is digonal dominate. also usng the Gershgorin Circle theorem. you can locate the eigenvalue. all in the Left half plane. might be zeros.


Reply by Tristram Armour on April 5, 2011 at 3:05pm

Thanks for your comments - indeed we are looking particularly hard at computing the exponential of A, "19 dubious ways of calculating matrix exponentials" is particularly useful, and we are currently trying to get some code from Nick's work that creates the block triangular form of A.


Reply by Paul Dellar on April 5, 2011 at 4:06pm

You can exponentiate A analytically from its Jordan canonical form, as in section 16 of the reference below. This won't work numerically because the Jordan canonical forms for exactly equal and almost equal eigenvalues are completely different. Numerically, scaling and squaring (to use the Taylor series definition of exp(tA)) is probably as good as anything. That shouldn't care whether the marix is diagonalisable. If A is large and sparse, there is a package called expokit for computing exp(tA) v for a specific vector v using Krylov-space methods. These days "large" probable means "bigger than 1000 by 1000".

Moler & Van Loan (2003) Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later, SIAM Rev. 45, 3-49,

No comments yet. Be the first.

Most read articles

Maintenance inspection design software for full-scale industrial systems

Asset inspection and maintenance is an important activity. The ability to use data from...

Special Interest Group in Environmental Risk Management

The Industrial Mathematics Special Interest Group (SIG) in Environmental Risk Management was...

Green buildings-thermal comfort and energy consumption

The objective of this Industrial CASE PhD project is to develop novel mathematical models for...