Gas desorption from liquids

GAS DESORPTION

FROM LIQUIDS

W. PASIUK-BRONIKOWSKA and K. J. RUDZIhSKI Instituteof PhysicalChemistryof the Polish Academyof Sciences, 01-224Warszawa,Poland (Received10 February 1980; accepted 20 Nooember 19i30) Ah&act-The formal methodologyof modelingthe simultaneousquiescentand bubbledesorptionof gases from liquidsis presentedby settingup a learningmodelof single gas desorptiondue to chemicalreaction.It is based on the deterministicpopulationbalance of noncoalescingand nonbreakablesphericalbubbles in a perfectly mixed volume of liquid. Simple approximateformulasare used to describe the rates of componentprocesses, i.e. of quiescent desorption and nucleation, growth and escape of bubbles. The method and exemplary results of estimatingthe model parameters from experimental data are shown, which are restrictedto the size independent growthof bubbles. _

INTRODUCTTON

Gas desorption from liquids is a system action to destroy gas supersaturation which is created within the liquid phase by some other process, e.g. chemical reaction, depressurization, changes of temperature, cavitation or salting out. It is generally recognized that desorption follows two distinct mechanisms. One of them is phase transition of dissolved gas at the liquid free surface (quiescent desorption), while the other is nucleation and growth of gas bubbles (bubble desorption). In most practical cases either the two mechanisms are employed simultaneously, or the quiescent one is employed alone. Danckwerts [1] neatly explained that bubble desorption is very different from absorption (and hence from quiescent desorption), “in which the area of surface available [for mass transfer] is determined by external factors and not by the (. . .) process itself”. This paper is basically concerned with gns desorpfion due to chemical reaction. which is quite common in industrial practice. Exemplary cases are: CO2 desorption in fermentation-process study and control[2] and flotation of solids[3], CO2 and Nz desorption in microbial denitrification of waste waters[4], N2 desorption in complex reactions of atmospheric pollution control and its use for agitation and transport of liquids or suspensions [5], coal pyrolysis[6]. Nevertheless, further extension of presented methodology is possible SO, as to cover other significant cases: desolption due to chemical and depressurization (iron and steel reaction manufacture[7-91, refining of glass[lO], regeneration of spent absorbents in tlash drums or reboilers [I, 11, 121); desorption due to physical or chemical absorption (chemical technology[l21_ a review,[l3]), pollution control [14,15]); gas desorption due lo depressuriztltion (degassing of oils and fuels[l6,17], regeneration of spent of manufacture foamed absorbents[l, 121, elastomers [19,20], iron and steel manufacture [7,9,21], floatation in deep fermenters [22], decompression disease[23]). The selected cases are representative for simultaneous bubble and quiescent desorption (only few of them may employ the quiescent process alone).

Within most of them the unsteady state conditions are of major significance. No formal mathematical treatment was reported which covered both mechanisms of desorption simultaneously. Up to now desorption was studied in parts rather than as the entire process. Quiescent desorption was investigated with the use of gas absorption methods, i.e. film and penetration theories[l, 1I-131. These theories were also applied to develop criteria of supersaturation, and hence of bubble nucleation, in desorption with chemical reaction[ti] and in simultaneous desorption and physical[24,25] or chemical [25,26] absorption. However, neither of the theories can be used to describe the kinetics of gas desorption at bubbling conditions [ 1,12,26,27]. The experimental attempt of develop the supersaturation criin the terion was also reported for depressurization string-of-spheres and stirred vessel desorbers[28,29]. Bubble desorption was extensively studied as growth of single bubbles in supersaturated solutions ([7,20,21,3Ol_a review[31], and many other papers). Few very simplified[l6,18] or semi_empirical[l7] models are available, which, however, are not capable of simulating exactly the desorption. Recently a “sink model” was suggested for gas desorption into stationary arrays of bubbles in gfassmelts[32] (nucleation was not considered within this model). Modeling of the bubble gas desorption at unsteady-state conditions was also reported [33], which was based on the bubble population balance. The population balance approach was strongly recommended by Resnick and GaLOr[34] for general study of gas-liquid systems, but further applications considered time invariant bubble populations only. This paper extended the population balance approach over simultaneous bubble and quiescent desorption. To introduce the methodology a simple case was selectedsingle gas desorption due to chemical reaction, from well mixed liquid into one-component gaseous phase. For this case a learning model (in the sense of Shinnar[JS]) was developed, general restrictions of which are: neglected evaporation of solvent, perfectly mixed liquid phase,

1153

1154

W. PASNK-BRONIKOWSKA and K.

constant saturation concentration, negligible effects of capillary overpressure, noncoalescing and nonbreakable spherical bubbles, absence of bubble feed or outflow streams (i.e. nucleation and escape are allowed only), size independent growth of bubbles. KmMAL

FQUATtONs

The process of gas desorption can be described as the development of gas supersaturation in liquid, due to the difference between the rate of gas production in liquid and the rate of gas desorption. For a perfectly mixed volume of liquid we may write: dz/dt = r(z, t) - F(z, t)/ V(t) z(0) = z”.

N&z, f)U, +

I

kl.,

0

Ni,(z, f, L)d2$(L, r) dL.

L(fo)=O,

to E to.t1.

(4b) (4c)

To avoid confusion in eqns (3) and (4) the dependence of nucleation and growth rates on gas supersaturation r(t) was not marked. If the particular form of z(t) is necessary for the solution procedure, then eqns (4) have to be solved simultaneously with eqns (I) and (2). Once the bubble population density + has been determined, the parameters of a bubble dispersion can immediately be calculated, e.g. number, surface area and volume of bubbles per unit volume of liquid: nb =

=- rb_(L,1)dL I0 =ln,, ab = EL*$(L, t) dL I0 =,, lJ*= o (&)L3$(L, t) dL I

(la)

[F,(.Z,t)+ Fb(Z, r)l/W =

dLJdt = G(L,‘t)

(1)

When desorption is caused by a step change of some quantity, the rate of gas production should be withdrawn from eqn (I) and the initial condition should be properly adjusted. Gas supersaturation, z, is defined as the difference between actual (bulk) and saturation concentrations of the dissolved gas. Generally, time dependence of the latter should be specified. Within the present case this concentration is constant. We assume that the functional forms of r(z, t) and V(t) are defined and known. Rate of desorption is a sum of the quiescent and the bubble rates, eqn (2). Generally the bubble rate can be calculated only with use of the distribution function of bubble diameters, since it contains the bubble interfacial area implicitly. This is ‘so because the unit rate of bubble desorption may depend upon the bubble diameter. flv=

1. RUDZINSKI

(2)

(5) (6) (7)

where L max= L( r, to = 0). The integration variable in eqns (S)-(7) can be changed to more convenient to, if only dL/dtO is derivable from eqns (4). When the growth of bubbles is size independent, eqn (3) reduces to the previously reported form, which has a general explicit solution[33). Moreover, the unit rate of desorption Nb becomes size independent as well. It can be separated from the integrand in eqn (2) to introduce explicitly the area ab in there. NWCLEATION, GROWTHAND ESCAPEOF BUBBLES

In our previous paper[33] we suggested simple rate Unit rates of desorption, N4 and Nb, are discussed with expressions for the physicochemical processes which the bubble growth in the next section. The distribution $ accompany gas desorption. Given no comment on, these is recognized as solution to the bubble population expressions looked somehow arbitrary. Now we introduce them again as the suitable “first touch” apbalance. For the noncoalescing and nonbreakable spherical bubbles which can leave the liquid phase only proximations to the exact yet undefined rate functions. Nucleation of bubbles was extensively studied within by escape through free surfaces (rate of escape is disthe theory of boiling. Blander provided an excellent cussed in the next section), the balance is: review and “state of art” paper[37]. However, it should be realiid that nucleation in reacting solutions of elec01 ++G(L, t)v(t)$(L, t)] trolytes may require different treatment. To approach the problem we assume that the rate of nucleation is = - f(L 0 VuML, 0 differentiable function of gas supersaturation. This Oa) function can be expanded into Maclaurin series which is m, 0)= 0 later truncated at the first nonzero term: @(a,t) = N(NG(O,I). W x= BZZ. (8) Standard analytical treatment[36] converts eqn (3) into the simple parametrical system: Generally the constant parameter BZ may depend on properties of a liquid solution, e.g. on the ionic strength. Formula (8) is consistent with the diffusion controlled rate of nucleation which was derived by Blander, Hengstenberg and Katz[37,38] and applied by Attar[Ll to + G,(L{t’, to), f’)] dt’ viscous reacting solutions. 1 Growth of spherical bubbles was also investigated t=o +(L,t)=o within the theory of boiling. Many of the reported

$[V(t)Q(L,

1155

Gas deserption from liquids

developments are readily applicable to gas desorption from liquids, as shown by Tao[31]. We are not aware of any report on the growth of bubbles in a concentration field that varied independently with time. For simplicity and convenience we suggest to assume the size independent growth-G: = G(z), and to employ the truncated series expansion as for nucleation. Since in solutions to some related problems a square root of a driving force is observed[39,40], it may be reasonable to use a suitable substitution, which would insert noninteger powers into the series: l/2, 1, 3j2, 2, 512,. . . . For the ideal gas the bubble growth rate and the unit rate of bubble desorption differ by a constant factor only, so we may write: Nb =B,%+)

(9)

G = (2MIp)Nb = (2Mlp)B,X42).

(10)

In spite of differences in geometries and movements of gas-liquids interfaces we suggest the same treatment for the unit rate of quiescent desorption. Escape of bubbles from the liquid phase is the most obscure of the considered phenomena. In agitated liquids it probably is more of random than of deterministic nature. Being interested in the phenomenological description of desorption we leave out the details of the problem and we modify the bubble population density to express directly the rate of escape:

where B, = B,sT(~B,M/~)~B~ Bs = 2B,B, M/p. The parameters B,, Bd and B5 can be estimated by

simple experiments. Here we present main concepts of the estimation, while details are available it@]. Unsteady state experiments were performed at constant temperature (298 K) and pressure (760 mm Hg). A perfectly mixed, batch reactor (stirred vessel) was used, of constant liquid volume (0.5 dm3). Dissolved reactants (sulphamic acid and nitrite) were brought together in the reactor, and since that moment the rate of gas nitrogen desorption F(t) was measured along time with the automatic bubble soap flowmeter. All solutions were saturated with nitrogen before the experiment, so the evolution of bubbles was immediate. This observation is consistent with[28,29]. Time dependence of the rate of nitrogen production in chemical reaction, r(t), was determined by initial concentrations of reactants (the only variable parameters of the experiments), according

to the equations in the Appendix. Gas supersaturation in liquid, r(l), was calculated directly from the rates of gas production and desorption (these three quantities will be referred forth as “experimental”). When determined, gas supersaturation and rate of desorption were used for parameter

estimation.

The complete

set of necessary

experimental data was: F(t)/ V-time dependence of the rate of gas desorption t,,,-experimental time at which 57 V has a maximum z(+time dependence of gas supersaturation in liquid (dz/df),-time derivative of r(t) at t,,,.

[ ~~~f~~~~]=,~~~L,t)~~i,t)~~)

Then we introduce a rather deterministic approximation for the modifying function:

f(L, 0 =

B3L.

(12)

Hence the rate of escape is proportional to the number and to the diameter of bubbles. A systematic experimental study is required to prove whether the approximations are sufficiently exact. It should be taken into account that any complexity in the form of the rate expressions may bring about signiiicant ditliculties in the empirical estimation of parameters.

The usual form of the data is sets of discrete values. If the sets are too scarce to allow for numerical integration, they should be made dense by a suitable interpolation. This, of course, may result in reduced accuracy of estimation. To find the value of B5 we assume that the “free surface” component in eqn (13) (i.e. BId(z) has no significant influence on the localization of the maximal rate of desorption along the time axis. Hence the time

derivative of the remaining component should vanish at t, (existance condition for maxima). This results in the equation with one unknown variable, Bs:

JMlMATIONOFfAUAME~

When eqns (8H12) are inserted into the model and

one expression is used for unit rates of desorption (this easily relaxable assumption simplifies further calculations), the total rate of desorption in eqn (1) becomes: where L” =

eml&t’)) I ‘0

dt’.

The primitive “trials and errors” method is efficient enough to supply the positive root of eqn (14) with the

1156

W. PASIUK-BRONIKOWSKA and IL I.

five digits accuracy. Although no existence or uniqueness of the solution was proved, the numerical variation of the right hand member of eqn (14) rather negated the possibility of mare positive roots. Eventually other parameters (B, and Blaf) are estimated by fitting eqn (13) by least squares to the experimental data. If the liquid free surface is known (e.g. being the cross section of a vessel) the original parameters B,, Bz and & are readily evaluated. Note, that in contradiction to the earlier case [33], only one set of these parameters exists for one set of experimental data If, however, the weak restriction is released (i.e. two distinct functions are used to describe unit rates of desorption), then the “random” character of the model is immediately recovered.

RurxznLun

The results of experimental tests are presented on Figs. 1 and 2. The variable process parameter (&= modified directly the rates of gas production (r in 8,) Fii. 1, for eqns see Appendix). Consequently gas supcrsaturation (Fig. 2) and rate of desorption (R V in Fii. 1) are influenced according to model eqns (1) and (13). Figures 1 and 2 indicate satisfactory agreement of experimental data (0) and theoretical predictions (-) (note that ‘Wing” refers to rates of gas desorption which were calculated from experimental supenaturations, while “simulation” refers to numerical solutions of eqns (1) and (13). which are entirely based on the estimated parameters). Pi 3 shows simulated values of specific interfacial area in both experiments, and Fig. 4-the exemplary bubble size distribution (experiment I).

J

0 0

100

50 t,

5

0

50

Km

t,

*

simulatiotL Fig. I. Rates of nitrogen production (r) and desotption (fl V). 0, experiment; ---, fitting; Expelimen~ I: c$\p = coNp= 0.0105 mol/dm’ B, = 4.0026 x 10+, Bz= 10.806x lff, ES= 55.153. Experiment II: &, = coHp= O&W mol/dm’ B, = 1.9873 x 10-4, Bz = 7.1723 X 106, & = 28.238.

Fig. 2. Nitrogen supersaturation in liquid. 0, experiment; -,

simulation. Experiments I and II as in Fii. 1.

1157

08s desorption from liquids

mental effort is required to correlate the model parameters and process variables, the presented methodology can already be recommended as the convenient tool in the fundamental studies of industrial cases of gas desorptioo. Reaction of sulphamic acid with nitrite is an excellent model reaction for desorpthn. since it has flexible, acid catalysed kinetics [Appendix]. To give the final closure, we would like to point out a design convenience, which arises from the use of population balance technique. The hardly measurable quantities, like number, surface area or volume of bubbles, can readily be calculated from the bubble population density, which is derived from the measurable quantities like the rate of gas desorption or gas supersaturation in liquid.

10

7

I ; 5

4

cI

NQTATION

-

0

Ion

50

ab

t_* %

Fii 3. Simulation of gas-liquid in&facial

areafor experimentsI

and n

ai

Experimental d+erminatioa of both quantities is desirable to cross-check the model, but the measuriug tech-

nique is hardly available. Simulations were performed with the explicit fourth order Runge-Kutta algorithm and Simpson or trapezoidal integration. coPJcLuWNs In the previous sections the learning model of single gas desorption was set up, which considered chemical reaction, quiescetit desorption and nucleation, growth and escape of bubbles. Simple approximate formulas were used to express the rates of component processes, excluding the rate of chemical reaction which followed the experimental kinetic equation. Experimental tests showed that this treatment was efficient in modeling the relevant case of gas desorption. It was also concluded that within this case the quiescent &sorption could not be neglected (compare[33]). Although further experi-

interfacial area of a bubble dispersion, per unit volume of liquid, l/dm area of liquid free surface per unit volume of liquid. lldm constant in eqn (9), mol’~/(dm’n~ s) constant in eqn (S), l/(mol s) constant in eqn (12), I/(dm s) constant in eqn (13), dmgn/(mo13”. s’) constant in eqn (13), dm3’V(molu2- s3 molar concentration, molldm’ formal (analytical) concentration, mol/dm’ modifying function in eqn (ll), l/s rate of gas desorption, mol/s rate of bubble gas desorption, mol/s rate of quiescent gas desorption, mol/s growth rate of bubbles ( = dL./dt), dmls aG(L, t)/aL l/s ionic strength, mol/dm’ kinetic constant, dmY(molz .s) dissociation constant, mol/dm’ thermodynamic dissociation constant, mol/dm’ bubble diameter, dm maximal bubble diameter of the dispersion, dm molar mass of a gas, Wmol

/

100 t.

‘0 0

L.102, Fig.

3

2

1

4. SimuJatfon of

drn

the bubble populationdensityfor experimentI.

6

w.

11%

P~-B~NIKOWSKA

number of bubbles per unit volume of liquid, l/dm3

unit rate of gas desorption (i.e. rate per unit interfacial area), mol/(dm’- s) Nb unit rate of bubble gas desorption, mol/(dm’ - s) of quiescent gas dcsorption, N, unit rate mol/(dm2 * s) x effective rate of bubble nucleation, l/(dm3 * s) r rate of gas production in chemical reaction, mol/(dm’ - s) t time, 9 experimental time at which R V has a maximun, s L b parameter in eqn (4) (moment of bubble nucleation), s N

T

temperature,

K

volume of bubbles per unit volume of liquid, dimensionless V volume of liquid, dm3 I gas supersaturation hl liquid (=d#erence between actual bulk and saturation concentrations), mol/dm3

ub

Greek symbols p gas density, kg/dm’ pw water density, hg/dm3

& distribution function of bubble diameters (bubble population density), l/dm4 Subscripts for reagents B bisulphate ton, HSO,-

H hydrogen ion, H+ HN nitrous acid, HNO, N nitrite ion, NW SA sulphamic acid, HaNSOsH W water, Hz0 REFsR8NcEs

Danckwerts P. V., Gas-Liquid Reactions Chap. 11. McGraw-Hill, New York 1970. Aiba S., Humphrey A. E. and Millis N. F., Biochemical Engineer. 2nd Edn, (Polish translation by WNT, Warszawa 1977). University of Tokyo Press, Tokyo 1973. Ebbinghaus L., Ericsson hf. and Lindblom M., Paper 578 presented at 27th IUPAC Congress, Helsiuki 1979. Parker A. L., Siiora L. J. and Hughes R. R, A.LChX.I 1976 22 851. Rudzinski K. J., Ph.D. Thesis, Inst. Physical Chem. Polish Academy Sci., Warszawa 1980 (in Polish). Attar A., A.i.Ch.E.J. 197% 24 106. Szekelv J. and Martins G. P.. Trans. Met. Sot. AIME 1969 245 1721. I91Kinsman G. J. M., Hazeldean G. S. F. and Davies M. W., 1 Iron Steel Inst. I%9 207 1463. [91Bodsworth C. and Bell H. B.. Physical Chemists of Iron and Steel Manufactvre, 2nd Edn. Lcmgman, Loution- 1972. ]lOl Weinberg M. C., Onorato P. I. K. and Uhfman D. R., X Am. Cerum. Sot. 1980 63 175. Ml Astarita G. and Savage D. W., Chem. Engng Sci. 198035 649. U21 Shah Y. T. and Sharma M. M., Trans. Inst. &hem. Engrs 1976 54 1. 1131Joosten G. E. M., Maatman H., Prins W. and Stamhuis E. J., Chem. Engng. Sci. 1980 38 223. (141 Sada E., Kumazawa H., Kudo I. and Koado T., Chem. Engng Sci. 1979 34 719. 1151Sada E.. Kumazawa H. and Hoshino T., Chem. Engng J. 1979 18 125.

and K. I. BuuztQsat

[la] Schweitzer P. hf. and SzebekelyV. G., J. Appl. Phys. 1950 21 1218. [17] Burrows G. and Preece F. H., Tmns. Inst. Chem. Engrs 1954 32 99. [18] Hunt E. B. and Berry V. J., A.1.Ch.E.J. 1956 2 560. 1191 Gent A. N. and Tompkins D. A., J. Appl. Phys. I%9 40 [20] gzer D. E. and Epstein M., Chent. Engng Sci. 1972 27 69. [21] Szekcly J. and Martins 0. P., Trans. Met. Sot. AJME 1969 245 629. [22] Jackson M. L. and Shen Ch.-Ch., A.LCh.EJ. 1978 24 63. [23]Horst A., The Pathological Physiology, p. 404. Polish Medical Publishers, Warszawa 1978 (in Polish). [24)Brian P. L. T., Vivian I. E. and Matiatos D. C., C/rem. Engng Sci. 1%7 22 7. [25] Shah Y. T., Juvekar V. A. and Sharma M. M., Chem. Engng Sci. 1976 31671. [26] Shah Y. T., Pangarkar V. G. and Sharma hf. M., Chem. Engng Sci. 1974 29 1601. [27]CharpentierJ. C., Chemic~IReuclionEngineering RwicwsHouston. ACS Symp. Ser. 72, (Edited by Luss D. and Weekman V. W.), pp. 221-261. ACS, Washington 1978. [28]Tkuy L. T. and Weilaad R. H., Ind. Engng Gem. Fundls 1976 15 286. [29]WeilandR. H., Thuy L. T. and Liveris A. N., Ind. Engng Chem. Fundls 1971 16 332. [30] Bankoff S. G.. Aduan. Chem Engng 1966 6 1. [31]Tao L. N., J. Chem. Phys. 1978 69 4189. 1321 Zak hf. and Weinberg M., J. Non-Cryst. Solids 1980 38-39

533. [33]Pasiuk-BronikowskaW. and RudzinskiK.

I., Chem. Engng sci. 1980 35 5 12. [34] Resnick W. and Gal--Or B., Aduan. Chem. Engng 1968 7 _.‘Y,. [351 Shinnar R., Chemical Reaction Engineerirrg ReuiewsHouston. ACS Symp. Ser. 72. (Edited by Luss D. and Weekman V. W.), pp. l-36. ACS,Washington1978. [36] Aris R. and Amundson N. R., First Order PDEs wirh Applotions. Prentice-Hall, Eoglewood Cliffs, New Jersey 1973. [37] Blander M., Advan. Colloid Inferf Sci. I979 10 1. [38] Blander hf.. Heugstenberg D. and Katz I. L., J. Phys. Chem. 197175 3613. [39] Striven L. E., C/tern. Er~gngSci. 1959 10 I. [40]. Astarita G. and Marrucci C., Ind. Engng Chem. Fundls 1963 2 4. [41] Hughes hf. N., J. Chem. Sot. A 1967 902. [42] Kiug E. J. and King G. W., J. Am. Chem. Sot. 1952 74 1212. [43] Freier R. K., Aqueous Solutions. D&r for Lorganic and Organic Compounds. Walter de Gruyter, Berlin 1976. [44] JXring C. and Gehlen H., Z. anog. al/g. Chem. I%1 312 32. 1451 Hamer W. J., I. Am. C/tern. Sot. 1934 86 869. [46] Doyle G. 1. and Davidson N., I. Am. C&m. Sot. 1949 71 3491. [47] Chan Y. N. I., Birnuaum I. and Lapidus L., Ind. Engng Chem. Fundls 1978 17 133. [48] Vhladsen J. and Michclsen M. L., Solution of Di,,erential Equation Models by Polynomial Approximation, Chap. 8. Prentice-Hall, Englewood Cliffs, New Jersey 1978. APPENDtx

Reaction of sulphamic acid with nitrite HsNSO,- t H+ t NOI- -+Nzt

t SO,‘- + H+ •t HsO.

Hughes[41] thoroughly investigated chemical kinetics of this reaction. Basing on his results (for cus 0.248 mobdm”) we derived the system of formal equations for the reaction at variable pH (amount of HsSO, molecules assumed negligible): dcsAddt =- kX,,CH2C~6sAd~(KnN ’ C&‘&t cnp = c”+~+~+_-6 w

HN

CHQAP Rw H KSA+CH CH

crr)1

(Ala) (Alb)

Gas desorption from liquids %F - cb

= &F-

CSAP

k=(1/60)exp(31.4947-M)81.2/T).

WC)

& - cNF= c”sp - cm

(Al&

c, = cow

(A*e)

1159 (A3)

Reported values of the dissociationconstants are:

c#w= CL, cS.u = &w

t=o

(AID

cm = c$. cEu= CtP According to stoichiometric relations the rate of nitrogen productionis:

r(t) = -dc&dt.

= X122-3792.8/T-O.WlfT;

I&.=7x

[42]

(A4)

lo-‘: 2uT, 1431

KaN = 5.1x 10?

(‘45)

2X, small I, [44]

(A61

log I&,%,,,)= 1.15612 log(7’)- 1387.6/T- 1.355x 1O”T

- 3.8182x IO-‘P + 3.27632; WC,

[4S] (A7)

los(Ka&@) = 2*026V(N(1t d(/tn) -0.1371; 25’C. [46] (AS) (A2)

The temperature dependence of the kinetic constant t can be derived from Hughes’originaldata:

CBSVol.36.No.1-D

log (I&/m)

It is recommendedto solve the system of eqns (Al) with semiimplicit Ruage-Kutta algorithms[47,48],since it may exhibit stiffness.