15 Reactive Systems—Finite Element Analysis CRAIG DOUGLAS and DAVID R O Y L A N C E Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, M A 02139
T h i s paper o u t l i n e s the equations which govern the nonisothermal flow of r e a c t i v e f l u i d s , and d e s c r i b e s the means by which finite element a n a l y s i s can be used t o s o l v e these equations f o r the s o r t of arbitrary boundary c o n d i t i o n s encountered i n industrial p r a c t i c e . The performance of the computer code is i l l u s t r a t e d by s e v e r a l trial problems, s e l e c t e d more f o r t h e i r value in p r o v i d i n g i n s i g h t t o polymer processing flows than as p r a c t i c a l production problems. Although a good d e a l remains t o be learned as t o the performance and proper use of t h i s numerical technique, it is undeniably u s e f u l i n p r o v i d i n g b e t t e r understanding of today's complicated polymer processing problems.
F i n i t e element a n a l y s i s o f f e r s great promise f o r reducing the empiricism o f t e n found i n polymer p r o c e s s i n g design, s i n c e i t i s well suited f o r modeling the complicated boundary c o n d i t i o n s and m a t e r i a l p r o p e r t i e s encountered in industrial practice. Although the method i s now w e l l accepted i n s t r u c t u r a l s t r e s s analysis, its use in fluid transport situations i s less widespread. We have sought t o e x p l o i t some of the many advantages of the f i n i t e element method i n polymer f l u i d p r o c e s s i n g a n a l y s i s , and t h i s paper d e s c r i b e s some of our work i n chemorheology. Reactive flows have s e v e r a l advantages i n polymer p r o c e s s i n g i n comparison with more t r a d i t i o n a l melt-flow techniques. Perhaps the most s i g n i f i c a n t of these i s the energy savings which r e s u l t from the e l i m i n a t i o n of the s e v e r a l melting stages necessary i n such t e c h n o l o g i e s as e x t r u s i o n or i n j e c t i o n molding. However, these advantages are o f f s e t t o some degree by the complexity of the process, which renders the e m p i r i c a l approach t o process development d i f f i c u l t i n the extreme. The flow v e l o c i t i e s i n such processes are governed by the f l u i d v i s c o s i t y , which i s a strong
0097-6156/83/0227-0251$06.00/0 © 1983 American Chemical Society
252
CHEMORHEOLOGY O F THERMOSETTING
POLYMERS
function of temperature and molecular weight. The temperature i n turn i s affected by the viscous dissipation and the heat released or consumed by the reaction, and the reaction rate i s also a strong function of temperature. A l l of these variables interact and change in such a way as to make an intuitive grasp of the process almost impossible, and there i s obviously an advantage to being able to provide some sort of mathematical or numerical simulation of the process. The f i n i t e element scheme to be described below i s very useful for obtaining numerical solutions for reactive flows with arbitrary boundaries, and such a technique i s well suited for detailed analysis of real industrial processes. However, we argue that the greatest value of the method may not be in these detailed calculations, but in the degree to which the method can enhance the process designer's i n t u i t i o n . For this latter purpose, i t i s often sufficient to run only very small and inexpensive t r i a l problems; these can elucidate the manner i n which the various problem parameters interact, so that the designer develops a much improved "feel" for the problem. It i s l i k e l y that most real advances i n processing technology w i l l continue to be made more by Edisonian innovation than by detailed mathematical calculations. But since today's processes have become so complicated, such a technique as the f i n i t e element method can be a powerful adjunct to intelligent i n t u i t i o n . Theoretical background The equations which govern the nonisothermal flow of a reactive f l u i d are derived i n several texts on transport phenomena and polymer processing (e.g. refs. 1,2). Regarding velocity, temperature, and concentration of unreacted species as the fundamental variables, the governing equations can be written as: ptà\x/bt
v
v a
+ u^u] = - p + ç u ν pc[bT/bt + υ · τ ] = Q + k T [èC/dt + u^c] = R + D C V2
V2
(A l i s t of nomenclature i s attached.) The similarity of these equations i s clear. In a l l cases, the time rate of change of the transported variable (velocity, temperature, or concentration) i s balanced by the convective or flow transport terms (e.g. \ι· Ο, the diffusive transport (e.g. D C ) , and a generation term (e.g. R) · The analyst seeks expressions for the space- and timedependent velocities, temperatures and concentrations which sat isfy these equations and also the problem's boundary conditions. It i s generally the boundary conditions which make real problems intractable: even i f one were able to describe the boundaries mathematically, the resulting expressions would not l i k e l y be ame nable to closed-form solution. In addition, many of the ν
V2
15.
Finite Element Analysis
DOUGLAS A N D R O Y L A N C E
253
"constants" i n the above equations a r e o f t e n n o n l i n e a r f u n c t i o n s of the problem v a r i a b l e s . In r e a c t i v e polymer p r o c e s s i n g , one might encounter such expressions as the f o l l o w i n g : η = ç ^ " e x p [ E / R g T ] e x p ( ^ p ) (mw - ) n
1
3
0
4
1
Q = (η/2)(%:%)
+ R(AH)
R = -k exp(E /R T)C m
2
m
0
It i s c l e a r t h a t a l l o f these expressions, taken together, c o n s t i t u t e a mathematical s i t u a t i o n which must be approached with caution. Even though i t i s not o v e r l y d i f f i c u l t t o i n c o r p o r a t e them i n t o a numerical scheme, which we have done, i t i s important t o proceed slowly enough t o develop the proper experience i n the computer code's behavior before t a c k l i n g f u l l - b l o w n problems. Computer Model We have sought t o develop a f i n i t e element code which i s a b l e t o p r e d i c t polymer f l u i d v e l o c i t i e s , s t r e s s e s , temperatures, and degrees o f chemical conversion i n a v a r i e t y o f flow geometries and f o r a v a r i e t y of f l u i d m a t e r i a l p r o p e r t i e s . Space l i m i t a t i o n s prohibit our l i s t i n g here the f u l l d e r i v a t i o n o f t h e f i n i t e e l e ment equations from the above d i f f e r e n t i a l equations, but there e x i s t s e v e r a l w e l l developed means by which t h i s may be done. The reader i s d i r e c t e d t o standard t e x t s ( 3 , 4 ) f o r a more complete d e s c r i p t i o n , and we w i l l just s t a t e here t h a t we employ the G a l e r k i n method o f weighted r e s i d u a l s , together with i s o p a r a m e t r i c mapping and i n t e r p o l a t i o n , t o r e p l a c e the d i f f e r e n t i a l s by i n t e g r a l s which can be evaluated n u m e r i c a l l y over s m a l l subregions ("elements"). The r e s u l t s of these numerical i n t e g r a t i o n s a r e then assembled i n t o a s e t of simultaneous a l g e b r a i c equations which can be solved n u m e r i c a l l y . The salient f e a t u r e s of our code can be l i s t e d b r i e f l y as follows: ( 1 )Velocity, temperature, and chemical conversion a r e taken as nodal unknowns, so t h a t coupled incompressible v i s c o u s flow and d i f f u s i v e - c o n v e c t i v e heat and mass t r a n s p o r t may be mod eled. ( 2 ) Incompressibility i s enforced by a " p e n a l t y " formu lation employing selective reduced i n t e g r a t i o n . T h i s approach requires the use of double p r e c i s i o n a r i t h m e t i c with a concommitant r e d u c t i o n i n the amount of a v a i l a b l e c o r e , but i t has s e v e r a l programming advantages which u s u a l l y outweigh t h i s draw back. ( 3 ) The code i s developed primarily f o r plane and axisymmetric flows. We have coded a three-dimensional c a p a b i l i t y , but g e n e r a l l y f e e l t h a t the expense o f running three-dimensional problems i s not j u s t i f i e d f o r most of our modeling r e s e a r c h . ( 4 ) The code i n c l u d e s s e v e r a l models f o r the e f f e c t o f shear r a t e , temperature, pressure, and chemical conversion on the f l u i d viscosity. These n o n l i n e a r models have not y e t been researched e x t e n s i v e l y , however, and the e x p l o r a t i o n o f the i n t e r a t i v e schemes needed f o r t h e i r proper use c o n s t i t u t e s a major g o a l f o r
254
CHEMORHEOLOGY O F THERMOSETTING
POLYMERS
future research. We have also coded a capability for viscoelastic fluid effects (ref. 5 ) , but currently feel that this d i f f i c u l t aspect of polymer flow i s being researched satisfactorily by other workers. (5) Convective transport of heat or chemical species can be handled either by a conventional Galerkin treatment or by the convenient but s t i l l controversial "optimal upwinding" approach. (6) The code can treat transient problems by means a two-point "theta-method" time-stepping algorithm. The dynamic algorithm is also useful in nonlinear problems, in which the f i n a l fluid state may be approached dynamically from an estimated i n i t i a l state. (7) The code is capable of a variety of iterative treatments of nonlinear problems, including Newton-Raphson iteration and incremental load methods. Some additional discussion is warranted concerning the treatment of convective effects beyond what has been mentioned in item (5) above. Momentum convection (/m» u) i s generally negligible in comparison with the viscous terms due to the high viscosities of polymer f l u i d s , but the convective terms tend to dominate the energy and mass transport equations due to the genera l l y low thermal and mass d i f f u s i v i t i e s . The programming of the convective terms presents no special problems in the Galerkin approach beyond the need to store and solve unsymmetric matrices, but i t is well known that the presence of strong convective terms tends to create oscillations in the f i n a l solutions which can be large enough to destroy their value. This i n s t a b i l i t y i s related to the tendency of convection to produce large downstream gradients which the finite element grid cannot resolve. A largely ad hoc procedure known as "upwinding" has been used in both finite element and finite difference work which seems to alleviate this problem by providing a greater numerical weight to the upstream portion of the element. Hughes (6) has published a very convenient means of upwinding, in which the sampling points in the numerical integration scheme are simply moved upstream an appropriate distance. We have made extensive use of the Hughes upwinding technique, but the reader i s cautioned that this method is regarded as controversial by many workers. A provocative paper by Gresho (7) details many of the possible p i t f a l l s in upwinding, and states a strong preference for grid refinement as the appropriate cure for convection-induced i n s t a b i l i t i e s . 7
Selected Numerical Result The operation of the numerical model w i l l be illustrated by means of three computer experiments, chosen more for their value in visualizing reactive flow than as detailed numerical simulations of real industrial processes. Simple fluid property models w i l l be used in which material parameters such as viscosity and reaction order are not allowed to vary during the process, although real situations are often much more complicated than this. The numerical model does have the capability of performing
15.
DOUGLAS A N D ROYLANCE
Finite Element Analysis
255
a variety of i t e r a t i v e schemes t o model r e a l s i t u a t i o n s i n which n o n l i n e a r and t i m e - v a r y i n g f l u i d p r o p e r t i e s are p r e s e n t , but such complicated simulaions are often not as r e v e a l i n g as simple t r i a l problems, at l e a s t i n the e a r l i e r stages of r e s e a r c h . Nonreactive entry flow. F i g u r e 1 shows the streamlines f o r a 4:1 entry flow. Here a g r i d of 100 four-node l i n e a r elements was used to model the upper symmetric h a l f of a plane c a p i l l a r y , and a f u l l y - d e v e l o p e d P o i s e u i l l e v e l o c i t y was imposed on the r e s e r v o i r entry as a boundary c o n d i t i o n . The streamlines are i d e n t i cal with published experimental and numerical r e s u l t s , although the grid used here was not intended to be f i n e enough to capture the weak r e c i r c u l a t i o n which develops i n the stagnant corner of the r e s e r v o i r . The temperature contours f o r c o n v e c t i o n l e s s flow are shown in figure 2, which shows a hot region at the entrance of the capillary due to the combination of high v i s c o u s energy d i s s i pation there and i t s d i s t a n c e from c o o l boundaries to which heat may be conducted. These isotherms are normalized on the maximum centerline temperature expected for Poiseuille flow in the capillary. The importance of thermal convection r e l a t i v e to conduction is given approximately by the P e c l e t number Pe = UL/?c/k, where U and L are a c h a r a c t e r i s t i c v e l o c i t y and l e n g t h . Figure 3 p l o t s the v a r i a t i o n of temperature along the c e n t e r l i n e f o r v a r i o u s v a l ues of the P e c l e t number, and i t can be seen t h a t the e f f e c t of increased thermal convection i s to sweep the c o o l e r upstream flow particles into the c a p i l l a r y , with a r e s u l t i n g lowering of the temperatures o v e r a l l and a s h i f t downstream of the hot spot near the throat. The r e l a t i v e l y coarse grid used i n t h i s problem produced unstable G a l e r k i n r e s u l t s f o r P e c l e t numbers higher than approximately ten, and the higher degrees of thermal convection were computed using the upwinding f o r m u l a t i o n . Further t e s t s with refined grids should be completed t o assess the accuracy of the upwinded s o l u t i o n s , although the p l o t s i n f i g u r e 3 appear reason able. One-dimensional r e a c t i v e flow. As a p r e l i m i n a r y t r i a l prob lem i n our computations of r e a c t i v e flow, we have studied a simple situation i n which a f l u i d obeying f i r s t - o r d e r chemical k i n e t i c s moves at constant velocity and temperature i n the p o s i t i v e χ direction. Here only the mass-transport equation i s o p e r a t i v e , and i t takes the simple form: 2
2
u(dC/dx) = -KC + D ( d C / d x ) T h i s equation i s solved e a s i l y , and f o r nonzero values of the d i f fusion c o e f f i c i e n t D two boundary values f o r C must be s p e c i f i e d . One of these i s the i n i t i a l c o n c e n t r a t i o n at the i n l e t , and the other requires a consideration of the o u t l e t c o n d i t i o n s . Here s e v e r a l p o s s i b i l i t i e s e x i s t , and we have studied the case i n which
256
CHEMORHEOLOGY O F THERMOSETTING
POLYMERS
Figure 2. Contours of constant temperature f o r convectionless entry flow, with heat generation by viscous d i s s i p a t i o n o n l y .
15.
DOUGLAS A N D ROYLANCE
Finite Element Analysis
257
the c o n c e n t r a t i o n of the o u t l e t r e s e r v o i r i s allowed t o r i s e to meet that s u p p l i e d by the flow; t h i s i s equivalent to s p e c i f y i n g a zero c o n c e n t r a t i o n gradient at the o u t l e t . For the case of n e g l i gible diffusion (D-*0), the second-order term vanishes from the above equation and the downstream boundary c o n d i t i o n cannot be specified. The s o l u t i o n i s then a simple e x p o n e n t i a l , i n which the reactive species vanish according t o f i r s t - o r d e r k i n e t i c s as they are c a r r i e d downstream at constant v e l o c i t y . Figure 4 shows the computed and exact p r e d i c t i o n s f o r r e a c t i v e group c o n c e n t r a t i o n as a f u n c t i o n of d i s t a n c e along the chan nel. The G a l e r k i n values are n e a r l y e x a c t , but i t i s c l e a r that upwinding l e a d s to erroneous r e s u l t s i n the s m a l l - d i f f u s i o n c a s e . The upwinding has introduced an a r t i f i c i a l l y high d i f f u s i v i t y and a zero c o n c e n t r a t i o n g r a d i e n t at the o u t l e t , and such a r t i f a c t s must be considered as p o s s i b i l i t i e s when upwinding i s used. Two-dimensional nonisothermal r e a c t i v e flow. Figure 5 shows the contours of constant conversion f o r a two-dimensional analog of the flow d i s c u s s e d i n the previous s e c t i o n . Here a g a i n , a sim ple uncoupled problem i s t r e a t e d i n which the m a t e r i a l parameters are taken t o be independent of the s o l u t i o n v a r i a b l e s , and i n which the v e l o c i t y c o n d i t i o n s are p r e s c r i b e d . The c o n c e n t r a t i o n i s taken to have a f i x e d value at the i n l e t and a zero g r a d i e n t at the outlet, as b e f o r e . The two-dimensionality of the problem i s contained i n two f e a t u r e s : the v e l o c i t y i s taken to be p a r a b o l i c , ranging from a maximum a t the c e n t e r l i n e to zero at the w a l l s (a Poiseuille flow); and now d i f f u s i v e heat and mass t r a n s p o r t can occur i n both the χ and y d i r e c t i o n . For the low d i f f u s i v i t i e s shown, mass d i f f u s i o n i n the χ d i r e c t i o n i s n e g l i g i b l e , as was demonstrated i n the previous s e c t i o n . However, the c o n c e n t r a t i o n gradients i n the y d i r e c t i o n are s u b s t a n t i a l , so that d i f f u s i v e transport i n that d i r e c t i o n i s a p p r e c i a b l e even at D = 0.01. At D = 0.001, even the y-diffusion is n e g l i g i b l e , so the concen tration contours simply represent a fluid which moves i n the x-direction while r e a c t i n g by f i r s t - o r d e r k i n e t i c s . The concen t r a t i o n s along the c e n t e r l i n e are then i d e n t i c a l with the D = 0.01 curve of f i g u r e 4. Figure 6 shows the contours of constant temperature which result from t h i s flow (with D = . 0 0 1 ) , where the temperature boundary c o n d i t i o n s were set to zero at the entry and along the top and bottom s u r f a c e s . The temperature g r a d i e n t at the o u t l e t was allowed to become z e r o , s i m i l a r to the c o n c e n t r a t i o n g r a d i e n t . The results obtained for the temperature f i e l d are of course dependent on the values chosen f o r f l u i d p r o p e r t i e s . To avoid using space here f o r a d e t a i l e d d i s c u s s i o n of the dimensional a n a l y s i s used f o r s e l e c t i n g these parameters, we w i l l s t a t e simply that in figure 6 the v i s c o u s d i s s i p a t i o n and r e a c t i o n heat make approximately equal contributions to the internal heating (Brinkman and Damkoehler numbers both equal to t h r e e ) .
CHEMORHEOLOGY O F THERMOSETTING
258
1.5
POLYMERS
Γ
-Π
-2
0
h
2
DISTANCE FROM ENTRY Figure 3. numbers.
* 0. 0
Entry flow c e n t e r l i n e
0.2
0.4
temperatures at various
0.6
0. θ
Peclet
1.0
D i eianct, x/L Figure 4. Degree of conversion along channel i n one-dimensional r e a c t i v e flow.
15.
DOUGLAS A N D ROYLANCE
F i g u r e f l o w
a t
5. two
I s o c o n v e r s i o n d i f f e r e n t
259
Finite Element Analysis
c o n t o u r s
mass
i n t w o - d i m e n s i o n a l
d i f f u s i v i t i e s ,
G a l e r k i n
r e a c t i v e c a l c u l a t i o n s .
260
F i g u r e
CHEMORHEOLOGY OF THERMOSETTING
6·
I s o t h e r m a l
G a l e r k i n
c a l c u l a t i o n ,
r e a c t i o n
h e a t i n g .
c o n t o u r s Heat
i n
t w o - d i m e n s i o n a l
g e n e r a t i o n
by
v i s c o u s
r e a c t i v e
POLYMERS
f l o w ,
d i s s i p a t i o n
and
15.
DOUGLAS A N D R O Y L A N C E
Finite Element Analysis
261
Conclusions The numerical model described above has a s i g n i f i c a n t p r e s ent a b i l i t y to simulate a wide range of problems i n polymer p r o c essing. At the same time, it i s small enough to permit easy implementation i n even r a t h e r small processing f a c i l i t i e s , and f o r quick f a m i l i a r i z a t i o n by process d e s i g n e r s . We f e e l such a capa bility would lead to a s i g n i f i c a n t advance i n i n d u s t r y c a p a b i l i t y f o r process development and o p t i m i z a t i o n . Acknowledgments The authors g r a t e f u l l y acknowledge the support of t h i s work by the Army M a t e r i a l s and Mechanics Research C e n t e r , the National Aeronautics and Space A d m i n i s t r a t i o n through MIT's M a t e r i a l s P r o c e s s i n g C e n t e r , and the Hysol-Dexter C o r p o r a t i o n . We a l s o acknowl edge the i n v a l u a b l e a s s i s t a n c e of C o l i n Freese on many aspects of t h i s work, and f o r h i s b e a u t i f u l p r e - and p o s t - p r o c e s s i n g graphics routines. Nomenclature c
c D Ei E
2
k Κ m mw η Ρ
Q R
τ u
fi %
10 Ρ ΔΗ
s p e c i f i c heat concentration of r e a c t i v e species diffusion coefficient a c t i v a t i o n energy f o r v i s c o s i t y a c t i v a t i o n energy f o r r e a c t i o n c o e f f i c i e n t of thermal conduction r e a c t i o n rate preexponential f a c t o r o v e r a l l reaction rate r e a c t i o n order molecular weight power-law exponent f o r v i s c o s i t y pressure i n t e r n a l heat generation i n t e r n a l species generation gas constant temperature v e l o c i t y or v e l o c i t y vector f a c t o r f o r pressure dependency shear rate viscosity viscosity coefficient density heat of r e a c t i o n g r a d i e n t operator
Literature Cited 1. B i r d , R . B . , W . E . Stewart, and E . N . L i g h t f o o t : "Transport Phenomena," John Wiley & Sons, I n c . , New York, 1960.
CHEMORHEOLOGY OF THERMOSETTING
262
POLYMERS
2. Middleman, S.: "Fundamentals of Polymer Processing," Mc Graw-Hill Co., New York, 1977. 3. Zienkiewicz, O.C.: "The Finite Element Method," McGraw-Hill Co., London, 1977. 4. Huebner, K.H.: "The Finite Element Method for Engineers," John Wiley & Sons, Inc., New York, 1975. 5. Collins, B.R., S.M. Thesis, MIT, Cambridge, 1981. 6. Hughes, T.R.J., W.K. Liu, and A. Brooks, J . Comp. Physics, vol. 30, pp. 1-59, 1979. 7. Gresho, P.M., and R.L. Lee, Computers and Fluids, vol. 9, pp. 223-253, 1981. RECEIVED
June 3, 1983
