- Email: [email protected]
Solution of partial differential equations systems by the moving finite element method
Compurers them. Engng,Vol. 16, No. 6, pp. 583-592, 1992 Printedin Great Britain.All rightsreserved
0X%-1354/92 $5.00 + 0.00
Copyright0 1992 F’ergamonPress Ltd
SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS SYSTEMS BY THE MOVING FINITE ELEMENT METHOD C. SERENO’~,
A. RODRIGUES’$and J. VILLADSEN~
’ Department of Chemical Engineering, University of Porto, 4099 Porte Codex, Portugal, r Departmentof Biotechnology, Technical University of Denmark, 2800 Lyngby, Denmark (Received 3 December
199O;finaI revision received 2 December 1991; received for publication 10 December 1991)
Abstrati-Partial differential equations (PDE) systems are solved by the moving finite element method (MFEM) using independent systems of finite elements to each PDE, with polynomial approximations of arbitrary degree in each of the finite elements. The method can be used for any type of linear boundary conditions. A computer code is developed to illustrate the method. This code is modified in the first example in order to solve one PDE coupled with one ordinary differential equation (ODE) in the simulation of a catalytic batch reactor. The code is applied to the following problems: wave progagation, flame propagation and a nonisothermal tubular catalytic reactor described by a pseudo-homogeneous axial dispersion model.
MFEM with polynomial approximations of arbitrary degree in each of the finite-elements using Lagrange polynomials, with the position of the interpolation points optimized (Sereno, 1989) as in the orthogonal collocation method. Every PDE has its own system of finite elements.
INTRODUCHON Modeling of systems and, in particular, of chemical engineering problems often leads to systems of partial differential equations (SPDE) whose solutions contain sharp profiles moving with different velocities. One way to solve these models is to use the moving finite element method (MFEM). After the initial works of Miller and coworkers in the field of MFEM (Miller and Miller, 1981; Miller, 1981; Alexander et al., 1979) related to one PDE in one-dimension, this method was applied in the solution of SPDEs (Gelinas et al., 1981; Djomehri, 1983), using only one system of finite elements and piecewise linear approximation in each of the finite elements. In order to avoid an overdetermined system of ordinary differential equations (SODE), the equations related to the nodes positions have to be weighted in some way. More recently Glasser (1989) developed a general-purpose l-D1 code based on MFEM, still using piecewise linear approximation in every finite element. Baines and Wathen (1985) suggested an MFEM where each of the PDEs has its own system of finite elements; however they still use a piecewise linear approximation in every finite element. We are not always able to solve SPDE by piecewise linear approximation and so we apply here the
APPLICATIONOF THE MFEM TO SPDE Let us consider equation is: !$=f
the system of nPDEs
whose mth
t 7 2 a2Ym ,_( ,x, .a;)~+&+x,K$),
(1)
. . . , y,}‘= vector of the depenwhere j;=(yt,yr, dent variables in which y, = y,,,(x, t) when t > 0 and a c x < b, under the initial condition:
Y, = g,(x ) for t=O
and
and the linear boundary
a
for x=a
thesent address: University Ago&rho Neto, Department of Chemical Engineering8 P.O. Box 1756, Luanda, Angola. #To whom all correspondence should be addressed.
and/or
x-b
and
t >O.
To solve this SPDE we shall use the MFEM polynomial
approximations
of arbitrary
each of the finite elements (Sereno et nI., 1991). 583
with
degree
in
c.
584
fhEN0
There are two ways to define the finite elements: 1. Only one system of moving finite elements for all the differential equations. 2. One independent system of moving finite elements for each differential equation. In the first situation if each of the PDEs leads to fronts moving with different velocities we need more finite elements in order to give a good description of all fronts. In this case the dimension of the implicit SODE will increase as well as the effort to solve it, although the information on some points of the grid related to the equations where their profiles are almost flat is useless. In this case the coefficient matrix of the system and its Jacobian are banded, although the band still has many zero elements. On the other hand, if there are fronts moving in opposite directions, the possibility of nodes crossing will also increase. In the second case we have to use in each system of finite elements the number and the degree of the polynomial approximation just necessary to give an accurate solution of the PDE, and so the dimension of the SODE will be smaller than in case 1, although the Jacobian will be a sparse matrix that will change as the integration proceeds. In spite of these drawbacks we believe that case 2 will perform better and will be followed here. IMPLEMENTATION
OF THE MFEM TO SPDE
Once each PDE has its own system of finite elements, one can apply the MFEM independently to each PDE just as if there was only one PDE. The resulting SODE is the set of ODES obtained by application of the MFEM to each one of the PDEs of the system. We shall only present the general set of ODES obtained by application of the MFEM to the m th PDE. The system of finite elements connected to the mth PDE is obtained by dividing the region a 6 x
t) if XzCx
(2) CXF,,.
(3)
If we define vm.L as: (4)
et d.
and yi is a n,,= - 1 polynomial obtained by Lagrange interpolation through n,,,‘ (after this point we shall always use N instead of nmr and N, instead of n,,, ,) interpolation points whose relative position in the finite element is v;“~ and where y;5; is y$,$, then:
where 8’y‘ is the ith Lagrange interpolation function defined in terms of v”‘vL. On the other hand the derivatives of yi within the finite element L are given by:
where z;l= Xm L+I - X;l is the length element L and
of the finite
The ODES we get from the m th PDE are the following: 1. For Y;,~, where Y;,~ is the approximation of y, at the kth interior node of the finite element L: GzL=O
(10)
2. For Y:.~, where yf,, is the approximation of y,,, at the separation’node Xy+ 1, i.e. where yL,,, is ‘ y,.,=y$‘:
3. For the separation
x (z?G:”
node XT,,:
- z’;, ,GI;‘L+‘)
Solution of PDEs by moving finite element method
13)
14)
585
we use the subroutines JCOBI, DFOPR, RADAU and INTRP developed by Villadsen and Michelsen (1978) to obtain the nodes inside each of the finite elements, the A$ and B$‘ given by equations (8) and (9), respectively, the quadrature points and their weights to evaluate Z$ and Fh, given by equations (14) and (17), respectively, and the values of the Lagrange interpolation functions at these quadrature points. The implicit SODE that results from application of the MFEM is solved by the package LSODI developed at the Lawrence Livermore National Laboratory (Hindmarsh, 1980) with the Jacobian of the SODE calculated numerically.
NUMERICAL
(17)
(18)
Catalytic batch reactor This model was developed by Leitlo (1987) and describes a well-mixed, isothermal, catalytic batch reactor with no film mass transfer limitations. The mass balance to the fluid phase leads to:
dx
.(1+&J em.L=(-&+c+(l+-&y.
dr
*d& =
au
Iu - ,
(21)
VhO- at
-
(19)
(20)
When we need to calculate yf at a given point of the domain a
Considering pore diffusion and catalytic reaction the mass balance for the catalyst particle becomes:
ay
4e$+62 au - ~~~~~ pppy
-= at
&pTd
- Knc(1
CODE
To solve by MFEM any general model defined by an SPDE of the type presented in equation (1)
% - &&,%Y2
(22)
Ep
where x and y are the fluid and solid phase concentrations scaled by the initial concentration c,,, r is the time, 7d is the diffusion time constant, M is the mass of catalyst, a is the rate of change of the fluid volume from the initial volume V,,,,, p_, is the apparent density of the catalyst support, X, is the mass of catalyst per mass of the support and K,,,, and KmM are rate constants. The initial conditions are: x=1 y=O and the boundary
COMPUTING
EXAMPLES
y=x
at
at
r=O,O
condition at
t = 0,
120
is: for
u=l,
once the symmetry condition is expressed by changing the dimensionless radius r to u = cz.
C. SsttsNo et al.
586
0.6
a4
0.2
0.0
I
0.0
0.2
I
I
I
0.4
0.6
,
0.6
I
1.0
E
Fig. 1. Concentration protiles in the solid phase of the batch reactor. Four finite elements with I x (4) + 2 x (0) + 1 x (4); c! = 10k5, cl = 0, c: = lo-‘, c: = 10e2, cl = 10e5 and CA= 10-r; initial condition given by the straight lines joining the pairs (u,~): (0, 0), (0.99985,0), (0.9999, lo-‘), (0.99995.3 x lo-‘), (1.1); r4E[O, 11; l separation nodes.
In this model
the values
of the parameters
are the
following: td = 0.0185 min; pap = 0.27 g cm-r; V,, = 500 cm3; a = 0.4167 cm3 min-r; cp = 0.85;
application of the MFEM to equation (22) we have one more ODE, equation (21). The integral that appears in equation (17) is calculated with seven interior quadrature points in all the finite elements. In Fig. 1 the concentration profiles in the solid phase are represented and in Fig. 2 the evolution of the concentration in the fluid phase. In Fig. 3 the trajectories of the separation nodes are presented.
KmM= 3.57 x 105cm3 g-r min-r; x,
Wave propagation
= 5.4 x 10-3;
K,,,, = 1.55 x 109 cm6 g-r min-’
The propagation of two waves moving in opposite directions presented by Vet-wet et al., (1989) is given by the following model:
mol-‘;
c, = 6.99 x 10e6 mol crn3;
m = 0.01 g.
au at
au 1oouv, ax
A=---
In this problem pore diffusion is the controlling mechanism. The code developed here was modified to take into account that besides the equations resulting from the
i
I
(23)
aV2Ll*” at
.
ax
0.6
Y
--
.
0
50
t
I&
Ii0
Fig. 2. Evolution of the concentration in the fluid phase of the batch reactor.
J
0.0
0.2
Fig. 3. Trajectories
0.4
5
0.6
of separation reactor.
0.6
nodes
1.0
for the batch
Solution
Boundary
of PDEs by moving finite element method
conditions:
u=O
for
u= 0
for
x=-O.5 x = 0.5
and and
t20,
0(x, 0) =
In Fig. 4 the amplitudes of the waves and in Fig. 5 the trajectories of the separation nodes are represented.
t 2 0.
Initial conditions:
u(x, 0) =
587
Flame propagation
0
-0.5
d x < -0.3,
0.5[1 +cos(loxX)]
-0.3
d X d -0.1,
0
-0.1
6 x GO.5,
0
-0.5
6 X < -0.1,
0.5[1 + c0s(10zx)]
0.1 B x < 0.3,
0
0.3 G x < 0.5.
The integral that appears in equation (17) is calculated with six interior quadrature points in all the finite elements.
This problem was also proposed by Venver et al. (1989) and the model obtained after the application of the mass and energy balances is the following: :’
I au 3% - - - 3.52 x lO%e-““, z-ax* (24)
I
2
$ = $+
3.52 x 106ue-““,
subjected to the initial conditions tr = 1, u = 0.2 throughout the space variable domain and to the
I.09
(cl 0.6 0.6u,v 0.4. 0.2 0.0-1 -0.24 -0.5
0.6
-0.3
-0.1
x
0.1
0.3
0.5
-0.3
-0.1
x
0.1
0.3
0.5
0.6
u,v
u,v
4.:o,,,,,5
x
a.21 -0.5
Fig. 4. Amplitudes of the waves at I: 0, 0.1, 0.25, 0.3 and 0.4. For both equations--9
finite elements
withlx(3)+2x(0)+3x(3)+2x(O)+lx(3);c;=10-5,c~=O,c~=10-3,c;=l0-1,c;=l0-~and CT = 1O-5: initial position of the separation nodes for II: -0.50; -0.30005, -0.3, -0.29995, -0.23333, -0.16667, -0.10005, -0.1, -0.09995,0.5; initial position oftheseparationnodeaforv: -0.50,0.09995, O,l, 0.10005, 0.16667. 0.23333, 0.29995, 0.3, 0.30005, 0.5; XE[-0.5,OJk + separation nodes for u; 0 separation nodes for 0.
588
C. Smmo
et al.
1.2 1.0 J
0.8 " 0.6 0.4 0.2
“._
-0.2-I -OS
6.3
-0.1
x
0.1
0.3
0.5
Fig. 5. Trajectories of separation nodes for wave propagation.
following
boundary
conditions:
au
-=?=O ax
for
ax
au
ax =O
x=0
I 2 0,
t
x = 1
O< t $0.0002,
and
o=1.2forx=1
This model does not allow us to use directly the actual code to solve SPDE in the interval 0 < t & 0.0002 because the boundary condition makes u change with time for x = 1. hence in this time interval we have defined a new variable y given by the ODE:
for
at
condition x=1
0.0
02
0.4
Pe X2
ay
(26)
t = 0.
for v is: O
and
and the code was modified to take into account this new ODE just as in the first example. For t > 0.0002 we use the general code.
-0.d
x
110
ax -=%?a,_D~uexp[-y(~-l)] au T -=_ ---++Daxv de zhl ay
under the initial condition:
v=y
018
This model was taken from Ferreira (1988). It simulates a tubular catalytic reactor through the pseudo-homogeneous model with axial dispersion in the fluid concentration. The mass and energy balances arc:
80
The new boundary
016
NON-ISOTHERMAL TUBULAR REACTOR
andt~0.0002.
y = 0.2
x
The integral that appears in equation (17) is calculated with 33 interior quadrature points in all the finite elements. In Fig. 6 the concentration profiles u, in Fig. 7 the temperature profiles a, in Fig. 8 the trajectories of the separation nodes of II and in Fig. 9 the trajectories of the separation nodes of v are represented.
U = 0.2 + 0.0002 for
014
ta0,
and
x = 1 and
for
0:2
0.0
Fig. 7. Temperature profiles for flame propagation. Four finite elements with 1 x (5) + 2 x (10) + 1 x (5); CT = 0, CT = 0, c; = lo-“, CT = IO-s, cy = 10-s and c;: = 10e5; initial condition given by u = 0.2 including in z = 1; initial position of the separation nodes: 0, 0.994, 0.996, 0.998, 1; x E[O, 11; + separation nodes.
0.6
0.8
1.0
Fig. 6. Concentration profiles for flame propagation. Initial condition given by u = 1 including in z = 1; + separation nodes (all the other elements are the same as in 0).
where x, v centration, temperature temperature.
0.0
and u, are respectively the fluid the fluid temperature and the scaled by the feed concentration < and 0 are the space variable and
0.2
0.4
x
0.6
0.8
conwall and time
1.0
Fig. 8. Trajectories of separation nodes for concentration in flame propagation.
Solution of PDEs by moving finite element method
1.2
589
1 -
U.”
0:2
0.0
0:4
x
0:s
Oi
Fig. 9. Trajectories of separation nodes for temperature in flame propagation.
scaled by the length of the bed L and space time 7. Pe is the Peclet number, Da is the Damkiihler number, y is the Arrhenius number, fl is the adiabatic temperature rise, 2, is the time constant for the thermal wave propagation and NW,,is the number of transfer units for heat transfer at the wall. The initial conditions are: x=0
at
f?=O
and
O<[
v=l
at
0=0
and
O
and the boundary
conditions
dn -=Pe(y-1) X
for
are: [=O
and
020,
0.4
0.2
0.0
1:o
for
[=I
and
for
In this problem values:
C=O
and
the parameters
T/Q, = 2.08 x 10-4; Da=0.7; NWh= 33.7;
920. have the following
y = 21.8;
/l =0.7; Pe=
1.0
ci = 0, c: = IO-), c: = 10e2, cl = IO-’ and c: = 10m5;initial condition given by v = 1 including in 5 = 0; initial position of the separation nodes; 0, 5 x IO-‘, IO-‘l, 1; C E[O, I]; + separation nodes.
The concentration wave is traveling 104 faster than the temperature wave and hence the concentration first reaches values equal to those of an isothermal reactor at the feed temperature. Only much later the influence of the changing temperature is felt. In Fig. 10 the concentration profiles, in Fig. 11 the temperature profiles, in Fig. 12 the concentration histories, in Fig. 13 the temperature histories, in Fig. 14 the trajectories of the separation nodes of x and in Fig. 15 the trajectories of the separation nodes of v are represented.
AND REMARKS
030,
a[ v=l
0.8
Fig. Il. Temperature profiles for the tubular reactor. Three finite elements with 1 x (10) f2 x (5); ct = 10-5,
CONCLUSIONS
ax -=0
0.6 5
104. .
The conclusions we have obtained by solving SPDEs resulting as models of dynamic systems by the MFEM with polynomial approximation of arbitrary degree in each finite element, are just similar to those already obtained for only one PDE (Sereno et al., 1991). The use of independent systems of finite elements for each PDE has proved to be a powerful way to decrease the possibility of nodes crossing which allows us to use much smaller values of the
0.8
1 0.0
0.2
0.4
r
4 0.6
0.8
1.0
0.0
0.2
0.4
5
0.6
0.8
Fig. 10. Concentration protiles for the tubular reactor. Twelve finite elements with 1 x (5) + 10 x (0) + 1 x (5); initial condition given by x = 0 including in [ = 0; initial position of the separation nodes: 0, 5 x 10-5, IO-“, 1.5 x 10-4, 2 x 10-4, 2.5 x 10-4 3 x 10-4, 3.5 x 10-4, 4 x lo-‘, 4.5 x IO--‘, 5 x 10-4, 5.5 x 10m4, 1; + separation nodes (all the other elements are the same as in a).
1.0
C. SnrzrJo et al.
0.2
_J
0.0 i
;,
m
Fig. 12. Concentration
-em histories
1200
for the tubular
1500
reactor.
1.016
1.012
v
1.006
1.004
1.000
0
600
300
Fig. 13. Temperature
histories
for the tubular
1500
1200
900
reactor.
0.6
Fig. 14. Trajectories
of separation
nodes for concentration
in the tubular
reactor.
Solution
of PDEs by moving
finite clement method
591
S,, = “Spring function”
in finite element L of the mth system of finite elements I = Time variable T, = Feed temperature U = Global heat transfer coefficient at the reactor wall u, = Interstitial velocity at the reactor entrance umJ = Space variable in finite element L of the mth system of finite elements u = Temperature scaled by T, o, = Wall temperature scaled by T, I+~ = Space variable in finite element L of the mth system of finite elements normaliid to [0, l] V, =Initial volume of the fluid phase in batch reactor x = Space variable; fluid phase concentration scaled by c0 X, = Mass of catalyst per mass of the support L+, = Separation node between finite Xm elements L and L + 1 of the mth system of finite elements y = Solid phase concentration scaled by co in batch reactor yf;;= Approximation of the mth dependent variable in finite element L of its own system of finite elements ym = m th dependent variable of the SPDE 2; = Length of fmite element L of the mth system of finite elements
1.016-
5
Fig. 15. Trajectories
r
nodes for temperaturein the tubular reactor. of separation
coefficients of the penalty functions. However the Jacobian of the right-hand side of the system resulting from the application of the MFEM to equation (1) is a sparse matrix that is changing as the integration proceeds. In order to improve the efficiency of the method a special effort has to be done in improving the algebraic part of the integrator. NOTATION a,
b = End points of the spaa variable domain A$ = Parameter defined by equation (8) B$‘ = Parameter defined by equation (9) c0 = Initial concentration in the batch reactor; feed concentration in the tubular reactor cy = Parameters in the definition of crnL
Greek symbols (-AH)c,
B= ~
and LL
cmr= Specific heat of the fluid
k(TO) = __ r = Damkiihler number
E y = = Arrhenius R T,
E
E = -AH = k, K,, K,,,M = L = m = N E n,., =
N, =?u+I
2iJT *uh = G,, Pe = g
Activation energy Heat of reaction Rate constants Reactor length Mass of catalyst Total number of interpolation points in finite element L of the rnth system of finite elements = Total number of interpolation points in finite element L + 1 of the mth system of finite elements
Number or transfer units to p the heat transfer at the wall
= Peclet number = Q t q,,, 0 Number of finite elements in the system associated to the mth PDE R = Particle radius; ideal gas constant &, = Reactor radius
temperature
rise
6 = Ekd porosity Ep= Particle porosity %I,~= “Viscosity function” in finite element L of the mth system of finite elements
CL= S&cific heat of the solid 0, = Axial dispersion coefficient II)&= Effective diffusivity Da
= Adiabatic
TOP&
pap = pr = p, = 5=
Apparent density Fluid density Solid density Space time
R2 q, = = Diffusion &I %I =
8 PrCps+ (1 -
8 P&r
elP,cp.
number
time constant
Time constant for the thermal ’ = wave propagation 0 = Time scaled by 5 C = Particle radius variable scaled by R C = Distance from the tubular reactor entrance scaled by L
RJEFERENCFS
Alexander R., P. Mansclli and K. Miller, Moving Finite Elemmta for the Stefa Problem in Two Dimensions, Rend. Accad. Naz. Lincci (Rome), Serie VIII, Vol. LXVII, fast. 1-2, pp. 57-61 (1979).
592
c.
%?RENO
Baines M. J. and A. J. Wathen, Moving finite element modelline of comoressible flow. Reoort 4/85. Universitv of Readiig, U.K.-(1985). _ _ Djomehri M. D., Moving finite element solution of systems of partial differential equations in one-dimension. Ph.D. Thesis, University of California, Berkeley (1983). Ferreira R. M., Contribui9So para o estudo de reactores cataliticos em leito fixo: efeito da convec@o em catalisadores de poros largos e cases de catalisadores bidispersos. Ph.D. Thesis, University of Porto, Portugal (1988). Gelinas R. J., S. K. Doss and K. Miller, The moving finite element method: applications to general partial differential equations with multiple large gradients. J. Comput. f’hys. 40, 202-249 (1981). Glasser A. H., A moving finite element model of the high density Z-pinch. J. Comput. Phys. 85, 159-209 (1989). Hindmarsh A. C., LSODE and LSODI, two new initial value ordinary differential equation solvers. ACMSIGNUM Newslett. 15, IO-11 (1980).
et al.
L&Ho A., Estudo do processo MEROX. Ph.D. Thesis, University of Porto, Portugal (1987). Miller K., Moving finite elements, Part II. SIAM J. Namer. Anal. 18, 1033-1057 (1981). Miller K. and R. N. Miller, Moving finite elements, Part I. SIAM J. Numer. Anal. 18, 1019-1032 (1981). Serene C., MCtodo dos elementos finitos mbveis: applica9&s em engenharia qulmica. Ph.D. Thesis, University of Porto, Portugal (1989). Sereno C., A. Rodrigues and J. Villadsen, The moving finite element method with polynomial approximation of any degree. Computers them. Engng 15, 25-33 (1991). Verwer J. G., J. G. Blom and J. M. Sanz-Serna, An adaptive moving grid method for one dimensional systems of partial differential equations. J. Comput. Pbys. 82, 454486. Villadsen J. and M. L. Michelsen, Solution of Dtjkrentiaf Equations Models by Polynomial Approximation. PrenticeHall, Englewood Cliffs, NJ (1978).
0X%-1354/92 $5.00 + 0.00
Copyright0 1992 F’ergamonPress Ltd
SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS SYSTEMS BY THE MOVING FINITE ELEMENT METHOD C. SERENO’~,
A. RODRIGUES’$and J. VILLADSEN~
’ Department of Chemical Engineering, University of Porto, 4099 Porte Codex, Portugal, r Departmentof Biotechnology, Technical University of Denmark, 2800 Lyngby, Denmark (Received 3 December
199O;finaI revision received 2 December 1991; received for publication 10 December 1991)
Abstrati-Partial differential equations (PDE) systems are solved by the moving finite element method (MFEM) using independent systems of finite elements to each PDE, with polynomial approximations of arbitrary degree in each of the finite elements. The method can be used for any type of linear boundary conditions. A computer code is developed to illustrate the method. This code is modified in the first example in order to solve one PDE coupled with one ordinary differential equation (ODE) in the simulation of a catalytic batch reactor. The code is applied to the following problems: wave progagation, flame propagation and a nonisothermal tubular catalytic reactor described by a pseudo-homogeneous axial dispersion model.
MFEM with polynomial approximations of arbitrary degree in each of the finite-elements using Lagrange polynomials, with the position of the interpolation points optimized (Sereno, 1989) as in the orthogonal collocation method. Every PDE has its own system of finite elements.
INTRODUCHON Modeling of systems and, in particular, of chemical engineering problems often leads to systems of partial differential equations (SPDE) whose solutions contain sharp profiles moving with different velocities. One way to solve these models is to use the moving finite element method (MFEM). After the initial works of Miller and coworkers in the field of MFEM (Miller and Miller, 1981; Miller, 1981; Alexander et al., 1979) related to one PDE in one-dimension, this method was applied in the solution of SPDEs (Gelinas et al., 1981; Djomehri, 1983), using only one system of finite elements and piecewise linear approximation in each of the finite elements. In order to avoid an overdetermined system of ordinary differential equations (SODE), the equations related to the nodes positions have to be weighted in some way. More recently Glasser (1989) developed a general-purpose l-D1 code based on MFEM, still using piecewise linear approximation in every finite element. Baines and Wathen (1985) suggested an MFEM where each of the PDEs has its own system of finite elements; however they still use a piecewise linear approximation in every finite element. We are not always able to solve SPDE by piecewise linear approximation and so we apply here the
APPLICATIONOF THE MFEM TO SPDE Let us consider equation is: !$=f
the system of nPDEs
whose mth
t 7 2 a2Ym ,_( ,x, .a;)~+&+x,K$),
(1)
. . . , y,}‘= vector of the depenwhere j;=(yt,yr, dent variables in which y, = y,,,(x, t) when t > 0 and a c x < b, under the initial condition:
Y, = g,(x ) for t=O
and
and the linear boundary
a
for x=a
thesent address: University Ago&rho Neto, Department of Chemical Engineering8 P.O. Box 1756, Luanda, Angola. #To whom all correspondence should be addressed.
and/or
x-b
and
t >O.
To solve this SPDE we shall use the MFEM polynomial
approximations
of arbitrary
each of the finite elements (Sereno et nI., 1991). 583
with
degree
in
c.
584
fhEN0
There are two ways to define the finite elements: 1. Only one system of moving finite elements for all the differential equations. 2. One independent system of moving finite elements for each differential equation. In the first situation if each of the PDEs leads to fronts moving with different velocities we need more finite elements in order to give a good description of all fronts. In this case the dimension of the implicit SODE will increase as well as the effort to solve it, although the information on some points of the grid related to the equations where their profiles are almost flat is useless. In this case the coefficient matrix of the system and its Jacobian are banded, although the band still has many zero elements. On the other hand, if there are fronts moving in opposite directions, the possibility of nodes crossing will also increase. In the second case we have to use in each system of finite elements the number and the degree of the polynomial approximation just necessary to give an accurate solution of the PDE, and so the dimension of the SODE will be smaller than in case 1, although the Jacobian will be a sparse matrix that will change as the integration proceeds. In spite of these drawbacks we believe that case 2 will perform better and will be followed here. IMPLEMENTATION
OF THE MFEM TO SPDE
Once each PDE has its own system of finite elements, one can apply the MFEM independently to each PDE just as if there was only one PDE. The resulting SODE is the set of ODES obtained by application of the MFEM to each one of the PDEs of the system. We shall only present the general set of ODES obtained by application of the MFEM to the m th PDE. The system of finite elements connected to the mth PDE is obtained by dividing the region a 6 x
t) if XzCx
(2) CXF,,.
(3)
If we define vm.L as: (4)
et d.
and yi is a n,,= - 1 polynomial obtained by Lagrange interpolation through n,,,‘ (after this point we shall always use N instead of nmr and N, instead of n,,, ,) interpolation points whose relative position in the finite element is v;“~ and where y;5; is y$,$, then:
where 8’y‘ is the ith Lagrange interpolation function defined in terms of v”‘vL. On the other hand the derivatives of yi within the finite element L are given by:
where z;l= Xm L+I - X;l is the length element L and
of the finite
The ODES we get from the m th PDE are the following: 1. For Y;,~, where Y;,~ is the approximation of y, at the kth interior node of the finite element L: GzL=O
(10)
2. For Y:.~, where yf,, is the approximation of y,,, at the separation’node Xy+ 1, i.e. where yL,,, is ‘ y,.,=y$‘:
3. For the separation
x (z?G:”
node XT,,:
- z’;, ,GI;‘L+‘)
Solution of PDEs by moving finite element method
13)
14)
585
we use the subroutines JCOBI, DFOPR, RADAU and INTRP developed by Villadsen and Michelsen (1978) to obtain the nodes inside each of the finite elements, the A$ and B$‘ given by equations (8) and (9), respectively, the quadrature points and their weights to evaluate Z$ and Fh, given by equations (14) and (17), respectively, and the values of the Lagrange interpolation functions at these quadrature points. The implicit SODE that results from application of the MFEM is solved by the package LSODI developed at the Lawrence Livermore National Laboratory (Hindmarsh, 1980) with the Jacobian of the SODE calculated numerically.
NUMERICAL
(17)
(18)
Catalytic batch reactor This model was developed by Leitlo (1987) and describes a well-mixed, isothermal, catalytic batch reactor with no film mass transfer limitations. The mass balance to the fluid phase leads to:
dx
.(1+&J em.L=(-&+c+(l+-&y.
dr
*d& =
au
Iu - ,
(21)
VhO- at
-
(19)
(20)
When we need to calculate yf at a given point of the domain a
Considering pore diffusion and catalytic reaction the mass balance for the catalyst particle becomes:
ay
4e$+62 au - ~~~~~ pppy
-= at
&pTd
- Knc(1
CODE
To solve by MFEM any general model defined by an SPDE of the type presented in equation (1)
% - &&,%Y2
(22)
Ep
where x and y are the fluid and solid phase concentrations scaled by the initial concentration c,,, r is the time, 7d is the diffusion time constant, M is the mass of catalyst, a is the rate of change of the fluid volume from the initial volume V,,,,, p_, is the apparent density of the catalyst support, X, is the mass of catalyst per mass of the support and K,,,, and KmM are rate constants. The initial conditions are: x=1 y=O and the boundary
COMPUTING
EXAMPLES
y=x
at
at
r=O,O
condition at
t = 0,
120
is: for
u=l,
once the symmetry condition is expressed by changing the dimensionless radius r to u = cz.
C. SsttsNo et al.
586
0.6
a4
0.2
0.0
I
0.0
0.2
I
I
I
0.4
0.6
,
0.6
I
1.0
E
Fig. 1. Concentration protiles in the solid phase of the batch reactor. Four finite elements with I x (4) + 2 x (0) + 1 x (4); c! = 10k5, cl = 0, c: = lo-‘, c: = 10e2, cl = 10e5 and CA= 10-r; initial condition given by the straight lines joining the pairs (u,~): (0, 0), (0.99985,0), (0.9999, lo-‘), (0.99995.3 x lo-‘), (1.1); r4E[O, 11; l separation nodes.
In this model
the values
of the parameters
are the
following: td = 0.0185 min; pap = 0.27 g cm-r; V,, = 500 cm3; a = 0.4167 cm3 min-r; cp = 0.85;
application of the MFEM to equation (22) we have one more ODE, equation (21). The integral that appears in equation (17) is calculated with seven interior quadrature points in all the finite elements. In Fig. 1 the concentration profiles in the solid phase are represented and in Fig. 2 the evolution of the concentration in the fluid phase. In Fig. 3 the trajectories of the separation nodes are presented.
KmM= 3.57 x 105cm3 g-r min-r; x,
Wave propagation
= 5.4 x 10-3;
K,,,, = 1.55 x 109 cm6 g-r min-’
The propagation of two waves moving in opposite directions presented by Vet-wet et al., (1989) is given by the following model:
mol-‘;
c, = 6.99 x 10e6 mol crn3;
m = 0.01 g.
au at
au 1oouv, ax
A=---
In this problem pore diffusion is the controlling mechanism. The code developed here was modified to take into account that besides the equations resulting from the
i
I
(23)
aV2Ll*” at
.
ax
0.6
Y
--
.
0
50
t
I&
Ii0
Fig. 2. Evolution of the concentration in the fluid phase of the batch reactor.
J
0.0
0.2
Fig. 3. Trajectories
0.4
5
0.6
of separation reactor.
0.6
nodes
1.0
for the batch
Solution
Boundary
of PDEs by moving finite element method
conditions:
u=O
for
u= 0
for
x=-O.5 x = 0.5
and and
t20,
0(x, 0) =
In Fig. 4 the amplitudes of the waves and in Fig. 5 the trajectories of the separation nodes are represented.
t 2 0.
Initial conditions:
u(x, 0) =
587
Flame propagation
0
-0.5
d x < -0.3,
0.5[1 +cos(loxX)]
-0.3
d X d -0.1,
0
-0.1
6 x GO.5,
0
-0.5
6 X < -0.1,
0.5[1 + c0s(10zx)]
0.1 B x < 0.3,
0
0.3 G x < 0.5.
The integral that appears in equation (17) is calculated with six interior quadrature points in all the finite elements.
This problem was also proposed by Venver et al. (1989) and the model obtained after the application of the mass and energy balances is the following: :’
I au 3% - - - 3.52 x lO%e-““, z-ax* (24)
I
2
$ = $+
3.52 x 106ue-““,
subjected to the initial conditions tr = 1, u = 0.2 throughout the space variable domain and to the
I.09
(cl 0.6 0.6u,v 0.4. 0.2 0.0-1 -0.24 -0.5
0.6
-0.3
-0.1
x
0.1
0.3
0.5
-0.3
-0.1
x
0.1
0.3
0.5
0.6
u,v
u,v
4.:o,,,,,5
x
a.21 -0.5
Fig. 4. Amplitudes of the waves at I: 0, 0.1, 0.25, 0.3 and 0.4. For both equations--9
finite elements
withlx(3)+2x(0)+3x(3)+2x(O)+lx(3);c;=10-5,c~=O,c~=10-3,c;=l0-1,c;=l0-~and CT = 1O-5: initial position of the separation nodes for II: -0.50; -0.30005, -0.3, -0.29995, -0.23333, -0.16667, -0.10005, -0.1, -0.09995,0.5; initial position oftheseparationnodeaforv: -0.50,0.09995, O,l, 0.10005, 0.16667. 0.23333, 0.29995, 0.3, 0.30005, 0.5; XE[-0.5,OJk + separation nodes for u; 0 separation nodes for 0.
588
C. Smmo
et al.
1.2 1.0 J
0.8 " 0.6 0.4 0.2
“._
-0.2-I -OS
6.3
-0.1
x
0.1
0.3
0.5
Fig. 5. Trajectories of separation nodes for wave propagation.
following
boundary
conditions:
au
-=?=O ax
for
ax
au
ax =O
x=0
I 2 0,
t
x = 1
O< t $0.0002,
and
o=1.2forx=1
This model does not allow us to use directly the actual code to solve SPDE in the interval 0 < t & 0.0002 because the boundary condition makes u change with time for x = 1. hence in this time interval we have defined a new variable y given by the ODE:
for
at
condition x=1
0.0
02
0.4
Pe X2
ay
(26)
t = 0.
for v is: O
and
and the code was modified to take into account this new ODE just as in the first example. For t > 0.0002 we use the general code.
-0.d
x
110
ax -=%?a,_D~uexp[-y(~-l)] au T -=_ ---++Daxv de zhl ay
under the initial condition:
v=y
018
This model was taken from Ferreira (1988). It simulates a tubular catalytic reactor through the pseudo-homogeneous model with axial dispersion in the fluid concentration. The mass and energy balances arc:
80
The new boundary
016
NON-ISOTHERMAL TUBULAR REACTOR
andt~0.0002.
y = 0.2
x
The integral that appears in equation (17) is calculated with 33 interior quadrature points in all the finite elements. In Fig. 6 the concentration profiles u, in Fig. 7 the temperature profiles a, in Fig. 8 the trajectories of the separation nodes of II and in Fig. 9 the trajectories of the separation nodes of v are represented.
U = 0.2 + 0.0002 for
014
ta0,
and
x = 1 and
for
0:2
0.0
Fig. 7. Temperature profiles for flame propagation. Four finite elements with 1 x (5) + 2 x (10) + 1 x (5); CT = 0, CT = 0, c; = lo-“, CT = IO-s, cy = 10-s and c;: = 10e5; initial condition given by u = 0.2 including in z = 1; initial position of the separation nodes: 0, 0.994, 0.996, 0.998, 1; x E[O, 11; + separation nodes.
0.6
0.8
1.0
Fig. 6. Concentration profiles for flame propagation. Initial condition given by u = 1 including in z = 1; + separation nodes (all the other elements are the same as in 0).
where x, v centration, temperature temperature.
0.0
and u, are respectively the fluid the fluid temperature and the scaled by the feed concentration < and 0 are the space variable and
0.2
0.4
x
0.6
0.8
conwall and time
1.0
Fig. 8. Trajectories of separation nodes for concentration in flame propagation.
Solution of PDEs by moving finite element method
1.2
589
1 -
U.”
0:2
0.0
0:4
x
0:s
Oi
Fig. 9. Trajectories of separation nodes for temperature in flame propagation.
scaled by the length of the bed L and space time 7. Pe is the Peclet number, Da is the Damkiihler number, y is the Arrhenius number, fl is the adiabatic temperature rise, 2, is the time constant for the thermal wave propagation and NW,,is the number of transfer units for heat transfer at the wall. The initial conditions are: x=0
at
f?=O
and
O<[
v=l
at
0=0
and
O
and the boundary
conditions
dn -=Pe(y-1) X
for
are: [=O
and
020,
0.4
0.2
0.0
1:o
for
[=I
and
for
In this problem values:
C=O
and
the parameters
T/Q, = 2.08 x 10-4; Da=0.7; NWh= 33.7;
920. have the following
y = 21.8;
/l =0.7; Pe=
1.0
ci = 0, c: = IO-), c: = 10e2, cl = IO-’ and c: = 10m5;initial condition given by v = 1 including in 5 = 0; initial position of the separation nodes; 0, 5 x IO-‘, IO-‘l, 1; C E[O, I]; + separation nodes.
The concentration wave is traveling 104 faster than the temperature wave and hence the concentration first reaches values equal to those of an isothermal reactor at the feed temperature. Only much later the influence of the changing temperature is felt. In Fig. 10 the concentration profiles, in Fig. 11 the temperature profiles, in Fig. 12 the concentration histories, in Fig. 13 the temperature histories, in Fig. 14 the trajectories of the separation nodes of x and in Fig. 15 the trajectories of the separation nodes of v are represented.
AND REMARKS
030,
a[ v=l
0.8
Fig. Il. Temperature profiles for the tubular reactor. Three finite elements with 1 x (10) f2 x (5); ct = 10-5,
CONCLUSIONS
ax -=0
0.6 5
104. .
The conclusions we have obtained by solving SPDEs resulting as models of dynamic systems by the MFEM with polynomial approximation of arbitrary degree in each finite element, are just similar to those already obtained for only one PDE (Sereno et al., 1991). The use of independent systems of finite elements for each PDE has proved to be a powerful way to decrease the possibility of nodes crossing which allows us to use much smaller values of the
0.8
1 0.0
0.2
0.4
r
4 0.6
0.8
1.0
0.0
0.2
0.4
5
0.6
0.8
Fig. 10. Concentration protiles for the tubular reactor. Twelve finite elements with 1 x (5) + 10 x (0) + 1 x (5); initial condition given by x = 0 including in [ = 0; initial position of the separation nodes: 0, 5 x 10-5, IO-“, 1.5 x 10-4, 2 x 10-4, 2.5 x 10-4 3 x 10-4, 3.5 x 10-4, 4 x lo-‘, 4.5 x IO--‘, 5 x 10-4, 5.5 x 10m4, 1; + separation nodes (all the other elements are the same as in a).
1.0
C. SnrzrJo et al.
0.2
_J
0.0 i
;,
m
Fig. 12. Concentration
-em histories
1200
for the tubular
1500
reactor.
1.016
1.012
v
1.006
1.004
1.000
0
600
300
Fig. 13. Temperature
histories
for the tubular
1500
1200
900
reactor.
0.6
Fig. 14. Trajectories
of separation
nodes for concentration
in the tubular
reactor.
Solution
of PDEs by moving
finite clement method
591
S,, = “Spring function”
in finite element L of the mth system of finite elements I = Time variable T, = Feed temperature U = Global heat transfer coefficient at the reactor wall u, = Interstitial velocity at the reactor entrance umJ = Space variable in finite element L of the mth system of finite elements u = Temperature scaled by T, o, = Wall temperature scaled by T, I+~ = Space variable in finite element L of the mth system of finite elements normaliid to [0, l] V, =Initial volume of the fluid phase in batch reactor x = Space variable; fluid phase concentration scaled by c0 X, = Mass of catalyst per mass of the support L+, = Separation node between finite Xm elements L and L + 1 of the mth system of finite elements y = Solid phase concentration scaled by co in batch reactor yf;;= Approximation of the mth dependent variable in finite element L of its own system of finite elements ym = m th dependent variable of the SPDE 2; = Length of fmite element L of the mth system of finite elements
1.016-
5
Fig. 15. Trajectories
r
nodes for temperaturein the tubular reactor. of separation
coefficients of the penalty functions. However the Jacobian of the right-hand side of the system resulting from the application of the MFEM to equation (1) is a sparse matrix that is changing as the integration proceeds. In order to improve the efficiency of the method a special effort has to be done in improving the algebraic part of the integrator. NOTATION a,
b = End points of the spaa variable domain A$ = Parameter defined by equation (8) B$‘ = Parameter defined by equation (9) c0 = Initial concentration in the batch reactor; feed concentration in the tubular reactor cy = Parameters in the definition of crnL
Greek symbols (-AH)c,
B= ~
and LL
cmr= Specific heat of the fluid
k(TO) = __ r = Damkiihler number
E y = = Arrhenius R T,
E
E = -AH = k, K,, K,,,M = L = m = N E n,., =
N, =?u+I
2iJT *uh = G,, Pe = g
Activation energy Heat of reaction Rate constants Reactor length Mass of catalyst Total number of interpolation points in finite element L of the rnth system of finite elements = Total number of interpolation points in finite element L + 1 of the mth system of finite elements
Number or transfer units to p the heat transfer at the wall
= Peclet number = Q t q,,, 0 Number of finite elements in the system associated to the mth PDE R = Particle radius; ideal gas constant &, = Reactor radius
temperature
rise
6 = Ekd porosity Ep= Particle porosity %I,~= “Viscosity function” in finite element L of the mth system of finite elements
CL= S&cific heat of the solid 0, = Axial dispersion coefficient II)&= Effective diffusivity Da
= Adiabatic
TOP&
pap = pr = p, = 5=
Apparent density Fluid density Solid density Space time
R2 q, = = Diffusion &I %I =
8 PrCps+ (1 -
8 P&r
elP,cp.
number
time constant
Time constant for the thermal ’ = wave propagation 0 = Time scaled by 5 C = Particle radius variable scaled by R C = Distance from the tubular reactor entrance scaled by L
RJEFERENCFS
Alexander R., P. Mansclli and K. Miller, Moving Finite Elemmta for the Stefa Problem in Two Dimensions, Rend. Accad. Naz. Lincci (Rome), Serie VIII, Vol. LXVII, fast. 1-2, pp. 57-61 (1979).
592
c.
%?RENO
Baines M. J. and A. J. Wathen, Moving finite element modelline of comoressible flow. Reoort 4/85. Universitv of Readiig, U.K.-(1985). _ _ Djomehri M. D., Moving finite element solution of systems of partial differential equations in one-dimension. Ph.D. Thesis, University of California, Berkeley (1983). Ferreira R. M., Contribui9So para o estudo de reactores cataliticos em leito fixo: efeito da convec@o em catalisadores de poros largos e cases de catalisadores bidispersos. Ph.D. Thesis, University of Porto, Portugal (1988). Gelinas R. J., S. K. Doss and K. Miller, The moving finite element method: applications to general partial differential equations with multiple large gradients. J. Comput. f’hys. 40, 202-249 (1981). Glasser A. H., A moving finite element model of the high density Z-pinch. J. Comput. Phys. 85, 159-209 (1989). Hindmarsh A. C., LSODE and LSODI, two new initial value ordinary differential equation solvers. ACMSIGNUM Newslett. 15, IO-11 (1980).
et al.
L&Ho A., Estudo do processo MEROX. Ph.D. Thesis, University of Porto, Portugal (1987). Miller K., Moving finite elements, Part II. SIAM J. Namer. Anal. 18, 1033-1057 (1981). Miller K. and R. N. Miller, Moving finite elements, Part I. SIAM J. Numer. Anal. 18, 1019-1032 (1981). Serene C., MCtodo dos elementos finitos mbveis: applica9&s em engenharia qulmica. Ph.D. Thesis, University of Porto, Portugal (1989). Sereno C., A. Rodrigues and J. Villadsen, The moving finite element method with polynomial approximation of any degree. Computers them. Engng 15, 25-33 (1991). Verwer J. G., J. G. Blom and J. M. Sanz-Serna, An adaptive moving grid method for one dimensional systems of partial differential equations. J. Comput. Pbys. 82, 454486. Villadsen J. and M. L. Michelsen, Solution of Dtjkrentiaf Equations Models by Polynomial Approximation. PrenticeHall, Englewood Cliffs, NJ (1978).