- Email: [email protected]
Studies of finite element procedures—the use of ADIANA-F in fluid flow analyses
Compulm & Slructure~ Vol. 32. No. 3/4. pp. 499-516, 1989 Printed in Greet Britain.
0
1989 Mawell
004.5.7949/89 s3.M) + 0.00 Pqamon Maanillan plc
STUDIES OF FINITE ELEMENT PROCEDURESTHE USE OF ADINA-F IN FLUID FLOW ANALYSES? KLAUS-J~RGENBATHE$ and JIANDON~ $Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A. §ADINA R & D, Inc., Watertown, MA 02172, U.S.A.
Abstract-We present in this paper some experiences obtained with the program ADINA-F in fluid flow analyses. The theory and numerical methods used in the program for laminar and turbulent flow, and flow with heat transfer, are summarized and then the solutions of various problems are presented. These problem solutions comprise the analyses of turbulent flow in a pipe and natural convection in a porous medium, the prediction of flow around a cylinder, the conjugate heat transfer in flow through a pipe and in the cooling of electronic equipment, and the 3D analysis of natural convection in a cavity.
I.INTRODUCTION
flow with heat transfer are, using index notation and the usual summation convention:
The solution of practical fluid flow problems by finite element methods has obtained increasing attention during recent years [l, 21. In particular, the solution of incompressible fluid flow phenomena with or without heat transfer can now be performed efficiently in many practical cases. In an earlier paper, we presented the basic formulation and finite element procedures that we are using for fluid flow analysis [3]. These techniques have been implemented in the computer program ADINA-F, which can be employed for the solution of 2D and 3D fluid flow problems with or without heat transfer. We have continued with the development of this computer program, and the objective of this paper is to present some recent experiences obtained with ADINA-F in the modeling of fluid flow problems. In the paper, we first briefly summarize the basic continuum mechanics equations solved and finite element procedures used. This presentation emphasizes those aspects of the theory that are important when modeling fluid flow problems and interpreting the calculated results. However, the main emphasis of the paper is on the solution of some complex fluid flow problems with ADINA-F. These solutions relate to modeling laminar and turbulent flow, and flow with heat transfer, and demonstrate how the program can be applied in engineering practice.
momentum:
constitutive:
(2) continuity: (3)
(g+vie,i )
(4)
heat transfer:
Pep
=(ke,i),j+q8
where velocity of fluid flow in direction xi mass density, constant components of stress tensor components of body force vector pressure Kronecker delta fluid viscosity, temperature-dependent components of velocity strain tensor = I(vij + Oj,i> specific heat at constant pressure, temperaturedependent temperature thermal conductivity, temperature-dependent rate of heat generated per unit volume, temperature-dependent (frictional effects are included as in Example 3.5)
2. FLUID FLOW EQUATIONS The objective in this section is to briefly summarize the fluid flow equations considered for solution, with an emphasis on the physical and practical aspects that are important when establishing finite element models. 2.1. Continuum mechanics equations Using a Cartesian reference frame (xi, i = 1, 2, 3) the fluid flow equations considered for incompressible
Oi,i =o
and in natural convection problems
t Original figures were generated with a colour producing
fB=dl
terminal and submitted in colour. 499
-iw -43
KLAUS-JURGEN BATHE and JUN DONG
500
where gi
B
0,
acceleration of gravity into direction i coefficient describing the change in mass density per unit change in temperature reference temperature.
The boundary (l)-(4) are:
conditions
corresponding
to eqns
vective terms in eqns (1) and (4), the convection and radiation boundary conditions in eqn (8) and the temperature-dependent material properties. In addition to the nonlinearities listed above, we also include those that arise due to the use of turbulence models (see Sec. 2.6) and the use of non-Newtonian constitutive relations. For solution of eqn (1) it is effective to substitute eqns (2) and (3) into eqn (1) to obtain
prescribed fluid velocities, 4, on the surface S,, ui =
prescribed tractions, ff,
cilS,
p(d) + (3
on the surface S,,
7ijnj=_ffls2
(6)
where nj are the components of the unit normal to the fluid surface and f f are the components of the traction vector; prescribed temperature,
=
-P,i
+ff
+
/iui,#ii.
(12)
Note that the incompressibility condition is now contained in eqn (12). To establish effective finite element models it can also be expedient to consider the nondimensionalized governing equations. If we use the characteristic length L, velocity V, reference temperature &,, reference temperature increment Af3, and
x:
v: = Vi/V,
=x,/L,
P * =p/pv,
8, on the surface S,, 6 =&,
Uj”i,j)
f”
e*,e-eO
(7)
--s-
t* = tV/L
=f tgL/pV qBL2
q*=k
(13)
*e
prescribed heat flux into the surface S,,
(8)
we obtain as the governing momentum equations, ti: +u,?vit,=
where qs is the heat flux input to the body. We note that for eqns (S)-(8) S = S,, fl= S,nS2 and S = m, 8 = S,nS,, where S is the total surface area of the fluid body. The heat flow input in eqn (8) comprises, firstly, the effect of actually applied distributed heat flow; secondly, the effect of convection heat transfer, in which case qS=h(e,-eS)
(9)
where h is the convection coefficient, and 6, and fIs are the environmental and body surface temperatures; and thirdly, the effect of radiation heat transfer, in which case qs =
K(e, - es)
K = h,[(w2 + (eS)%e, + 6’)
(10) (11)
and h, is the radiation coefficient and 6, is the temperature of the radiating body [4]. The value of h, is determined from the Stefan-Boltzmann constant, the emissivity of the radiating and absorbing materials and the geometric view factors. The above equations are the standard Navier-Stokes equations governing the motion of a viscous incompressible fluid in laminar flow with heat transfer. Inherent nonlinearities are due to the con-
-p:
and energy
+fY +(lIRe)v&
(14)
and
Pe(B*+O:e,*,)=e:i+q*
(15)
where the Reynolds number is
Re.!%; v=!.! V
P
and the Peclet number is
pe=yL; c,=k. a
PCP
The use of the nondimensionalized equations in the finite element solution is achieved by simply using the following material properties: mass density = unit value 1 . v1scosny = G heat capacity = Pe conductivity = unit value.
(16)
Of course, in the above equations we implicitly assumed that the fluid properties are constant. If these properties vary with temperature, eqns (14) and
ADINA-F in fluid fiow analysis
(15) are still applicable but the properties to be used in these equations are then tem~ratu~~e~dent. Note that we may also interpret eqns (14) and (15) as the governing fluid flow equations with the material properties in eqn (16) in any consistent set of units. When using eqns (13), (14) and (15) for solution of a fluid flow problem, it is of course also necessary to impose the velocity, traction and heat transfer boundary conditions in nondimensionalized form. These boundary conditions are a ~statement of eqns (5+0-(l) using the nondimensionaiized variables given above. For the momentum equation we use v’t to impose boundary velocities, and to impose tractions we recognize that
zij*=
-p*a,+j$v~j+v$)*(17)
For the energy equation we use the values 8* to prescribe temperatures and qs* to impose heat flux conditions,
B
s+ _
qSL
-z-z
The appropriate boundary convection and radiation conditions are also easily derived.
501
heat transfer:
PCJ@ + v,B,,)68 dk’ +
k6,, 60,, d ?’ sV
+
S@K(~, - es) dS
(21)
s&
where S,, S, and S, are part of S, and correspond, respectively, to the surfaces subjected to heat flow input, convection and radiation conditions. Consider next a single finite element for solution of eqns (19x21). The finite element equations for an assemblage of elements are obtained by the usual direct addition of element matrices [S]. For the finite element discretization we use for 22) (planar and axisymmetric) solutions the 9/4 node element shown in Fig. 1: the nine nodes are used to interpolate the geometry, the velocities and the temperature
2.2. Finite element equations We presented the finite element formulation used in detail in [3]. Here we summarize briefly the major ingredients of the formulation and resulting equations and point out some important features regarding the practical application of these equations. The finite element solution of the continuum mechanics equations given above is obtained by establishing a weak form of the equations using the Galerkin procedure [5]. The momentum equations are weighted with the velocities, the continuity equation is weighted with pressure and the energy equation is weighted with the temperature. Using integration by parts gives the variational equations to be discretized by finite element inte~lations: momentum:
=S;f:d~jdV+~sJ’6v;dS
e= i
h&e&;
(24)
k-1
the four nodes are used to interpolate
P= i
k-l
the pressure
$Pk.
Correspondingly, we use in 3D analysis the 27/8 node element (see Fig. 1) for a triquadratic interpolation of the geometry, velocities and temperature and a trilinear interpolation of the pressure. These interpolations provide stable elements as expressed by the inf-sup condition of Brezzi and Babugka [6]. The finite element governing equations are now established as
(19)
M,.* + (K,,,+ Km. Iv + K, p = R, + Rs where 6 means ‘variation in’ (or a virtual quantity); continuity:
K&v=0 Ce+K,e+K,~,~v=Q,+Q,+Q,+Q,
s
Spv,,,dV=O
V
(25)
(20)
(26) (27)
(281
where the coefficient matrices and load vectors are dire&y obtained by substituting the finite element
502
KLAU-JORGEN Bnrns and
1
W&
Fig. 1. 2D and 3D elements available in ADINA-F.
interpolations [see eqns (22)-(25)] into the variational eqns (19H21) with eqns (2) and (3). The details of the matrices in eqns (26x28) are given in [3]. 2.3. Boundary conditions in finite element solutions
The boundary conditions to a problem governed by eqns (l)-(4) are those listed in eqns (5x8). Here we note that eqns (5) and (7) are the essential boundary conditions and eqns (6) and (8) are the natural boundary conditions. The essential boundary conditions are imposed in the finite element analysis by directly setting the solution variables in eqns (26) and (28) equal to the given values. Note that the velocities may be prescribed in a coordinate system different from the Cartesian system, which then requires a standard transformation [S]. The natural boundary conditions correspond to prescribed tractions ff and prescribed heat fluxes qs (which may be given in terms of convection or radiation conditions). As for prescribed velocities, the imposed tractions can be defined in any coordinate system, which again only requires a standard transformation Considering the finite element solution, the natural boundary conditions of eqn (8) are imposed by the right-hand-side vectors in eqn (28). However, as for imposing the natural boundary conditions of eqn (6), consider the vector R, given for an element side as RS=
HSTfs dS
s .%
(29)
where HS is the standard element surface displacement interpolation matrix and B is the vector of the imposed force effects [3]. In the solution F is given
JIAN DONG
and substituted into eqn (29). Note that in general this does not mean that the pressure on the element side is prescribed. Further, we recall that since eqn (26) contains the continuity equation, the components in ‘fs are in general not given by eqn (6) [3]. The physical tractions would be used in F if the continuity equation were not imposed onto the momentum equation. In analyses there is little difference in the results using either approach unless an actual traction is prescribed (and for this case we have in ADINA-F a special option of solution). However, we have also found that using eqn (26) the convergence in the iterative solution is sometimes better than when not imposing the continuity condition onto the momentum equation. Since we are solving incompressible flow problems, particular attention needs to be given to the boundary condition on the pressure. Namely, the pressure is undetermined within a constant unless prescribed at, at least, one point. In the finite element solution we have that: all pressure degrees of freedom in eqn (26) can be left free if the component in F normal to the boundary (which involves the pressure) is prescribed (to zero or any other value) at least at one nodal point on the boundary; if this component is not prescribed anywhere along the boundary, then it is necessary to arbitrarily assign a pressure value at one node in the finite element mesh. The need to prescribe a pressure value at one node in the finite element mesh arises therefore when closed fluid regions are analyzed for which the velocities are prescribed at all points of the boundary (see Sec. 3). 2.4. Time stepping and iterative equation solution With our solution procedures we consider steadystate analyses and transient analyses. In steady-state solutions, the terms M,ir and Cd in eqns (26) and (28) are neglected, and the time stepping merely corresponds to an iterative incremental analysis of the nonlinear response (51. The number of solution steps are selected to arrive at an efficient analysis and to obtain the solution response at specific levels of prescribed loadings and boundary conditions. In transient analyses, we use the a-method with f < c( < 1 with the assumptions ‘i.b’+~+A$_$,)
(30)
~++=(l
(31)
_._~~()‘v+~‘+*‘v
f+zAf(j='(,+A,
_'e)
At
r+xAq = (1 _ a)‘@+ E’+A’g
(32) (33)
ADINA-F in fluid flow analysis and ‘+*A’p= (1 - a)‘p + a’+Alp.
(34)
Usually, we use a = $ {trapezoidal rule) or a = 1 (Euler backward method). A transient analysis can be performed with variable time step sizes (as input by the user) and in steadystate and transient solutions an automatic-time-stepping (ATS) method can be used. In the ATS method, the program decreases automatically the load or time step sizes from those input by the user if convergent in the iterations is not reached for the current load or time considered or if certain accuracy criteria are not satisfied. Finally, we should mention that the incremental solution using eqns (26)--(28) is performed by the method of successive substitution or the full Newton-Raphson iteration (with or without line searches). In cases with strong nonlinea~ties, the Newton-Raphson iteration is more effective. 2.5. Calculated quantities-the stream faction values
use of vorticity and
The solution of the governing finite element equations [eqns (26)-(28)] yields directly the velocities, pressure and temperature over the fluid domain. All these quantities are continuous over the element boundaries. As additional quantities it can be useful to calculate the vorticity, w, of the fluid o=vxv.
(35)
In 2D analysis the stream function values, $, can also be useful. These are calculated by the line integral
$ -$a=
I
(v,h-vsW
(36)
where & is a constant. We note that the vorticity of the fluid is obtained by differentiation of the calculated velocity field, hence for coarse finite element models discontinuities across element boundaries will be observed. Therefore, we can use the isobar& of vorticity to identify whether a finite element mesh was fine enough to predict an unknown flow field accurately [7j. The stream function values are obtained by integrating the calculated velocity field. Here we note that the stream function at a point is only unique if the continuity condition, eqn (3), is satisfied. Hence, if the calculated solution does not satisfy the continuity ~ndition accurately, different values for $ will be obtained depending on which path of integration is chosen. This observation shows that the calculated values of the stream function can also be employed to assess whether a finite element mesh was fine enough to predict a flow field.
503
In the ADINA-F implementation we integrate from a starting point (a finite element node) along the element sides and around the element boundaries, always only along those element boundaries which attach to nodes for which rj has not yet been calculated. Hence, the integration is only performed along an element boundary once and only if it leads to the calculation of a yet unknown value of JI at a node on that boundary. Using this algorithm we can conclude that, on display of the calculated stream function values, nonsmooth lines indicate regions where the mesh was not fine enough to predict the fluid fIow accurately. 2.6. Turbulence modeling One of the most difficult aspects in many fluid flow problems is the modeling of the turbulent behavior of the fluid. The field of modeling turbulence has attracted much attention and various models have been proposed and are in current use [S, 91. We have implemented in ADINA-F the use of the mixing length model and the large-eddy-simulation (LES) model. For the mixing length model, the users can define by their own prog~mming the specific equations for the mixing length 1 to be used. The general form of the model is given by It,=P~21JI p-1
(37)
* -i dul+duI ax, axi
( )
where p, is the turbulence (eddy) viscosity to be added to the laminar viscosity, A, to obtain the total viscosity p. For the LES model we use
where k, is a dimensionless constant and A is taken to be the mean of the side lengths of the element for which P, is being computed. An application of the mixing length model is presented in Sec. 3.1.
3. EXAMPLE SOLUTIONS
The following solution results have been obtained with ADINA-F, where we should note that no specific attempt was made to use optimal meshes or optimal solution approaches with the program. The solution results are presented using the postprocessor ADINA-PLOT. 3.1. Analysis of turbulent flow in a pipe Figure 2(a) shows the pipe considered. The flow is to be solved for a high Reynolds number, so that turbulence effects are important. For the analysis a mesh of 20 x 1 axisymmetric elements was used with
KLAUS-JORGEN BATHE and
504
JUN
DONG
one layer to model a typical section of the fluid in the pipe. To include the turbulence effects, the mixing length model proposed by Nikuradse [8] was employed 1 = 0.14 - 0.08~~ - 0.06~~
(40)
where 2 = i/R, q = r/R, 1 is the Prandtl mixing length, R is the radius of the pipe and r is measured from the centerline of the pipe. This model is adequately used for the flow not very close to the wall. To include the near-wall effects we used for the element at the wall the Nikuradse formula with the van Driest correction 1 = (0.14 - 0.08r~’ - 0.06~~) x
1 -exp
-“U,
rl)R
(
0 NIKURADSE A AOINA-F. X ADINA-F.
EXPERIMENT. RE=3.24E4 RE=2.3E4
RE=2.3E4
(4,) )
where r7= (TO/~)“~, v = h/p, k is the laminar viscosity and 70 is the shear stress at the wall. We programmed the above relations into the usersupplied subroutine of ADINA-F to define the mixing length I. In this case, the application of the mixing length model gives for the total viscosity
Figure 2 shows the solution obtained for Reynolds numbers of 2.3 x lo4 and 3.24 x 104. These numbers are based on the average fluid velocity in the pipe.
-8.0
Heat genwatio~ in
ONIKURADSE A ADINA-F. X AOINA-F.
DISTANCE
EXPERIMENT RE-3.24E4 RE=2.3E4
porous medium
8
FROM WALL
Fig. 2. (a) Analysis of turbulent flow in a pipe. Twenty equal sized elements into radial direction are used to model the fluid. (b) Velocity distribution in the pipe. (c) Virtual viscosity distribution in the pipe.
w-0,
Insulated
c----l.0
Fig. 3(a)
qB
8-O
505
ADINA-F in fluid flow analysis
1EtlPERATURE TlnE 1.000 - 0.1400
lEtlPERATURE TlHE 1.000 - 0.1200
(b)
0.1200 0.0900
-
0.0900
-
0.0400
-
0.0100
- 0.0600
- 0.0300
- 0.0100
RamlO*
TEMPERATURE TIME 1.000 - 0.02500
TEMPERATURE TINE 1.000 - O.llOO
0.0900
-
0.0600
-
0.0300
-
0.0100
0.02000
- 0.01500
-
*Ra.lO'
RarlO'
Ra.105
(d
+ ADINA-F. RA=IDO - ANALYTICAL VERTICAL
CENTEhNE
2. COORDINATE
VERllCAL
CENlERLlNE
COORDlNAlE
-
Fig. 3. (a) Natural convection in a porous medium. (b) Streamlines for Merent values of the Darcy-Rayleigh number. (c) Temperature contours for different values of the Darcy-Rayleigh number. (d) Velocity along the vertkal centerline. (e) Temperature along the vertical centerline.
0.01000
- 0.00500
506
KLAUS-JORGEN
BATHE and
Figure 2(b) shows the velocity distribution and Fig. 2(c) gives the nondimensional&d virtual viscosity defined as
v*=--- P
pGR ’
The figures show good correspondence with the experimental results of Nikuradse. Note, however, that the velocity prediction near the wall is still quite approximate due to the very high velocity gradient. Also, we should note that the solution for the virtual viscosity is given at the element integration points. 3.2. Natural convection in a porous enclosure with internal heat generation
Figure 3(a) shows the problem considered. For this class of problems the governing continuum mechanics equations are vi., = 0
P
KVi =
-P,i
(PJp)fi +
+
JIAN LANG
3.3. Wake behind a cylinder at low Reynolh numbers Figure 4(a) shows the problem considered and schematically the twin vortices and trail in the wake of the cylinder. Our objective was to predict these vortices and the pressure distribution on the cylinder for a range of low Reynolds numbers. Figure 4(b) presents the two finite element meshes that we have used. The analyses with both meshes gave almost the same streamlines and vorticity results; however, mesh A was used to obtain a better pressure prediction. Figure 4(c) shows the streamlines calculated for Re = 40. We clearly identify the twin vortices, which are also shown in Fig. 4(d) by the plot of the vorticity. Figure 4(e) gives a comparison of the predicted wake lengths for different Reynolds numbers with the experimental results of Taneda [ll], and Fig. 4(f) shows the pressure distributions on the cylinder, as observed by Thorn [12] and in this finite element analysis. The nondimensionalized pressure is defined as
Pgi(l - BCe- eO))
$PcpVje,i
where C$is the porosity,
(l;81i),i + 4’
=
K
*_P-PO P -1 TPVO
where p. is the inflow pressure and V. is the inflow velocity.
is the permeability,
BEp= 4PCp + (1- d)(PCJS E=4k+(l--i$)k,
(46) (47)
and fi is the viscosity of the fluid accounting for porosity. The subscript s is used to identify the solid matrix properties, and properties without a subscript are those of the fluid. The above equations are directly obtained from eqns (l)-(4) by including Darcy’s law with the Boussinesq approximation, and hence are an application of the general equations solved in ADINA-F. This problem is characterized by a Darcy-Rayleigh number
3.4. Conjugate heat transfer in a pipe withJluidJow Figures 5(a) and (b) show the problem considered and the finite element discretization used for the analysis. The objective is to study the temperature variation in the fluid and the wall. As a first case, we solved the problem when the heat flux is imposed along the complete pipe section considered and the fluid flow is a fully developed laminar flow. The solution not near the entrance and the exit sections should be in accordance with the equation [ 131,
where Ob= bulk temperature where L is the width of the enclosure and 6 is the effective thermal diffusivity. A 20 x 40 element mesh (with equal sized elements) was used to model the domain of the enclosure. Figures 3(b-e) give the calculated steady-state solutions. These correspond to the results given by Haajizadeh et al. [lo]. We note that the solution for Ra = lo5 was obtained by adding the term p,,(vy + u,,,),, with ~,/p = 0.001 to the right-hand-side of eqn (44) in order to obtain faster convergence in the iterations. Figures 3(d) and (e) shows the vertical velocity and temperature distributions along the vertical centerline for the case of Ra = 100.
v = 2v,,,
j Ovr dr of fluid = Ivr dr
r*
( > 1- -
R*
c(= thermal diffusivity of fluid qi = interface heat flux V aYe = average velocity of fluid k, = fluid conductivity and z is the axial coordinate of the pipe. Figure 5(c) shows wall, pipe axis and bulk temperatures of the fluid (here we have tl = 0.2, qi = 1.l, u.,~ = 0.5, k,= l.O), and shows that the correct variations of temperature along the pipe are predicted.
507
ADINA-F in fluid flow analysis Next, the case of Fig. S(a) was considered, in which the heat flux is imposed on the pipe wall for a finite length. This type of problem is usually characterized by the Peclet number, Pe, and the ratio of the solid and fluid conductivities, K, where Pe = Re*Pr K = k,/k,
24, K Re (Reynolds number) = V
Pr (Prandtl number) = f
Reynolds number is used for each solution. The Eckert numbers based on the maximum temperature differences in the fluid are calculated to be 1.5 and 0.5, respectively. As expected, the figures show that, as the Peclet number increases, the heat transfer by convection becomes more predominant. The finer mesh (mesh B) gives a better streamline solution, and therefore a better temperature prediction. Note that in this analysis the solution for the velocities and temperatures were obtained simultaneously, where actually the solution could be calculated by first solving for the velocities and then for the temperatures. 3.6. Natural convection in a 30 cavity
v = fluid kinematic viscosity k, = solid conductivity k, = fluid conductivity.
The results of the bulk temperature distributions and interface heat fluxes along the pipe for different Peclet numbers and values of K are shown in Figs. S(d-g). 3.5. Conjugate heat transfer in cooling of electronic equipment The problem considered here was mainly chosen to investigate some of the interesting features available in our fluid flow analysis capabilities. We consider a channel of fluid flow with two equal obstructing bodies (the electronic equipment) in its flow that generate heat. The heat is transported in the fluid by conduction and convection. We have no comparison with experimental or other analytical results available, but want to show the predicted temperature distributions for different flow conditions. The problem and finite element meshes used are shown in Figs. 6(a) and (b). The meshes A and B were used in the solution for comparison purposes. We note that the flow is prescribed at the left entrance of the channel. The solid bodies generating heat are modeled with fluid elements that only have the temperature degree of freedom (velocities and pressure are set to zero). The conductivity of the solid is set to be the same as that of the the fluid. The governing equations of the fluid flow considered are those discussed in Sec. 2 with the heat generation due to friction included @ = 2pe,e,.
(49)
The effect of this term is discussed extensively in [8] and is characterized by the Eckert number, EC. The complete problem is characterized by the Peclet, Eckert and Reynolds numbers. Figures 6&-f) give the calculated temperature distributions for different Peclet numbers (10 and 100) and the streamlines of the flow corresponding to a Reynolds number 100, for meshes A and B. The same C.A.S. 32,34-B
The solution of a natural convection problem considering two-dimensional motion was presented by many authors (see also [3]). To study a threedimensional motion, we considered the problem shown in Fig. 7(a). No slip conditions were assumed at all surfaces. The temperature was prescribed at the side surfaces and adiabatic conditions were assumed at the other surfaces. The governing equations are those discussed in Sec. 2 but including the Boussinesq term, pg(l - /?(0 - I$,)), in the z-momentum equation. We required the solution for Ra = lo4 and Pr = 0.2. For the finite element analysis, we considered only one half of the domain due to symmetry conditions with respect to the plane x = 0. The mesh used is shown in Fig. 7(b). Figures 7(c) and (d) give velocity distributions at the planes y = 0.5 and z = 0.5, respectively, for different values of x = constant. As expected the y-velocity field shown in Fig. 7(c) is similar in shape to the velocity field in a two-dimensional analysis (31, but the magnitude of the velocities, of course, decreases as x increases. Similarly, for the z-velocity field given in Fig. 7(d). An interesting result is obtained for the x-velocity field, which shows recirculation near the edges. Figure 7(e) presents the calculated temperature field. It is seen that for this problem the temperature distribution is still very nearly two-dimensional. 4. CONCLUDING REMARKS
The objective in this paper was to summarize some of the latest developments we have undertaken for finite element analysis of fluid flow. We first summarized the theory and numerical methods used, and then presented some example solutions. It can be concluded that ADINA-F is a powerful tool for finite element solution of a variety of fluid flow problems. However, considering analyses of incompressible fluid flows in general, there are of course still a number of important topics for research and development; these include the development of a robust algorithm to predict free surfaces, the development of new more efficient iterative solution schemes,
508
KLAUS-J~~RQEN BATHE and JIAN DONO
t. w
I
T v.v, _, -
z-
(0)
(b)
STREAH FUHCT IOH TIME 1.000 1.000
Fig. 4(a-c)
-
-0.500
-
-1 .ooo
OMEGA-YZ
TIRE i.OOC
(d)
OEXPERIHENTAL BY TANEDA A ADINA-F. tlESH A X ADINA-F. HEW B
REYNOLDS
1956
NUMBER
*. i-
0.
DEGREES
I
40.
FRDII
FD:iARG
*
S~~NA,*O~O~N~
Fig. 4. (a) Flow past a cylinder, Re = pV,d/~, schematic representation. (b) Top, mesh A, bottom, mesh B. (c) Streamlines at Re = 40, mesh A. (d) Isobands of vorticity at Re = 40, mesh A. (e} Lengths of vortioes for different Reynolds numbers. (f) Pressure distribution around cylinder. 509
KLAUS-J~RGENBATHEand JMN DONG
510
t= 0.1
(a)
RI 1.0
FIVE FLUID RADIALLY
ONE SOL IO RADIALLY
+ WALL 0 BULK X PIPE
“07
Id.
5.
AXIAL
Ii.
ELEMENTS
ELEMENT
AXIS
0 PE-5 A PE=S + PE-5
K=I K=5 K-10
2d.
DISTANCE
AXIAL
Fig. 5(a-d)
DISTANCE
511
ADINA-F in fluid flow analysis
(e)
0 PE=50 A PE-50 + PE=SO
AXIAL
0PE=5 A PE-5 + PE-5
K=lO K-50 K=lOO
DISTANCE
AXIAL
0 PE-50 A PE-50 + PE=50
AXIAL
K=l K=5 K=IO
DISTANCE
K=IO K-50 K=lOO
DISTANCE
Fig. 5. (a) Conjugate heat transfer in an “infinitely long” pipe. with fluid flow. A section of length 20 is modeled. (b) Mesh used. (c) Temperatures along the pipe for applied heat flux along complete pipe. (d) Interface heat flux, Pe = 5, K = 1,s and 10. (e) Interface heat flux, Pe = 50, K = 10,SO and 100. (f) Fluid bulk temperature, Pe = 5, K = 1.5 and 10. (g) Fluid bulk temperature, Pe = 50, K = lo,50 and 100. Note that for the case Pe = 50, a longer section pipe should be modeled.
512
KLAUS-JORGENBATHEand JUN DONG
:bl
TEMPERATURE TIME 1.000 ~5.000 -
4.000
--
3.000
--
2.000
--
I .ooo
Pe = 10 TEMPERATURE TIME 1.000 J_
I .700
T
I .500
0.900
t Pe = 100 Fig. 6(a-c)
0.600
513
ADINA-F in fluid flow analysis TEMPERATURE TIME 1.000 l.
(d)
5.000
4.000
I 3.000
2.000
I .ooo
Pe=
10
TEMPERATURE TIME 1.000 1
I .700
-
1 *so0
--
1.200
--
0.900
- - 0.600 Pe = 100 STREAM FUNCTlON TIME 1.000 - -0.1000
(e)
s
-
-0.3000
-
-0.5000
-
-0.7000
-0.8000
Fig. 6(d,e)
KLAUS-J~~RGEN BATHEand JIAN DONG
514
STREAM FUNCTION TIME 1.000 - -0.200
(f)
-
-0.500
-
-0.800
Fig. 6. (a) Analysis of conjugate heat transfer in cooling of electronic equipment. (b) Top, mesh bottom, mesh B. (c) Predicted temperatures using mesh A. (d) Predicted temperatures using mesh B. Streamlines using mesh A, Re = 100. (f) Streamlines using mesh B, Re = 100. top surface, insulated (a)
back surface, insulated
--.
cdd side surface,&
front surface, ins&ted
L, - L,d.O L, -0.4
hot side surface, oh
L
,’
/
/ /’
/
/
insulated
(b)
Fig. 7(a,b)
515
ADINA-F in fluid flow analysis
0 AT x-o.02 A
AT
+A1
x-o.10
XAT 0
+
AT
X-O.10
Q
X=0.06 x=0.14
AT
X=0.18
d 7
v; 7 0.0
0.2
:ObRDINATES . 0.6
0.8
I .o
0.8
I.0
2
d) v;
0 AT x-o.02
d 4’
AT
A
AT
X=0.06
+
AT
x=0.10
X
AT
x=0.14
0
AT
X=0.18
X.0.0
X
AT
X-O.14
0
AT
X=0.18
v; * Cm 8’0 _I-& y* it
0.0
0.2
0.4 Y
0.6
0.8
v; 7 0.0
I.0
COORDINATES
0.2
0.4 Y
Fig. 7(c,d)
0.6
COORDINAlES
KLAUS-JORGEN BATHEand JUN DUNG
516 IINA-F
ORIGINAL
:)
K942
XVtllN XVMAX YVMIN YVIIAX
-0.6261 0.000 -0.365$ 0.9494
2 X
$ Y
TEMPERATURE TIME 1
1.000 0.8000
- 0.6000
-- 0.4000
-- 0.2000
- - o.ooor!
Y
COORDINATES
Fig. 7. (a) Natural convection in a 3D cavity. (b) Mesh used (for half the cavity). (c) Velocity distributions in the plane y = 0.5. (d) Velocity distributions in the plane z = 0.5. (e) Temperature distribution; isobands, and in the plane z = 0.5 and x = 0.
and, just as in the analysis of solid mechanics problems, procedures for the automatic construction and refinement of finite element meshes [l]. These and additional topics will provide for some exciting research challenges for the years to come. Acknowledgement-We
appreciate the valuable discussions with our colleague Dr M. H. Wang regarding these developments.
REFERENCES 1. A. Noor (Editor), Computational Mechanics-Adoanees and Trends. AMD, Vol. 75, ASME (1986). 2. P. Bergan, K. J. Bathe and W. Wunderlich (Editors), Finite
Element
Methods
for
Nonlinear
Problems.
Springer, Berlin (1986). 3. K. J. Bathe and J. Dong, Solution of incompressible viscous fluid flow with heat transfer using ADINA-F. Comput. Struct. 26, 17-31 (1987). 4. W. Rohsenow and M. Choi, Heat, Mass and Momentum Transfer. Prentice-Hall, Englewood Cliffs, NJ (1961).
5. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 6. F. Brezzi and K. J. Bathe, Studies of finite element procedures-the inf-sup condition, equivalent forms and applications. In Reliability of Methods for Engineer ing Analysis (Edited by K. J. Bathe and D. R. J. Owen). Pineridge Press, Swansea (1986). I T. Sussman and K. J. Bathe, Studies of finite element procedures--stress band plots and the evaluation of fmiteelement meshes. Engng Comput. 3,178-191(1986). 8. H. Schlichting, Boundary Layer Theory, 7th Edn. McGraw-Hill, New York (1979). 9. W. G. Habashi (Editor), Computational Methoa!s in Viscous Flows, Vol. 3. Pineridge Press, Swansea (1984). 10. M. Haajizadeh, A. Ozguc and C. Tien, Natural convection in a vertical porous enclosure with internal heat generation. Int. J. Heat Mass Transfer 27, 1893-1902 (1984).
11. S. Taneda, Experimental investigation of the wakes behind cylinders and plates at low Reynolds numbers. J. Phys. Sot. Japan 11, 302-307 (1956). 12. G. Batchelor, An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge (1981). 13. J. H. Lienhard, A Heat Transfer Textbook. PrenticeHall, Englewood Cliffs, NJ (1981).
0
1989 Mawell
004.5.7949/89 s3.M) + 0.00 Pqamon Maanillan plc
STUDIES OF FINITE ELEMENT PROCEDURESTHE USE OF ADINA-F IN FLUID FLOW ANALYSES? KLAUS-J~RGENBATHE$ and JIANDON~ $Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A. §ADINA R & D, Inc., Watertown, MA 02172, U.S.A.
Abstract-We present in this paper some experiences obtained with the program ADINA-F in fluid flow analyses. The theory and numerical methods used in the program for laminar and turbulent flow, and flow with heat transfer, are summarized and then the solutions of various problems are presented. These problem solutions comprise the analyses of turbulent flow in a pipe and natural convection in a porous medium, the prediction of flow around a cylinder, the conjugate heat transfer in flow through a pipe and in the cooling of electronic equipment, and the 3D analysis of natural convection in a cavity.
I.INTRODUCTION
flow with heat transfer are, using index notation and the usual summation convention:
The solution of practical fluid flow problems by finite element methods has obtained increasing attention during recent years [l, 21. In particular, the solution of incompressible fluid flow phenomena with or without heat transfer can now be performed efficiently in many practical cases. In an earlier paper, we presented the basic formulation and finite element procedures that we are using for fluid flow analysis [3]. These techniques have been implemented in the computer program ADINA-F, which can be employed for the solution of 2D and 3D fluid flow problems with or without heat transfer. We have continued with the development of this computer program, and the objective of this paper is to present some recent experiences obtained with ADINA-F in the modeling of fluid flow problems. In the paper, we first briefly summarize the basic continuum mechanics equations solved and finite element procedures used. This presentation emphasizes those aspects of the theory that are important when modeling fluid flow problems and interpreting the calculated results. However, the main emphasis of the paper is on the solution of some complex fluid flow problems with ADINA-F. These solutions relate to modeling laminar and turbulent flow, and flow with heat transfer, and demonstrate how the program can be applied in engineering practice.
momentum:
constitutive:
(2) continuity: (3)
(g+vie,i )
(4)
heat transfer:
Pep
=(ke,i),j+q8
where velocity of fluid flow in direction xi mass density, constant components of stress tensor components of body force vector pressure Kronecker delta fluid viscosity, temperature-dependent components of velocity strain tensor = I(vij + Oj,i> specific heat at constant pressure, temperaturedependent temperature thermal conductivity, temperature-dependent rate of heat generated per unit volume, temperature-dependent (frictional effects are included as in Example 3.5)
2. FLUID FLOW EQUATIONS The objective in this section is to briefly summarize the fluid flow equations considered for solution, with an emphasis on the physical and practical aspects that are important when establishing finite element models. 2.1. Continuum mechanics equations Using a Cartesian reference frame (xi, i = 1, 2, 3) the fluid flow equations considered for incompressible
Oi,i =o
and in natural convection problems
t Original figures were generated with a colour producing
fB=dl
terminal and submitted in colour. 499
-iw -43
KLAUS-JURGEN BATHE and JUN DONG
500
where gi
B
0,
acceleration of gravity into direction i coefficient describing the change in mass density per unit change in temperature reference temperature.
The boundary (l)-(4) are:
conditions
corresponding
to eqns
vective terms in eqns (1) and (4), the convection and radiation boundary conditions in eqn (8) and the temperature-dependent material properties. In addition to the nonlinearities listed above, we also include those that arise due to the use of turbulence models (see Sec. 2.6) and the use of non-Newtonian constitutive relations. For solution of eqn (1) it is effective to substitute eqns (2) and (3) into eqn (1) to obtain
prescribed fluid velocities, 4, on the surface S,, ui =
prescribed tractions, ff,
cilS,
p(d) + (3
on the surface S,,
7ijnj=_ffls2
(6)
where nj are the components of the unit normal to the fluid surface and f f are the components of the traction vector; prescribed temperature,
=
-P,i
+ff
+
/iui,#ii.
(12)
Note that the incompressibility condition is now contained in eqn (12). To establish effective finite element models it can also be expedient to consider the nondimensionalized governing equations. If we use the characteristic length L, velocity V, reference temperature &,, reference temperature increment Af3, and
x:
v: = Vi/V,
=x,/L,
P * =p/pv,
8, on the surface S,, 6 =&,
Uj”i,j)
f”
e*,e-eO
(7)
--s-
t* = tV/L
=f tgL/pV qBL2
q*=k
(13)
*e
prescribed heat flux into the surface S,,
(8)
we obtain as the governing momentum equations, ti: +u,?vit,=
where qs is the heat flux input to the body. We note that for eqns (S)-(8) S = S,, fl= S,nS2 and S = m, 8 = S,nS,, where S is the total surface area of the fluid body. The heat flow input in eqn (8) comprises, firstly, the effect of actually applied distributed heat flow; secondly, the effect of convection heat transfer, in which case qS=h(e,-eS)
(9)
where h is the convection coefficient, and 6, and fIs are the environmental and body surface temperatures; and thirdly, the effect of radiation heat transfer, in which case qs =
K(e, - es)
K = h,[(w2 + (eS)%e, + 6’)
(10) (11)
and h, is the radiation coefficient and 6, is the temperature of the radiating body [4]. The value of h, is determined from the Stefan-Boltzmann constant, the emissivity of the radiating and absorbing materials and the geometric view factors. The above equations are the standard Navier-Stokes equations governing the motion of a viscous incompressible fluid in laminar flow with heat transfer. Inherent nonlinearities are due to the con-
-p:
and energy
+fY +(lIRe)v&
(14)
and
Pe(B*+O:e,*,)=e:i+q*
(15)
where the Reynolds number is
Re.!%; v=!.! V
P
and the Peclet number is
pe=yL; c,=k. a
PCP
The use of the nondimensionalized equations in the finite element solution is achieved by simply using the following material properties: mass density = unit value 1 . v1scosny = G heat capacity = Pe conductivity = unit value.
(16)
Of course, in the above equations we implicitly assumed that the fluid properties are constant. If these properties vary with temperature, eqns (14) and
ADINA-F in fluid fiow analysis
(15) are still applicable but the properties to be used in these equations are then tem~ratu~~e~dent. Note that we may also interpret eqns (14) and (15) as the governing fluid flow equations with the material properties in eqn (16) in any consistent set of units. When using eqns (13), (14) and (15) for solution of a fluid flow problem, it is of course also necessary to impose the velocity, traction and heat transfer boundary conditions in nondimensionalized form. These boundary conditions are a ~statement of eqns (5+0-(l) using the nondimensionaiized variables given above. For the momentum equation we use v’t to impose boundary velocities, and to impose tractions we recognize that
zij*=
-p*a,+j$v~j+v$)*(17)
For the energy equation we use the values 8* to prescribe temperatures and qs* to impose heat flux conditions,
B
s+ _
qSL
-z-z
The appropriate boundary convection and radiation conditions are also easily derived.
501
heat transfer:
PCJ@ + v,B,,)68 dk’ +
k6,, 60,, d ?’ sV
+
S@K(~, - es) dS
(21)
s&
where S,, S, and S, are part of S, and correspond, respectively, to the surfaces subjected to heat flow input, convection and radiation conditions. Consider next a single finite element for solution of eqns (19x21). The finite element equations for an assemblage of elements are obtained by the usual direct addition of element matrices [S]. For the finite element discretization we use for 22) (planar and axisymmetric) solutions the 9/4 node element shown in Fig. 1: the nine nodes are used to interpolate the geometry, the velocities and the temperature
2.2. Finite element equations We presented the finite element formulation used in detail in [3]. Here we summarize briefly the major ingredients of the formulation and resulting equations and point out some important features regarding the practical application of these equations. The finite element solution of the continuum mechanics equations given above is obtained by establishing a weak form of the equations using the Galerkin procedure [5]. The momentum equations are weighted with the velocities, the continuity equation is weighted with pressure and the energy equation is weighted with the temperature. Using integration by parts gives the variational equations to be discretized by finite element inte~lations: momentum:
=S;f:d~jdV+~sJ’6v;dS
e= i
h&e&;
(24)
k-1
the four nodes are used to interpolate
P= i
k-l
the pressure
$Pk.
Correspondingly, we use in 3D analysis the 27/8 node element (see Fig. 1) for a triquadratic interpolation of the geometry, velocities and temperature and a trilinear interpolation of the pressure. These interpolations provide stable elements as expressed by the inf-sup condition of Brezzi and Babugka [6]. The finite element governing equations are now established as
(19)
M,.* + (K,,,+ Km. Iv + K, p = R, + Rs where 6 means ‘variation in’ (or a virtual quantity); continuity:
K&v=0 Ce+K,e+K,~,~v=Q,+Q,+Q,+Q,
s
Spv,,,dV=O
V
(25)
(20)
(26) (27)
(281
where the coefficient matrices and load vectors are dire&y obtained by substituting the finite element
502
KLAU-JORGEN Bnrns and
1
W&
Fig. 1. 2D and 3D elements available in ADINA-F.
interpolations [see eqns (22)-(25)] into the variational eqns (19H21) with eqns (2) and (3). The details of the matrices in eqns (26x28) are given in [3]. 2.3. Boundary conditions in finite element solutions
The boundary conditions to a problem governed by eqns (l)-(4) are those listed in eqns (5x8). Here we note that eqns (5) and (7) are the essential boundary conditions and eqns (6) and (8) are the natural boundary conditions. The essential boundary conditions are imposed in the finite element analysis by directly setting the solution variables in eqns (26) and (28) equal to the given values. Note that the velocities may be prescribed in a coordinate system different from the Cartesian system, which then requires a standard transformation [S]. The natural boundary conditions correspond to prescribed tractions ff and prescribed heat fluxes qs (which may be given in terms of convection or radiation conditions). As for prescribed velocities, the imposed tractions can be defined in any coordinate system, which again only requires a standard transformation Considering the finite element solution, the natural boundary conditions of eqn (8) are imposed by the right-hand-side vectors in eqn (28). However, as for imposing the natural boundary conditions of eqn (6), consider the vector R, given for an element side as RS=
HSTfs dS
s .%
(29)
where HS is the standard element surface displacement interpolation matrix and B is the vector of the imposed force effects [3]. In the solution F is given
JIAN DONG
and substituted into eqn (29). Note that in general this does not mean that the pressure on the element side is prescribed. Further, we recall that since eqn (26) contains the continuity equation, the components in ‘fs are in general not given by eqn (6) [3]. The physical tractions would be used in F if the continuity equation were not imposed onto the momentum equation. In analyses there is little difference in the results using either approach unless an actual traction is prescribed (and for this case we have in ADINA-F a special option of solution). However, we have also found that using eqn (26) the convergence in the iterative solution is sometimes better than when not imposing the continuity condition onto the momentum equation. Since we are solving incompressible flow problems, particular attention needs to be given to the boundary condition on the pressure. Namely, the pressure is undetermined within a constant unless prescribed at, at least, one point. In the finite element solution we have that: all pressure degrees of freedom in eqn (26) can be left free if the component in F normal to the boundary (which involves the pressure) is prescribed (to zero or any other value) at least at one nodal point on the boundary; if this component is not prescribed anywhere along the boundary, then it is necessary to arbitrarily assign a pressure value at one node in the finite element mesh. The need to prescribe a pressure value at one node in the finite element mesh arises therefore when closed fluid regions are analyzed for which the velocities are prescribed at all points of the boundary (see Sec. 3). 2.4. Time stepping and iterative equation solution With our solution procedures we consider steadystate analyses and transient analyses. In steady-state solutions, the terms M,ir and Cd in eqns (26) and (28) are neglected, and the time stepping merely corresponds to an iterative incremental analysis of the nonlinear response (51. The number of solution steps are selected to arrive at an efficient analysis and to obtain the solution response at specific levels of prescribed loadings and boundary conditions. In transient analyses, we use the a-method with f < c( < 1 with the assumptions ‘i.b’+~+A$_$,)
(30)
~++=(l
(31)
_._~~()‘v+~‘+*‘v
f+zAf(j='(,+A,
_'e)
At
r+xAq = (1 _ a)‘@+ E’+A’g
(32) (33)
ADINA-F in fluid flow analysis and ‘+*A’p= (1 - a)‘p + a’+Alp.
(34)
Usually, we use a = $ {trapezoidal rule) or a = 1 (Euler backward method). A transient analysis can be performed with variable time step sizes (as input by the user) and in steadystate and transient solutions an automatic-time-stepping (ATS) method can be used. In the ATS method, the program decreases automatically the load or time step sizes from those input by the user if convergent in the iterations is not reached for the current load or time considered or if certain accuracy criteria are not satisfied. Finally, we should mention that the incremental solution using eqns (26)--(28) is performed by the method of successive substitution or the full Newton-Raphson iteration (with or without line searches). In cases with strong nonlinea~ties, the Newton-Raphson iteration is more effective. 2.5. Calculated quantities-the stream faction values
use of vorticity and
The solution of the governing finite element equations [eqns (26)-(28)] yields directly the velocities, pressure and temperature over the fluid domain. All these quantities are continuous over the element boundaries. As additional quantities it can be useful to calculate the vorticity, w, of the fluid o=vxv.
(35)
In 2D analysis the stream function values, $, can also be useful. These are calculated by the line integral
$ -$a=
I
(v,h-vsW
(36)
where & is a constant. We note that the vorticity of the fluid is obtained by differentiation of the calculated velocity field, hence for coarse finite element models discontinuities across element boundaries will be observed. Therefore, we can use the isobar& of vorticity to identify whether a finite element mesh was fine enough to predict an unknown flow field accurately [7j. The stream function values are obtained by integrating the calculated velocity field. Here we note that the stream function at a point is only unique if the continuity condition, eqn (3), is satisfied. Hence, if the calculated solution does not satisfy the continuity ~ndition accurately, different values for $ will be obtained depending on which path of integration is chosen. This observation shows that the calculated values of the stream function can also be employed to assess whether a finite element mesh was fine enough to predict a flow field.
503
In the ADINA-F implementation we integrate from a starting point (a finite element node) along the element sides and around the element boundaries, always only along those element boundaries which attach to nodes for which rj has not yet been calculated. Hence, the integration is only performed along an element boundary once and only if it leads to the calculation of a yet unknown value of JI at a node on that boundary. Using this algorithm we can conclude that, on display of the calculated stream function values, nonsmooth lines indicate regions where the mesh was not fine enough to predict the fluid fIow accurately. 2.6. Turbulence modeling One of the most difficult aspects in many fluid flow problems is the modeling of the turbulent behavior of the fluid. The field of modeling turbulence has attracted much attention and various models have been proposed and are in current use [S, 91. We have implemented in ADINA-F the use of the mixing length model and the large-eddy-simulation (LES) model. For the mixing length model, the users can define by their own prog~mming the specific equations for the mixing length 1 to be used. The general form of the model is given by It,=P~21JI p-1
(37)
* -i dul+duI ax, axi
( )
where p, is the turbulence (eddy) viscosity to be added to the laminar viscosity, A, to obtain the total viscosity p. For the LES model we use
where k, is a dimensionless constant and A is taken to be the mean of the side lengths of the element for which P, is being computed. An application of the mixing length model is presented in Sec. 3.1.
3. EXAMPLE SOLUTIONS
The following solution results have been obtained with ADINA-F, where we should note that no specific attempt was made to use optimal meshes or optimal solution approaches with the program. The solution results are presented using the postprocessor ADINA-PLOT. 3.1. Analysis of turbulent flow in a pipe Figure 2(a) shows the pipe considered. The flow is to be solved for a high Reynolds number, so that turbulence effects are important. For the analysis a mesh of 20 x 1 axisymmetric elements was used with
KLAUS-JORGEN BATHE and
504
JUN
DONG
one layer to model a typical section of the fluid in the pipe. To include the turbulence effects, the mixing length model proposed by Nikuradse [8] was employed 1 = 0.14 - 0.08~~ - 0.06~~
(40)
where 2 = i/R, q = r/R, 1 is the Prandtl mixing length, R is the radius of the pipe and r is measured from the centerline of the pipe. This model is adequately used for the flow not very close to the wall. To include the near-wall effects we used for the element at the wall the Nikuradse formula with the van Driest correction 1 = (0.14 - 0.08r~’ - 0.06~~) x
1 -exp
-“U,
rl)R
(
0 NIKURADSE A AOINA-F. X ADINA-F.
EXPERIMENT. RE=3.24E4 RE=2.3E4
RE=2.3E4
(4,) )
where r7= (TO/~)“~, v = h/p, k is the laminar viscosity and 70 is the shear stress at the wall. We programmed the above relations into the usersupplied subroutine of ADINA-F to define the mixing length I. In this case, the application of the mixing length model gives for the total viscosity
Figure 2 shows the solution obtained for Reynolds numbers of 2.3 x lo4 and 3.24 x 104. These numbers are based on the average fluid velocity in the pipe.
-8.0
Heat genwatio~ in
ONIKURADSE A ADINA-F. X AOINA-F.
DISTANCE
EXPERIMENT RE-3.24E4 RE=2.3E4
porous medium
8
FROM WALL
Fig. 2. (a) Analysis of turbulent flow in a pipe. Twenty equal sized elements into radial direction are used to model the fluid. (b) Velocity distribution in the pipe. (c) Virtual viscosity distribution in the pipe.
w-0,
Insulated
c----l.0
Fig. 3(a)
qB
8-O
505
ADINA-F in fluid flow analysis
1EtlPERATURE TlnE 1.000 - 0.1400
lEtlPERATURE TlHE 1.000 - 0.1200
(b)
0.1200 0.0900
-
0.0900
-
0.0400
-
0.0100
- 0.0600
- 0.0300
- 0.0100
RamlO*
TEMPERATURE TIME 1.000 - 0.02500
TEMPERATURE TINE 1.000 - O.llOO
0.0900
-
0.0600
-
0.0300
-
0.0100
0.02000
- 0.01500
-
*Ra.lO'
RarlO'
Ra.105
(d
+ ADINA-F. RA=IDO - ANALYTICAL VERTICAL
CENTEhNE
2. COORDINATE
VERllCAL
CENlERLlNE
COORDlNAlE
-
Fig. 3. (a) Natural convection in a porous medium. (b) Streamlines for Merent values of the Darcy-Rayleigh number. (c) Temperature contours for different values of the Darcy-Rayleigh number. (d) Velocity along the vertkal centerline. (e) Temperature along the vertical centerline.
0.01000
- 0.00500
506
KLAUS-JORGEN
BATHE and
Figure 2(b) shows the velocity distribution and Fig. 2(c) gives the nondimensional&d virtual viscosity defined as
v*=--- P
pGR ’
The figures show good correspondence with the experimental results of Nikuradse. Note, however, that the velocity prediction near the wall is still quite approximate due to the very high velocity gradient. Also, we should note that the solution for the virtual viscosity is given at the element integration points. 3.2. Natural convection in a porous enclosure with internal heat generation
Figure 3(a) shows the problem considered. For this class of problems the governing continuum mechanics equations are vi., = 0
P
KVi =
-P,i
(PJp)fi +
+
JIAN LANG
3.3. Wake behind a cylinder at low Reynolh numbers Figure 4(a) shows the problem considered and schematically the twin vortices and trail in the wake of the cylinder. Our objective was to predict these vortices and the pressure distribution on the cylinder for a range of low Reynolds numbers. Figure 4(b) presents the two finite element meshes that we have used. The analyses with both meshes gave almost the same streamlines and vorticity results; however, mesh A was used to obtain a better pressure prediction. Figure 4(c) shows the streamlines calculated for Re = 40. We clearly identify the twin vortices, which are also shown in Fig. 4(d) by the plot of the vorticity. Figure 4(e) gives a comparison of the predicted wake lengths for different Reynolds numbers with the experimental results of Taneda [ll], and Fig. 4(f) shows the pressure distributions on the cylinder, as observed by Thorn [12] and in this finite element analysis. The nondimensionalized pressure is defined as
Pgi(l - BCe- eO))
$PcpVje,i
where C$is the porosity,
(l;81i),i + 4’
=
K
*_P-PO P -1 TPVO
where p. is the inflow pressure and V. is the inflow velocity.
is the permeability,
BEp= 4PCp + (1- d)(PCJS E=4k+(l--i$)k,
(46) (47)
and fi is the viscosity of the fluid accounting for porosity. The subscript s is used to identify the solid matrix properties, and properties without a subscript are those of the fluid. The above equations are directly obtained from eqns (l)-(4) by including Darcy’s law with the Boussinesq approximation, and hence are an application of the general equations solved in ADINA-F. This problem is characterized by a Darcy-Rayleigh number
3.4. Conjugate heat transfer in a pipe withJluidJow Figures 5(a) and (b) show the problem considered and the finite element discretization used for the analysis. The objective is to study the temperature variation in the fluid and the wall. As a first case, we solved the problem when the heat flux is imposed along the complete pipe section considered and the fluid flow is a fully developed laminar flow. The solution not near the entrance and the exit sections should be in accordance with the equation [ 131,
where Ob= bulk temperature where L is the width of the enclosure and 6 is the effective thermal diffusivity. A 20 x 40 element mesh (with equal sized elements) was used to model the domain of the enclosure. Figures 3(b-e) give the calculated steady-state solutions. These correspond to the results given by Haajizadeh et al. [lo]. We note that the solution for Ra = lo5 was obtained by adding the term p,,(vy + u,,,),, with ~,/p = 0.001 to the right-hand-side of eqn (44) in order to obtain faster convergence in the iterations. Figures 3(d) and (e) shows the vertical velocity and temperature distributions along the vertical centerline for the case of Ra = 100.
v = 2v,,,
j Ovr dr of fluid = Ivr dr
r*
( > 1- -
R*
c(= thermal diffusivity of fluid qi = interface heat flux V aYe = average velocity of fluid k, = fluid conductivity and z is the axial coordinate of the pipe. Figure 5(c) shows wall, pipe axis and bulk temperatures of the fluid (here we have tl = 0.2, qi = 1.l, u.,~ = 0.5, k,= l.O), and shows that the correct variations of temperature along the pipe are predicted.
507
ADINA-F in fluid flow analysis Next, the case of Fig. S(a) was considered, in which the heat flux is imposed on the pipe wall for a finite length. This type of problem is usually characterized by the Peclet number, Pe, and the ratio of the solid and fluid conductivities, K, where Pe = Re*Pr K = k,/k,
24, K Re (Reynolds number) = V
Pr (Prandtl number) = f
Reynolds number is used for each solution. The Eckert numbers based on the maximum temperature differences in the fluid are calculated to be 1.5 and 0.5, respectively. As expected, the figures show that, as the Peclet number increases, the heat transfer by convection becomes more predominant. The finer mesh (mesh B) gives a better streamline solution, and therefore a better temperature prediction. Note that in this analysis the solution for the velocities and temperatures were obtained simultaneously, where actually the solution could be calculated by first solving for the velocities and then for the temperatures. 3.6. Natural convection in a 30 cavity
v = fluid kinematic viscosity k, = solid conductivity k, = fluid conductivity.
The results of the bulk temperature distributions and interface heat fluxes along the pipe for different Peclet numbers and values of K are shown in Figs. S(d-g). 3.5. Conjugate heat transfer in cooling of electronic equipment The problem considered here was mainly chosen to investigate some of the interesting features available in our fluid flow analysis capabilities. We consider a channel of fluid flow with two equal obstructing bodies (the electronic equipment) in its flow that generate heat. The heat is transported in the fluid by conduction and convection. We have no comparison with experimental or other analytical results available, but want to show the predicted temperature distributions for different flow conditions. The problem and finite element meshes used are shown in Figs. 6(a) and (b). The meshes A and B were used in the solution for comparison purposes. We note that the flow is prescribed at the left entrance of the channel. The solid bodies generating heat are modeled with fluid elements that only have the temperature degree of freedom (velocities and pressure are set to zero). The conductivity of the solid is set to be the same as that of the the fluid. The governing equations of the fluid flow considered are those discussed in Sec. 2 with the heat generation due to friction included @ = 2pe,e,.
(49)
The effect of this term is discussed extensively in [8] and is characterized by the Eckert number, EC. The complete problem is characterized by the Peclet, Eckert and Reynolds numbers. Figures 6&-f) give the calculated temperature distributions for different Peclet numbers (10 and 100) and the streamlines of the flow corresponding to a Reynolds number 100, for meshes A and B. The same C.A.S. 32,34-B
The solution of a natural convection problem considering two-dimensional motion was presented by many authors (see also [3]). To study a threedimensional motion, we considered the problem shown in Fig. 7(a). No slip conditions were assumed at all surfaces. The temperature was prescribed at the side surfaces and adiabatic conditions were assumed at the other surfaces. The governing equations are those discussed in Sec. 2 but including the Boussinesq term, pg(l - /?(0 - I$,)), in the z-momentum equation. We required the solution for Ra = lo4 and Pr = 0.2. For the finite element analysis, we considered only one half of the domain due to symmetry conditions with respect to the plane x = 0. The mesh used is shown in Fig. 7(b). Figures 7(c) and (d) give velocity distributions at the planes y = 0.5 and z = 0.5, respectively, for different values of x = constant. As expected the y-velocity field shown in Fig. 7(c) is similar in shape to the velocity field in a two-dimensional analysis (31, but the magnitude of the velocities, of course, decreases as x increases. Similarly, for the z-velocity field given in Fig. 7(d). An interesting result is obtained for the x-velocity field, which shows recirculation near the edges. Figure 7(e) presents the calculated temperature field. It is seen that for this problem the temperature distribution is still very nearly two-dimensional. 4. CONCLUDING REMARKS
The objective in this paper was to summarize some of the latest developments we have undertaken for finite element analysis of fluid flow. We first summarized the theory and numerical methods used, and then presented some example solutions. It can be concluded that ADINA-F is a powerful tool for finite element solution of a variety of fluid flow problems. However, considering analyses of incompressible fluid flows in general, there are of course still a number of important topics for research and development; these include the development of a robust algorithm to predict free surfaces, the development of new more efficient iterative solution schemes,
508
KLAUS-J~~RQEN BATHE and JIAN DONO
t. w
I
T v.v, _, -
z-
(0)
(b)
STREAH FUHCT IOH TIME 1.000 1.000
Fig. 4(a-c)
-
-0.500
-
-1 .ooo
OMEGA-YZ
TIRE i.OOC
(d)
OEXPERIHENTAL BY TANEDA A ADINA-F. tlESH A X ADINA-F. HEW B
REYNOLDS
1956
NUMBER
*. i-
0.
DEGREES
I
40.
FRDII
FD:iARG
*
S~~NA,*O~O~N~
Fig. 4. (a) Flow past a cylinder, Re = pV,d/~, schematic representation. (b) Top, mesh A, bottom, mesh B. (c) Streamlines at Re = 40, mesh A. (d) Isobands of vorticity at Re = 40, mesh A. (e} Lengths of vortioes for different Reynolds numbers. (f) Pressure distribution around cylinder. 509
KLAUS-J~RGENBATHEand JMN DONG
510
t= 0.1
(a)
RI 1.0
FIVE FLUID RADIALLY
ONE SOL IO RADIALLY
+ WALL 0 BULK X PIPE
“07
Id.
5.
AXIAL
Ii.
ELEMENTS
ELEMENT
AXIS
0 PE-5 A PE=S + PE-5
K=I K=5 K-10
2d.
DISTANCE
AXIAL
Fig. 5(a-d)
DISTANCE
511
ADINA-F in fluid flow analysis
(e)
0 PE=50 A PE-50 + PE=SO
AXIAL
0PE=5 A PE-5 + PE-5
K=lO K-50 K=lOO
DISTANCE
AXIAL
0 PE-50 A PE-50 + PE=50
AXIAL
K=l K=5 K=IO
DISTANCE
K=IO K-50 K=lOO
DISTANCE
Fig. 5. (a) Conjugate heat transfer in an “infinitely long” pipe. with fluid flow. A section of length 20 is modeled. (b) Mesh used. (c) Temperatures along the pipe for applied heat flux along complete pipe. (d) Interface heat flux, Pe = 5, K = 1,s and 10. (e) Interface heat flux, Pe = 50, K = 10,SO and 100. (f) Fluid bulk temperature, Pe = 5, K = 1.5 and 10. (g) Fluid bulk temperature, Pe = 50, K = lo,50 and 100. Note that for the case Pe = 50, a longer section pipe should be modeled.
512
KLAUS-JORGENBATHEand JUN DONG
:bl
TEMPERATURE TIME 1.000 ~5.000 -
4.000
--
3.000
--
2.000
--
I .ooo
Pe = 10 TEMPERATURE TIME 1.000 J_
I .700
T
I .500
0.900
t Pe = 100 Fig. 6(a-c)
0.600
513
ADINA-F in fluid flow analysis TEMPERATURE TIME 1.000 l.
(d)
5.000
4.000
I 3.000
2.000
I .ooo
Pe=
10
TEMPERATURE TIME 1.000 1
I .700
-
1 *so0
--
1.200
--
0.900
- - 0.600 Pe = 100 STREAM FUNCTlON TIME 1.000 - -0.1000
(e)
s
-
-0.3000
-
-0.5000
-
-0.7000
-0.8000
Fig. 6(d,e)
KLAUS-J~~RGEN BATHEand JIAN DONG
514
STREAM FUNCTION TIME 1.000 - -0.200
(f)
-
-0.500
-
-0.800
Fig. 6. (a) Analysis of conjugate heat transfer in cooling of electronic equipment. (b) Top, mesh bottom, mesh B. (c) Predicted temperatures using mesh A. (d) Predicted temperatures using mesh B. Streamlines using mesh A, Re = 100. (f) Streamlines using mesh B, Re = 100. top surface, insulated (a)
back surface, insulated
--.
cdd side surface,&
front surface, ins&ted
L, - L,d.O L, -0.4
hot side surface, oh
L
,’
/
/ /’
/
/
insulated
(b)
Fig. 7(a,b)
515
ADINA-F in fluid flow analysis
0 AT x-o.02 A
AT
+A1
x-o.10
XAT 0
+
AT
X-O.10
Q
X=0.06 x=0.14
AT
X=0.18
d 7
v; 7 0.0
0.2
:ObRDINATES . 0.6
0.8
I .o
0.8
I.0
2
d) v;
0 AT x-o.02
d 4’
AT
A
AT
X=0.06
+
AT
x=0.10
X
AT
x=0.14
0
AT
X=0.18
X.0.0
X
AT
X-O.14
0
AT
X=0.18
v; * Cm 8’0 _I-& y* it
0.0
0.2
0.4 Y
0.6
0.8
v; 7 0.0
I.0
COORDINATES
0.2
0.4 Y
Fig. 7(c,d)
0.6
COORDINAlES
KLAUS-JORGEN BATHEand JUN DUNG
516 IINA-F
ORIGINAL
:)
K942
XVtllN XVMAX YVMIN YVIIAX
-0.6261 0.000 -0.365$ 0.9494
2 X
$ Y
TEMPERATURE TIME 1
1.000 0.8000
- 0.6000
-- 0.4000
-- 0.2000
- - o.ooor!
Y
COORDINATES
Fig. 7. (a) Natural convection in a 3D cavity. (b) Mesh used (for half the cavity). (c) Velocity distributions in the plane y = 0.5. (d) Velocity distributions in the plane z = 0.5. (e) Temperature distribution; isobands, and in the plane z = 0.5 and x = 0.
and, just as in the analysis of solid mechanics problems, procedures for the automatic construction and refinement of finite element meshes [l]. These and additional topics will provide for some exciting research challenges for the years to come. Acknowledgement-We
appreciate the valuable discussions with our colleague Dr M. H. Wang regarding these developments.
REFERENCES 1. A. Noor (Editor), Computational Mechanics-Adoanees and Trends. AMD, Vol. 75, ASME (1986). 2. P. Bergan, K. J. Bathe and W. Wunderlich (Editors), Finite
Element
Methods
for
Nonlinear
Problems.
Springer, Berlin (1986). 3. K. J. Bathe and J. Dong, Solution of incompressible viscous fluid flow with heat transfer using ADINA-F. Comput. Struct. 26, 17-31 (1987). 4. W. Rohsenow and M. Choi, Heat, Mass and Momentum Transfer. Prentice-Hall, Englewood Cliffs, NJ (1961).
5. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 6. F. Brezzi and K. J. Bathe, Studies of finite element procedures-the inf-sup condition, equivalent forms and applications. In Reliability of Methods for Engineer ing Analysis (Edited by K. J. Bathe and D. R. J. Owen). Pineridge Press, Swansea (1986). I T. Sussman and K. J. Bathe, Studies of finite element procedures--stress band plots and the evaluation of fmiteelement meshes. Engng Comput. 3,178-191(1986). 8. H. Schlichting, Boundary Layer Theory, 7th Edn. McGraw-Hill, New York (1979). 9. W. G. Habashi (Editor), Computational Methoa!s in Viscous Flows, Vol. 3. Pineridge Press, Swansea (1984). 10. M. Haajizadeh, A. Ozguc and C. Tien, Natural convection in a vertical porous enclosure with internal heat generation. Int. J. Heat Mass Transfer 27, 1893-1902 (1984).
11. S. Taneda, Experimental investigation of the wakes behind cylinders and plates at low Reynolds numbers. J. Phys. Sot. Japan 11, 302-307 (1956). 12. G. Batchelor, An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge (1981). 13. J. H. Lienhard, A Heat Transfer Textbook. PrenticeHall, Englewood Cliffs, NJ (1981).