PROCESS SYSTEMS ENGINEERING AND PROCESS SAFETY Chinese Journal of Chemical Engineering, 21(4) 376—381 (2013) DOI: 10.1016/S1004-9541(13)60467-X
Synthesis of Multi-component Mass-exchange Networks* LIU Linlin (刘琳琳)1, DU Jian (都健)1,**, Mahmoud M. El-Halwagi2,3, José María Ponce-Ortega4 and YAO Pingjing (姚平经)1 1
Chemical Engineering Department, Dalian University of Technology, Dalian 116024, China Chemical Engineering Department, Texas A&M University, College Station, TX 77843, USA Adjunct Faculty at the Chemical and Materials Engineering Department, King Abdulaziz University, Jeddah, Saudi Arabia 4 Chemical Engineering Department, Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Mich 58060, Mexico 2 3
Abstract This paper presents a superstructure-based formulation for the synthesis of mass-exchange networks (MENs) considering multiple components. The superstructure is simplified by directly using the mass separation agents (MSA) from their sources, and therefore the automatic synthesis of the multi-component system involved in the MENs can be achieved without choosing a ‘key-component’ either for the whole process or the mass exchangers. A mathematical model is proposed to carry out the optimization process. The concentrations, flow rates, matches and unit operation displayed in the obtained network constitute the exact representation of the mass exchange process in terms of all species in the system. An example is used to illustrate and demonstrate the application of the proposed method. Keywords multiple components, mass exchange network, superstructure, optimization
1
INTRODUCTION
The needs for pollution prevention, environmental conservation and sustainable development are among the primary obligations of modern industries. Design techniques have been developed to address the objectives of pollution prevention and sustainability. In this context, an important design problem is the synthesis of mass-exchange networks introduced by El-Halwagi and Manousiouthakis [1]. Several techniques and categories have been developed to address this problem including mass pinch technology [1-3], mathematical programming approaches [4-7], combined mass- and heat- exchange networks [8-11], batch MEN synthesis [12], and property-based integration [13, 14]. The problem of synthesizing MENs with multiple components is an important problem that has not received sufficient attention. In a multi-component system, the phase-equilibrium relationships between solutes and solvents yield complex relationships with components of the MEN. A common approach for simplifying the multi-component problem is to select a “key component” representing the dominant or the most demanding species for which the targets of the other species are compatible [1]. El-Halwagi and Manousiouthakis [15] developed a multi-component superstructure that uses mass-exchange design equations to reconcile the transfer of the multiple components. Papalexandri et al. [16] solved the MEN problem with a hyperstructure-based representation, where the simultaneous minimization of the operating and capital costs was investigated. Alva-Argaez et al. [17] presented a utility targeting approach assuming the same mass transfer capacity for solutes and solvents.
The incompatible multi-component problem is discussed by Wang et al [18]. In their work, the MEN is synthesized by determining the minimum mass separation agent mass flow rate requirement. Chen and Hung [19] presented a MEN synthesis approach based on a superstructure in which the maximum-sized exchangers are selected based on the demands of different components. Nonetheless, the compatibility for all components in the network was not guaranteed. Du et al. [20] also proposed a superstructure-based method, and the composition difference is especially considered in their research. The objective of this work is to overcome some limitations of previous research efforts. The paper presents a superstructure-based approach to determine the optimal solution for the synthesis of multi-component MENs. Compatibility of transferring different components is addressed and a large-enough search space is generated to embed potential configurations of interest without complicating the mathematical formulation. 2
PROBLEM STATEMENT
The MEN synthesis problem addressed in this paper can be stated as follows. Given are a set of rich streams with the concentration targets, a set of MSAs with the flowrate and concentration limits, and a set of transferable components that need to be removed from the rich streams to the MSAs. In addition, the phase equilibrium relationships of the considered components are also given. The corresponding constraint for mass-transfer feasibility is yi , j , p ≥ yi*, j , p , where
Received 2012-08-01, accepted 2012-10-29. * Supported by the National Natural Science Foundation of China (20976022). ** To whom correspondence should be addressed. E-mail:
[email protected]
Chin. J. Chem. Eng., Vol. 21, No. 4, Apirl 2013
yi*, j , p = f ( xi , j , p ) is the equilibrium function for transferring component p from rich-stream i to lean-stream j. The assumption that mass is transferred independently among the species is considered. The objective is to synthesize a MEN that can simultaneously meet the mass transfer requirements and the process constrains of all considered components at the minimum total annual cost (TAC). 3 3.1
THE SYNTHESIS METHOD Method statement
The traditional superstructure represents the alternative network structures by means of providing split streams for all potential rich-MSA matches during network stages. Then, the splitters and the mixers that link the streams among stages are required. To avoid the massive iterative calculation between stages and components, the superstructure representation is simplified in this paper. The MSA streams are directly used from their resources, and their splitters and mixers between stages are moved away. The conceptual graph of a two-stage superstructure that involves two rich streams (R1, R2) and two MSAs (S1, S2) is shown in Fig. 1. In Fig. 1, eight counter-current mass exchangers are set for the four potential pairs of stream matches in each stage. YiS,p and Yi T,p indicate the overall inlet and outlet concentrations for the rich steams, respectively. Each rich stream in the figure is split into two “substreams” for its possible matches with S1 and S2 during Stages 1 and 2. The mixing and splitting operations for the riches are settled between the two stages.
Figure 1
377
Other from the riches, the MSA-split streams all directly source from their resources with the overall inlet concentration of X Sj , p . The overall outlet concentrations X Tj , p are received by mixing corresponding MSA-substream at the left terminal of Stage 1. In the multi-component system, there are two or more solutes removed simultaneously from the MSAs to the rich streams. Therefore, several rules are given to guarantee the feasibility for all species in the obtained network. Rule 1: One component at least is transferred in each mass exchanger; Rule 2: For each mass-transferred component, its phase equilibrium constraints must be satisfied at the both ends of each mass exchanger. Rule 3: For the no-mass-transferred component, concentrations of inlet and outlet are equal, while meeting the phase equilibrium constraints. As has been stated, no “key component” is required in the proposed method. The mass exchanges of all involved components are considered simultaneously. It is suitable to either the compatible problem or the incompatible problem, and the network solution will appear to be a real reflection of the mass exchange of all components that are referred. 3.2
Model formulation
A non-linear programming (NLP) model is established for the optimization of the minimum TAC solution. The existence of the mass exchange unit is determined by the mass flow rates of streams in the model. “Zero” means no corresponding units, and “nonzero” is on contrast. The process simulation is carried out also based on the continuous variables of the tray numbers and the stream flow rates. The detailed mathematical
The simplified superstructure
378
Chin. J. Chem. Eng., Vol. 21, No. 4, Apirl 2013
model is represented below. (1) The overall mass balances for all components in the overall MEN In MENs, the transferable components are removed from the riches to the MSAs. For rich streams, the overall mass balances are set to ensure the sufficient mass exchange against their specified requirement. For MSA streams, they are used to guarantee the enough capacity of receiving the transferred mass. Notice that the balance for each MSA involves a constant inlet concentration X Sj , p since the sub-MSAs come directly from the resource.
(
)
Gi (Yi ,1, p − Yi , N S +1, p ) = ∑∑ gi , j ,k yiin, j , k , p − yiout , j ,k , p , j
k
∀i ∈ N R ,
(
)
∀p ∈ N P
(1)
(
)
S L j X Tj , p − X Sj , p = ∑∑ li , j ,k xiout , j ,k , p − X j , p , i
k
∀j ∈ N L ,
∀p ∈ N P
(2)
(2) The mass flow rate balances for MSAs in the overall network MSAs are used directly from their resources, and then discharged directly to their destinations. Each sub-MSA stream flows through only one unit in this study. So the total mass flow rate requirement of the jth MSA should be represented as the sum of that assigned to the i − j ( ∀i ∈ N R ) match in the kth ( ∀k ∈ NS ) stage: L j = ∑∑ li , j ,k , i
∀j ∈ N L
(3)
k
(3) The mass balance of rich streams during each stage Rich streams flow through all stages proposed in the superstructure. Among the different stages, the total flow rates of each rich stream are the same as Gi ( ∀i ∈ N R ). Gi is equal to the sum of its tributaries during each stage as Eq. (4) indicates:
Gi = ∑ gi , j ,k ,
∀i ∈ N R ,
∀k ∈ NS
In addition to this, the mass balances of the rich streams by stage are also necessary.
(
)
Gi (Yi , k , p − Yi ,k +1, p ) = ∑ gi , j ,k yiin, j ,k , p − yiout , j ,k , p , ∀i ∈ N R ,
∀k ∈ NS ,
)
∀p ∈ N P
(
(5)
)
out S gi , j ,k yiin, j , k , p − yiout , j , k , p = li , j , k xi , j , k , p − X j , p ,
∀i ∈ N R ,
∀j ∈ N L ,
∀k ∈ NS ,
after the last stage. The given inlet concentration of rich stream, YiS,p , is assigned as the inlet concentration of the superstructure, Yi ,1, p , in the model. Then, the last stage composition Yi , NS +1, p is the overall outlet composition of stream Yi T,p . For each exchange unit, the inlet concentration of rich stream corresponds to its stage composition as stated in the following: Yi ,1, p = YiS, p ,
∀p ∈ N P
∀i ∈ N R ,
Yi , NS +1, p = Yi T, p , Yi ,k , p = yiin, j ,k , p ,
(5) The assignment of rich stream concentrations along the network
(7)
∀p ∈ N P
∀i ∈ N R ,
(8)
∀j ∈ N L ,
∀p ∈ N P
(9)
(6) The equipment sizing equations for mass exchangers Tray-column is the most commonly used type of mass exchanger unit in practice. To obtain the stream outlet concentration for an exchanger whose tray number is already given, the Kremser equation listed below are applied in condition of the operating and equilibrium lines are straight: if mi , j , p gi , j ,k / li , j , p ≠ 1 ,
( (
S ⎡⎛ y in ⎡ li , j , p ⎤ i , j ,k , p − f X j , p ⎢⎜ Ni , j , k ln ⎢ ⎥ = ln ⎢⎜ yiout − f X Sj , p ⎢⎣ mi , j , p gi , j , k ⎥⎦ ⎣⎝ , j , k , p
⎛ mi , j , p gi , j ,k ⎜⎜1 − li , j , p ⎝ ∀i ∈ N R ,
∀j ∈ N L ,
) ⎞⎟ × ) ⎟⎠
⎞ mi , j , p gi , j ,k ⎤ ⎥ ⎟⎟ + ⎥ l , , i j p ⎠ ⎦
∀k ∈ NS ,
∀p ∈ N P
(10)
if mi , j , p gi , j ,k / li , j , p = 1 , yiin, j ,k , p − yiout , j ,k , p
(
S yiout , j ,k , p − f X j , p
)
,
∀i ∈ N R ,
∀j ∈ N L ,
∀k ∈ NS , ∀p ∈ N P (11) (7) The feasibility constrains for mass exchange It is necessary to guarantee that mass is removed from the rich stream to the MSA stream in a mass exchanger. Therefore, the constraints for decreasing the composition are required: yiin, j ,k , p ≥ yiout , j ,k , p ,
(6)
∀p ∈ N P
∀i ∈ N R ,
∀k ∈ NS ,
Ni , j , k =
(4) The mass balance of components in each mass exchange unit Mass exchange occurs in each mass exchange unit in the MEN. The mass removed from the rich stream must be equal to that absorbed by the MSA:
(
tions before the stage, and Yi , NS +1, p is the one that is
(4)
j
j
NS stages are involved in a MEN. In this model, Yi ,1, p , Yi ,2, p , " , Yi , NS , p are set as the stage composi-
∀i ∈ N R ,
∀k ∈ NS , S xiout , j ,k , p ≥ X j , p ,
∀j ∈ N L ,
∀p ∈ N P
∀i ∈ N R ,
(12)
∀j ∈ N L ,
379
Chin. J. Chem. Eng., Vol. 21, No. 4, Apirl 2013
∀k ∈ NS ,
∀p ∈ N P
(13)
(8) The feasibility constrains based on equilibrium relationships Mass exchange happens when the composition of the rich stream is no less than the equilibrium composition of the MSA. Hence, the following equations are set to prevent the undesirable mass transfer from the MSAs to the rich streams:
(
)
yiin, j ,k , p ≥ f xiout , j ,k , p , ∀k ∈ NS ,
(
)
S yiout , j ,k , p ≥ f X j , p ,
∀i ∈ N R , ∀p ∈ N P ∀i ∈ N R ,
∀k ∈ NS ,
∀j ∈ N L , (14)
∀j ∈ N L ,
∀p ∈ N P
(15)
(9) Limits for variables The bounds on the mass flow rates of MSAs and the overall final concentrations of streams are set in accordance with the process requirements or environmental regulations: L j ≤ Lupj ,
∀j ∈ N L
CASE STUDY
An example problem on the sweetening of coke-oven gas (COG) is taken from El-Halwagi and Manousiouthakis [1]. It has been recently investigated by Wang et al. [18], Chen and Hung [19] and Du et al [20]. The basic objective of the process is to remove H2S and CO2 from COG. Two rich streams (R1, R2) and two MSA stream (S1, S2) are considered in the example. The data taken from the literature [1, 16] are listed in Tables 1 and 2, in which the supply and target concentrations of components, the mass flow rates of rich streams and the flow rate upper bounds of MSAs are all given. The overall outlet compositions of H2S and CO2 should not be higher than their target compositions in streams. The operation time is set as 8150 h·a−1. The tray column with the annualized cost of 4552 USD·stage−1·a−1 is applied as the mass exchange unit. During the range of concentration involved, the phase equilibrium relationship of H2S and CO2 in S1 and S2 are represented by the Eqs. (20)-(23), respectively:
(16)
yH2S = 1.45 xHS12S
(20)
Yi T, p ≤ Yi ,upp ,
∀i ∈ N R ,
∀p ∈ N P
(17)
yH2S = 0.26 xHS22S
(21)
X Tj , p ≤ X up j, p ,
∀j ∈ N L ,
∀p ∈ N P
(18)
S1 yCO2 = 0.35 xCO 2
(22)
S2 yCO2 = 0.58 xCO 2
(23)
(10) The objective function The objective function of this model is set as the minimum total annual cost of the MEN: min TAC =
∑ ∑ ∑
∀i∈N R ∀j∈N L ∀k∈NS
HY
Citray ,j
× Ni , j ,k +
× Lj ) ∑ ( C MSA j
This MEN synthesis problem is solved with the proposed method by setting two superstructure stages, and the optimal overall mass-exchange process is represented by Fig. 2.
(19) Table 1
∀j∈N L
where the first term is the annualized capital cost from the column tray used in the mass exchangers, and the latter is the operating cost of external MSAs. 3.3
4
Model solution
The NLP model proposed in this paper features strong nonlinearity since it brings in several nonlinear equations. To effectively solve the model, a hybrid genetic algorithm-simulated annealing algorithm (GA-SA) [21] is adopted. This algorithm works commendably at the aspects of global searching and local optimum breaking in solving large-scale and serious non-convex problems. Following the concept of the method, the substream flow rates and the tray numbers of exchange units are selected as the genes in each solution ‘DNA’ of the algorithm. All the genes are real-encoded. The operations of selection, crossover and mutation are especially provided to accelerate the optimization process. In this work, the overall expression of the model and the optimization based on GA-SA are carried out with the programming language C++.
Gi/kg·s−1
Data for the rich streams [1] YHS2S
YHT2S
S YCO 2
T YCO 2
/kg·kg−1
/ kg·kg−1
/ kg·kg−1
/ kg·kg−1
R1
0.9
0.07
0.0003
0.06
0.005
R2
0.1
0.051
0.0001
0.115
0.005
Table 2 Lj,max /kg·s−1
Data for the MSA streams [1,16]
X HS 2S /kg·kg
X HT2S −1
/ kg·kg
S X CO 2 −1
/ kg·kg
T X CO 2 −1
/ kg·kg
①
−1
Ann. cost /USD·kg−1
S1
2.3
0.0006
0.031
0
0.171
0.004
S2
∞
0.0002
0.0035
0
0.103
0.006
① Papalexandri et al. [16] used an operating time of 8150 h·a-1 but erroneously multiplied the already annualized operating cost by an annualization factor (β = 0.2 a−1). To ensure consistent comparison with Papalexandri et al. [16] and Wang et al. [18], Chen and Hung [19] and Du et al. [20], we have used the same reduced operating costs used in these works.
As shown in Fig. 2, the total flow rate of S1 is 2.192 kg·s−1, and 0.37 kg·s−1 for S2, and four mass exchange units are involved in the resulting MEN.
380
Chin. J. Chem. Eng., Vol. 21, No. 4, Apirl 2013
Figure 2 The minimum TAC featured MEN
The compositions along streams are the real values, which matches exactly with the referred mass exchange unit. These concentrations indicate that the mass removal of H2S and CO2 both exist between the rich streams and the MSAs in units U1, U2 and U4, while in U3, only mass transfer of H2S occurs. This indicates that the simultaneous consideration of different species during the synthesis work is indeed necessary, because the network required by different components is distinct. Furthermore, as the mass removal degree among the species is not usually the same, units used to remove specific compounds are required at the end of the network to meet the mass exchange target, in this case unit U3 is used for this purpose. The TAC of the obtained network is 413450 USD·a−1, where 91040 USD·a−1 corresponds to the cost of the 20 column trays, and 322410 USD·a−1 is for the consumption of MSAs. This cost solution is very close to that in Du et al. [20], and lower than Wang et al. [18], Chen and Hung [19]. It is worth noting that as the solution space of this simplified method is smaller compared with the regular superstructure-based method, the computational complexity is significantly reduced. 5
CONCLUSIONS
A simplified superstructure-based representation has been proposed for the multi-component MEN synthesis problem. The MSAs sourcing directly from the resources are split into several substreams for the matches with the rich streams during each superstructure stage. No ‘key-component’ is required in this work. Each component has been considered with its proper thermodynamic and mass-exchange considerations. To meet the minimum TAC target, a NLP program has been formulated. The optimization/synthesis work is restricted by the process and the environmental regulations from all components. Finally, the applicability for the proposed approach is shown through the solution of an example problem previously reported.
Moreover, the incentive of considering all components simultaneously has also been demonstrated. NOMENCLATURE Ctray
annual unit cost of each column stage, USD·a−1
C MSA j
annual unit cost of the jth MSA/lean stream, USD·kg−1
f(x)
phase equilibrium relationship
Gi
total mass flow rate of rich stream i, kg·s−1
gi,j,k
mass flow rate of rich stream i that matches with MSA/lean stream j in stage k, kg·s−1
HY
operating hours per year, h·a−1
Lj
total mass flow rate of MSA/lean stream j, kg·s−1
Lup j
upper bound mass flow rate of MSA/lean stream j, kg·s−1
li,j,k
mass flow rate of MSA/lean stream j that matches with rich stream i in stage k, kg·s−1
mi,j,p
solubility equilibrium coefficient for component p between match ( i, j)
Ni,j,k
tray number of the mass exchange unit (i, j, k)
NL
{j|j is a MSA/lean stream, j = 1, 2, " }
NP
{p|p is a component, p = 1, 2, " }
NR
{i|i is a rich stream, i = 1, 2, " }
NS
{k|k is a superstructure stage, k = 1, 2, " }
X Sj , p
initial concentration of p in MSA stream j, kg·kg−1
X Tj , p X up j, p xiout , j ,k , p
overall outlet concentration of p in MSA stream j, kg·kg−1 upper bound concentration of p in MSA stream j, kg·kg−1 outlet concentration of p for the MSA/lean stream j that matches with rich i in stage k, kg·kg−1
Yi,k,p
concentration of p in rich stream i before/after superstructure stage k, kg·kg−1
YiS,p
initial concentration of p in rich stream i, kg·kg−1
Yi T,p Yi up ,p yiin,j ,k , p
overall outlet concentration of p in rich stream i, kg·kg−1 upper bound concentration of p in MSA stream j, kg·kg−1 inlet concentration of p for the rich i that matches with MSA/ lean stream j in stage k, kg·kg−1
yiout , j ,k , p
outlet concentration of p for the rich i that matches with MSA/ lean stream j in stage k, kg·kg−1
Chin. J. Chem. Eng., Vol. 21, No. 4, Apirl 2013
Superscripts in out S T up *
inlet outlet supply target upper bound equilibrium concentration
Subscripts i j k p
rich stream, i = 1, 2, " MSA/lean stream, j = 1, 2, " superstructure stage, k = 1, 2, " components, p = 1, 2, "
REFERENCES 1 2 3 4 5 6 7 8 9
El-Halwagi, M.M., Manousiouthakis, V., “Synthesis of mass exchange networks”, AIChE J., 35 (18), 1233-1244 (1989). Hallale, N., Fraser, D.M., “Supertargeting for mass exchange networks Part I: Targeting and design techniques”, Chem. Eng. Res. Des., 78 (A2), 202-207 (2000). Hallale, N., Fraser, D.M., “Super targeting for mass exchange networks Part II: Applications”, Chem. Eng. Res. Des., 78 (A2), 208-216 (2000). El-Halwagi, M.M., Manousiouthakis, V., “Automatic synthesis of mass exchange networks with single-component targets”, Chem. Eng. Sci., 45 (9), 2813-2831 (1990). El-Halwagi, M.M., Manousiouthakis, V., “simultaneous synthesis of mass exchange and regeneration networks”, AIChE J., 36 (8), 1209-1219 (1990). El-Halwagi, M.M., Hamad, A.A., Garrison, G.W., “Synthesis of waste interception and allocation networks”, AIChE J., 42 (11), 3087-3101 (1996). Gabriel, F.B., El-Halwagi, M.M., “Simultaneous synthesis of waste interception and material reuse network: Problem reformulation for global optimization”, Env. Prog., 24 (2), 171-180 (2005). Srinivas, B.K., El-Halwagi, M.M., “Synthesis of combined heat and reactive mass exchange networks”, Chem. Eng. Sci., 49 (13), 2059-2074 (1994). Isafiade, A.J., Fraser, D.M., “Optimization of combined heat and
10 11
12 13
14
15 16 17
18 19 20 21
381
mass exchanger networks using pinch technology”, Asia-Pacif. J. Chem. Eng., 2 (6), 554-565 (2007). Isafiade, A.J., Fraser, D.M., “Interval based MINLP superstructure synthesis of combined heat and mass exchanger networks”, Chem. Eng. Res. Des., 87 (11), 1536-1542 (2009). Du, J., Li, X.F., Chen, L., Yao, P.J., “Simultaneous synthesis of combined mass and heat exchanger networks”, In: 5th International Symposium on Design, Operation and Control of Chemical Processes (PSE Asia 2010), Singapore (2010). Foo, C.Y., Manan, A.A., Yunus, R.M., Aziz, R.A., “Synthesis of mass exchange network for batch process—Part 1: Utility targeting”, Chem. Eng. Sci., 59 (5), 1009-1026 (2004). Ponce-Ortega, J.M., Hortua, A.C., El-Halwagi, M.M., JiménezGutiérrez, A, “A property-based optimization of direct recycle networks and wastewater treatment processes”, AIChE J., 55 (9), 2329-2344 (2009). Nápoles-Rivera, F., Ponce-Ortega, J.M., El-Halwagi, M.M., Jiménez-Gutiérrez, A., “Global optimization of mass and property integration networks with in-plant property interceptors”, Chem. Eng. Sci., 65 (15), 4363-4377 (2010). El-Halwagi, M.M., Manousiouthakis, V., “Design and analysis of mass exchange networks with multicomponent targets”, AIChE Ann. Meet., San Francisco, November (1989). Papalexandri, K.P., Pistikopoulos, E.N., Floudas, A., “Mass exchange networks for waste minimization: A simultaneous approach”, Chem. Eng. Res. Des. , 72 (A3), 279-294 (1994). Alva-Argaez, A., Vallianatos, A., Kokossis, A., “A multi-contaminant transshipment model for mass exchange networks and wastewater minimisation problems”, Comput. Chem. Eng., 23 (10), 1439-1453 (1999). Wang, J.F., Shen, J.Z., Li, Y.R., Hu, S.Y., “Incompatible multicomponent mass exchange network synthesis”, J. Chem. Ind. Eng. (China), 55 (2), 297-304 (2004). (in Chinese) Chen, C.L., Hung, P.S., “Simultaneous synthesis of mass exchange networks for waste minimization”, Comp. Chem. Eng., 29 (7), 1561-1576 (2005). Du, J., Gao, Z.H., Chen, L., Yao, P.J., “Mass exchange network design using simultaneous optimization of composition differences”, J. Chem. Ind. Eng. (China), 58 (7), 1768-1775 (2007). (in Chinese) Xiao, W., Yao, P.J., Luo, X., Roetzel, W., “A new and efficient NLP formulation for synthesizing large scale multi-stream heat exchanger networks”, J. Chin. Inst. Chem. Eng., 37 (4), 383-394 (2006).