CE55504 – Surface Water Quality Modeling



CE5504 – Surface Water Quality Modeling

Lab 1 Handout: Numerical Methods - Applications of Excel and Visual Basic Software

Reading: Chapra, Surface Water Quality Modeling, Chapter 7 (intro, 7.1, 7.4)

Numerical Integration

The field of mathematical modeling was developed on the backs of those who could derive solutions for the ordinary differential equations that constituted the models. These analytical solutions have limitations, however, as described by Chapra (1997, Ch. 7):

• non-idealized loading functions: analytical solutions require the idealized loading functions (impulse, step, etc.) developed in Chapter 3. Actual loads behave in a much more arbitrary fashion can only be accommodated numerically.

• variable parameters: an analytical solution requires coefficient values remain constant over the course of the simulation. This is certainly not the case, as the specific growth rate for microbial populations, for example, may vary over the course of the simulation due to changes in environmental forcing conditions (e.g. light, temperature and nutrients).

• multi-segment systems: the development of 1, 2, and 3-D models for requires spatial segmentation that is inefficient to handle with analytical solutions (see Chapra 1997, p. 92, Eq. 5.21).

• non-linear kinetics: while linear (e.g. first order) kinetics are quite useful, there are several water quality applications that require non-linear kinetics (e.g. the Monod model). Analytical solutions cannot be obtained for these.

Thus we turn to numerical methods, developing numerical or computer-based solutions. The basic idea here is the solutions are trivially simply, but exhaustively long were it not for computing power.

The Euler Method

In the Euler method, we transform the differential equation (here the exponential growth model) as follows:

[pic]

[pic]

such that we can calculate a change in X and then add it to the original value of X to get an updated value at the end of the time interval dt:

[pic]

This approach is implemented in a computer program in the form of a loop, calculating X over an interval 0 to tmax:

[pic]

The selection of a value for dt is critical to the success of the process (Spain 1982, Figure 5.1). We find as dt(0, the numerical solution approaches the analytical solution. The trade-off becomes one of computation time – as we decrease the time step (dt) the solution becomes more accurate, but the number of computations that must be performed increases. Alternative numerical methods (Chapra 1997, Chapter 7: Heun’s Method, Runge-Kutta Methods) can provide accurate results at significantly longer time steps, thus reducing computation time.

A. Setting up the workbook.

1. Open a new workbook.

2. Label four worksheet, one for each of the Lab 1 problems.

3. Establish and label input-output regions.

4. Open the Visual Basic editor

Tools > Macro > Visual Basic editor

5. In the empty programming work space, type –

Sub Integration

and hit Return. Your program is intended to fit between these lines.

6. Returning to the worksheet, add a program control –

View > Toolbars > Forms

Select the ‘bulls-eye’ form and drag it onto the worksheet. Move away from it and then return to activate. You can ‘right-click’ on the form to edit the text, shape and color.

7. Right click on the form and select –

Assign macro …

and then select

Sheet1.Integration

and click on ‘OK’. This links the form and the program, running the software every time you click on the form.

8. The worksheet is now ready to have the software added.

B. Writing the program.

Sub Integration()

Assign model inputs (C0, k, tmax, dt)

variablename = Cells(row,column)

Initialize state variable

variablename = variablename0

Begin outer loop

For i = 1 to tmax

Assign model output

Cells (row,column) = variablename

Begin inner loop

For j = 1 to 1/dt

Perform incremental calculation

dC = k * C * dt

Perform update

C = C + dC

Close loop inner loop

Next j

Close loop outer loop

Next i

EndSub

C. Applying the program

We will apply the software to examine four population growth models:

1. Exponential Model

[pic]

where: µmax = 1 d-1, X0 = 2 mg∙L-1 and tmax = 10 days. Calculate x = f(t), comparing the analytical solution to the numerical solution for the Euler, Heun and 4th Order Runge Kutta methods to evaluate the accuracy of the integrators for different time steps (e.g. values of dt ranging from 1 day to 0.0001 days). Consider the impact of this on run time and accuracy. Examine model sensitivity to values of X0 and µmax.

2. Logistic Model

[pic]

where: Xmax = 25000 mg∙L-1, tmax = 25 days, dt is as determined above for best performance with an integrator of your choice. Examine model sensitivity to changes in Xmax. Would you expect carrying capacity or resource availability to govern population growth in most ecosystems and world populations.

3. Monod Model

[pic]

where: S0 = 500 mg∙L-1, Ks = 5 mg∙L-1, Y is 0.5 mgX∙mgS-1 and tmax = 10 days; all other inputs as before. Note that this requires equations for both S and X:

[pic]

The coupled differential equations are solved simultaneously, carrying the two equations through the increment (dX = ) and update (X = X + …) steps together.

The biomass curve looks similar to that for the logistic model; how does it most differ? Why? How do the half-saturation constant and yield coefficient influence the maximum value of X? Why? Is this realistic? What does the half-saturation constant influence? Why is this? There is one common shortfall of both the logisitic and Monod models. What is this?

4. Batch Growth Model

[pic]

where: kd = 0.2 d-1, tmax = 25 days and all other inputs are as before. The equation for S remains as before, unless recycle is added,

[pic] without recycle

[pic] with recycle

What is the effect of recycle on model behavior? Can you explain this? Which mechanism is most realistic? Compare model sensitivity to variation in the half-saturation constant for both cases. Can you explain the resulting behavior?

5. Completely-Mixed Flow Reactor

[pic]

[pic]

where: V = 100 L, Q = 10 L∙d-1, and Sin = 200 mg∙L-1. Examine model sensitivity to the half-saturation constant and initial and final substrate concentrations. Explain the behavior in each case.

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download