Ind. Eng. Chem. Res. 2007, 46, 1519-1526
1519
PROCESS DESIGN AND CONTROL Simulation and Control of an Upflow Anaerobic Sludge Blanket (UASB) Reactor Using an ADM1-Based Distributed Parameter Model S. J. Mu,† Y. Zeng,† B. Tartakovsky,‡ and P. Wu*,† Institute of High Performance Computing, 1 Science Park Road, #01-01, The Capricorn, Singapore 117528, and Biotechnology Research Institute, NRC, 6100 Royalmount AVenue, Montre´ al, Que´ bec, Canada H4P 2R2
In this study, an ADM1-based distributed parameter model was used to evaluate the influence of organic loading and external recirculation on upflow anaerobic sludge blanket (UASB) reactor performance. A process control strategy based on the manipulation of the external recirculation-to-influent ratio was proposed. This strategy was compared with the traditional dilution rate-based method to control the effluent chemical oxygen demand (COD) concentration and methane gas flow rate. Furthermore, a simple rule, which suggested when and how to use the dilution rate and the recirculation-to-influent ratio in a multivariable control system, was proposed. The simulations showed that the use of a multivariable control system improved reactor performance and stability, in comparison to a traditional dilution rate-based control method. 1. Introduction Worldwide, a large number of anaerobic treatment systems are based on the concept of upflow anaerobic sludge blanket (UASB) reactor design.1 Because of the high sensitivity of the anaerobic microbial consortium to variations in organic load, alkalinity, and other external disturbances, the anaerobic digestion process is often prone to population shifts and process instability. The existence of significant axial pH, substrate, and biomass gradients in UASB reactors2-4 is an additional source of potential process instability. Therefore, the use of advanced process control is desirable to improve process stability, and it has been an active research area since the concept of UASB reactors was first introduced in the 1970s. In anaerobic reactors, the influent flow rate or dilution rate (defined as the influent flow rate divided by the reactor volume, or the reciprocal of the hydraulic retention time (HRT)) is used as the main manipulated variable to control the effluent chemical oxygen demand (COD) or volatile fatty acid (VFA) concentration.5-7 The dilution rate can also be used to control biogas production8,9 and reactor pH or alkalinity consumption.10 In practice, the efficiency of the dilution-rate-based control is limited by upstream factory discharges7 and wastewater storage capacity.11 This work presents a new control strategy, which uses the external recirculation rate as a manipulated variable. Indeed, experimental studies have demonstrated the significant impact of recirculation flow rate on the stability and performance of anaerobic reactors, although the impact is dependent on the reactor type, operational conditions, and wastewater composition.12 Sam-soon et al.2 found that recirculation could increase the effective alkalinity of feed and thus reduce the alkalinity * To whom correspondence should be addressed. Tel.: (65) 64191212. Fax: (65) 64191290. E-mail address:
[email protected]. † Institute of High Performance Computing. ‡ Biotechnology Research Institute, NRC.
requirement for a UASB reactor. It was also demonstrated that a higher recirculation-to-influent ratio improved reactor performance at moderate organic loads. Also, it was determined that an introduction of recirculation minimized total operational costs, because of considerable savings in alkali addition.13 Setiadi et al.14 suggested that higher recirculation rates improve the stability of methanogenic populations and reduce alkalinity requirements. However, Mshandete et al.15 observed that, at high organic loading rates (OLRs), the reactor with a lower recirculation rate performed better. The anaerobic digestion process is often simulated by assuming conditions of ideal mixing. However, ideal mixing models are not suitable for studying the impact of recirculation on reactor performance. In this work, process control simulations were performed using a recently developed ADM1-based distributed parameter model, ADM1d.16 This model uses the kinetic equations of the IWA anaerobic digestion model No. 1 (ADM1),17 which has been widely used both for laboratoryscale and full-scale anaerobic reactors.18,19 ADM1d combines ADM1 kinetic equations with the axial dispersion model of a UASB reactor and enables prediction of the component distribution along the reactor height. ADM1d, in many respects, is similar to other emerging distributed-parameter models,20,21 but uses position-dependent dispersion coefficients22 and hyperbolic tangential functions to model sludge distribution and sludge washout in a UASB reactor. 2. Materials and Methods 2.1. Process Model. The ADM1d distributed parameter model16 was used for process simulations. A brief model description is given below. Liquid and gas-phase components of the model are those of ADM1.17 Each liquid-phase component ci is described by a material balance,
(
)
∂ci (z, t) ∂ci ∂ ∂ Di(z, t) - (uupi(z, t)ci(z, t)) + ri(z, t) ) ∂t ∂z ∂z ∂z (1)
10.1021/ie060853l CCC: $37.00 © 2007 American Chemical Society Published on Web 02/07/2007
1520
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
with Danckwerts boundary conditions,
Di(z, t)
∂ci ) uupi(z, t)(ci(z)0) - ci,in) ∂z ∂ci )0 ∂z
(at z ) 0)
(at z ) H)
(2) (3)
where ci,in is the influent concentration of the ith component and H is the reactor height. The axial dispersion coefficients of the components were almost linear, relative to the upflow liquid velocity and exponent negative to reactor height.22 The gas-phase components are described by the expression
∂sgas,j(z, t) 1 ∂ ) rT,j(z, t) (q (z, t)sgas,j(z, t)) ∂t Ar ∂z gas
(4)
where j is the index of biogas components (H2, CH4, and CO2), sgas denotes the concentration of the jth compound, Ar is the reactor cross section, and qgas is the volumetric flow rate of biogas, which is defined as
qgas(z, t) ) RT Patm - pgas,H2O
∫0V
(
rT,H2(z, t) 16
+
rT,CH4(z, t) 64
)
+ rT,CO2 dV (5)
where Patm and pgas,H2O denote the total gas pressure and the partial pressure of water vapor, respectively; V is the reactor volume; T is the operating temperature; and R is the gas constant. The boundary conditions of the gas components are sgas,j ) 0 at z ) 0 and dsgas,j/dt ) 0 at z ) 1. The overall mass balances of the reactor are
qin + qrec ) qs + Q
(6)
qinci,in + qrecci,r ) qsci + Qci,0
(7)
where qin, qrec, and qs are flow rates of the influent, recirculation, and bypass streams, respectively; Q is the total input flow; and ci,in, ci,r, and ci,0 are the concentrations of the ith component in the corresponding flows. The sludge distribution in the reactor is modeled by a hyperbolic tangent function, which defines a maximal attainable biomass concentration (xat(z) ) ∑ xi ) at each reactor position:
xat(z) )
(
)
xmax - xmin [1 - tanh(β(h - Hs))] + xmin (8) 2
where xmax and xmin are the total attainable concentrations of microorganisms in the sludge bed and the liquid above the sludge bed, respectively; z is the position within the reactor; Hs is the sludge bed height; and β is a positive constant. The retention of biomass granules in the sludge bed is described by multiplying the transport term (u) in eq 1 by a washout factor (Rw), modeled as
Rw(z) )
tanh(xtotal(z) - xat(z)) + 1 2
(9)
where xtotal(z) ) ∑ xi(z) is the total biomass concentration at reactor position z, xat(z) is defined by eq 11, and xi(z) is the concentration of the ith microorganism. Equation 9 implies that, if the total biomass concentration exceeds the maximal attainable biomass concentration, then Rw ) 1, and the transport term in eq 1 is the same as for all other liquid-phase components.
Otherwise, Rw ) 0 and the transport term in eq 1 will be equal to zero (i.e., the biomass is retained in the reactor). This approach simulates biomass accumulation within the sludge bed and its washout in the liquid zone, whereas a transition from the sludge bed to the liquid zone is described by a monotonous function. The model was validated using experimental results obtained in a 10.4-L laboratory UASB reactor that was equipped with sensors for on-line COD and VFA measurements at several distances from the reactor base.4 Most of the ADM1d parameters were the same as those in the ADM1 model, whereas some parameters were re-estimated. The re-estimated parameters included the maximum specific uptake rates of butyrate/valerate, propionate, and acetate (km,c4 ) 2.3 d-1, km,pro ) 2.4 d-1, and km,ac ) 1.0 d-1, respectively) and the gas-liquid transfer coefficient (kLa ) 120 d-1). 2.2. Simulation of pH Control. To simulate the overall process control of a UASB reactor, a simulation of the pH control loop aimed at maintaining a constant pH within the reactor was required. As in ADM1, the pH value in ADM1d is calculated using six physicochemical processes, which describe the acid/base equilibria of CO2/HCO3+, NH4+/NH3, acetic acid/ acetate, propionic acid/propionate, butyric acid/butyrate, and valeric acid/valerate. The dynamic equation of the H+ concentration is shown below and can be found in the works of Batstone et al.17 and Rosen and Jeppsson:23
dSH+ θ 1 ) - + xθ2 + 4KW - SH+ dt 2 2
(10)
where KW is the equilibrium coefficient of H2O and θ is an ion balance variable:
θ ) Scat+ + (SIN - SNH3) - Shco-3 -
Sac- Spr- Sbu64 112 160 Sva- San- (11) 208
Both cations and anions were included (Scat+ and San-, respectively), to simulate the influence of strong bases and acids, respectively, in feed streams.17 In this study, only additional alkali (sodium hydroxide, NaOH) was considered to be present in the influent stream. In pH control simulations using ADM1d, the pH value at the reactor inlet is used as a controlled variable and the flow rate of additional sodium hydroxide is the manipulated variable. For the controller, the input is pH-calculated using the model and the output is the desired alkali flow rate. The calculation of the desired alkali flow rate is based on the desired concentration of OH- or Scat+ (eq 11) in the reactor. The additional influent alkali flow rate is calculated from the desired Scat+ value, using the following equation:
(
Falk ) Fin
)
Calk,reactor - Calk,COD Calk,in - Calk,reactor
(12)
where Fin is the influent wastewater flow rate; Calk,reactor is the desired controlled reactor alkalinity level Scat+; and Calk,COD and Calk,in are the alkalinity levels of influent wastewater and additional alkalinity flow, respectively. In the following simulations of control strategies, pH control is performed by manipulating the flow rate of a 0.5 N NaOH solution with the parameters shown in Table 1. 2.3. Numerical Methods. ADM1d was solved by a finitedifference (FD) method24 with an error of order h2. The central
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1521 Table 1. Parameters Used for UASB Reactor Control Simulation parameter
R-based COD control
D-based COD control
pH control
CV setpoint Kp Ki MV boundary DV boundary
0.2 kg COD/m3 0.05 0.001 [0.5 15] [1.67 3.75]
0.2 kg COD/m3 8.0 × 10-5 1.0 × 10-5 [1 8] [1.67 3.75]
7 0.02 [0 0.02]
finite difference approximation was used for the internal finite difference points, and the forward and the backward approximation was used for the first and second boundary conditions, respectively. By specifying N internal finite difference points, i.e., N + 1 equally spaced intervals, this method resulted in 36 × N ordinary differential equations (ODEs) and 36 × 2 algebraic equations. The set of ODEs was solved numerically by the integration function “ode15s” of MATLAB with relative and absolute tolerance parameters (“RelTol” and “AbsTol”) set to 10-6 and 10-8, respectively. 3. Results and Discussion 3.1. Simulation of Organic Overloads and Underloads. The study was started by ADM1d simulations showing the influence of the OLR and the recirculation-to-influent ratio (R) on UASB reactor performance. The UASB reactor operation was simulated for a period of 320 d. Under normal operating conditions, the following values were observed: OLR ) 6 kg m-3 d-1, influent flow rate ) 0.024 m3/d, and R ) 4. During the simulation, the influent wastewater flow rate and R remained unchanged, and no additional NaOH was fed into the system while the OLR was varied (no pH control). After the first 40 d of operation, the OLR was reduced to 4 kg m-3 d-1 and maintained for 40 d. Several step increases of the OLR then were imposed until an OLR of 16 kg m-3 d-1 was reached at day 240. This OLR value was maintained for another 40 d. At day 280, the OLR was decreased to a normal value of 6 kg m-3 d-1 and maintained for an additional 40 d (see Figure 1a). Figure 1 shows the influence of OLR variations on key process variables. During the first 240 d, the OLR increased but remained within the loading capacity of the reactor (OLR < 16 kg m-3 d-1). The effluent COD was maintained within an acceptable range (
