Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 519-526 UDK - UDC 536.2 Izvirni znanstveni clanek - Original scientific paper (1.01) Heat transfer distribution for a free/porous system with forced convection and heat generation - a numerical study Antonio C.M. Sousa 1,2* 1 1Departamento de Engenharia Mecânica; Universidade de Aveiro; Campus Universitário de Santiago; 3810-193 Aveiro-Portugal. 2Department of Mechanical Engineering; University of New Brunswick; Fredericton, NB, Canada E3B 5A3,*asousa@unb.ca Abstract This paper reports on a numerical study for steady flow and heat transfer distribution for a configuration relevant to Liquid Composite Molding, where a gap between a porous substrate and the solid boundary of a mold cavity yields an edge flow. The flow within the porous domain is modeled by the Brinkman-Forchheimer formulation, and the edge flow itself is described by the Navier-Stokes equations. The cure of the fluid (resin) is simulated as a volumetric heat generation. The predictions are obtained using a well-tested control-volume finite element method, however, a novel methodology had to be devised to define the interface between the free and porous system. The most relevant finding is the critical role of the gap upon the quality of the part. The presence of the gap can reduce substantially the average flow through the porous substrate, therefore yielding high temperature levels in this region. These temperatures may be sufficiently high to cause serious defects to the part being molded. Introduction Numerical and analytical investigations of flow and heat transfer distribution in composite systems containing simultaneously a porous and an open fluid domain are receiving renewed interest from the scientific and engineering community due to the demand for the development of models that can enhance mold design in the area of Liquid Composite Molding (LCM)[1,2]. LCM processes, such as Resin Transfer Molding (RTM) are preferred manufacturing processes for large structural components of complicated shape made out of polymer composites. These processes, in general, are characterized by their high cost-effectiveness, relatively simple tooling requirements, low cycle times, and net-shape production. They require the impregnation of a polymeric resin through a porous preform, which is placed in the mold cavity, and it can be composed of glass, carbon, or Kevlar fibres. The main difficulty with these processes is usually associated with the eventual presence of small clearances between the preform and the mold edges, which result from rough cutting, ill fitting, or deformation of the preform. The clearance yields a preferential flow path, which can disrupt the filling of the mold and the impregnation of the preform, resulting in poor quality of the part due to voids, residual stresses, and poor bonding between the fibres and the resin. This preferential flow, which involves the interface between a porous and non- porous medium, poses a major phenomenological challenge, and it has been the focus of intense research over the years. In the past, most models for edge flows used the Darcy’s Law, and they were formulated either on the basis of an analytical formulation for the interface [3], e.g. [1,2], or on the basis of an equivalent permeability for the edge [1,2,4]. Both approaches, however, do not take into account the transverse flow, which can be important, and based on analyses using a transverse flow factor, as discussed in [2], often, it cannot be neglected. Moreover, although the Darcy’s Law with appropriate modifications [5] can successfully replicate the flow features for random fibre mats as well as woven fibre mats, significant deviations from the Darcy’s Law predictions were documented [6]. To overcome these shortcomings a state-of-the-art model for the flow [7] was developed, and it is based on the Navier-Stokes equations combined with the Brinkman-Forchheimer equations. In the present work, the above-mentioned model [7] is further enhanced by taking into consideration thermal effects. For this purpose a heat transfer submodel is built in. The cure process [8], which yields an exothermic reaction, is modeled, although in a somewhat simplified form, as a steady, uniform volumetric heat generation in the fluid. Control of this exothermic reaction is critical in what concerns the component's final structural integrity and quality. The heat dissipation is investigated in terms of a non- 519 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 Nomenclature y cartesian co-ordinate, m Da Darcy number, ND Greek symbols F Forchheimer coefficient, ND a thermal diffusivity, m2s-1 H channel height (reference dimension), m Ct1 viscosity ratio, ND k thermal conductivity, Wm-1 *ij Kronecker delta, ND L channel length, m M temperature difference, K m mass flow rate, kg.s-1 e porosity, ND n unit normal vector, ND K permeability, m2 p pressure, N.m-2 M dynamic viscosity, kg.s-1.m-1 Pr Prandtl number, ND V kinematic viscosity, m2.s-1 q" heat flux, Wm-2 fiB Brinkman modified effective viscosity, kg.s-.m-1 q" heat generation, Wm-3 p density, kg.m-3 Rc ratio of thermal conductivities, ND a total stress, N.m-2 Re Reynolds number, ND t unit tangent vector, ND Subscripts T temperature, K B Brinkman value U cartesian velocity component in the ch channel x-direction, m.s-1 f "free" fluid domain u reference velocity, m.s-1 i interface y-direction, m.s-1 i,j cartesian subscripts ; tensorial notation v cartesian velocity component in the in inlet y-direction, m.s- fd free domain V volume, m3 pd porous domain V velocity vector, m.s-1 ref reference value * normalized variable s solid x cartesian co-ordinate, m * normalized variable dimensional heat generation, Reynolds number, Darcy number, and the fluid and porous medium properties. The influence of the gap, and in particular its dimensions, upon the heat dissipation from the core region of the porous medium is also investigated. The present model, in what concerns the flow modeling, as already mentioned, follows closely [7], and it can be succinctly described as follows: 1. The flow within the porous medium is governed by the Brinkman-Forchheimer model; 2. The Navier-Stokes equations are used in the modeling of the edge flow, i.e. the “free” medium. This approach, by combining the Navier-Stokes equations with the Brinkman-Forchheimer equations, avoids the well-known difficulties associated with linking the Navier-Stokes and the Darcy equations [9,10]; 3. At the “free”/porous medium interface, the only approximation required is to assume that the fluid fully supports the tangential and normal stresses; and 4. The cure is modeled by assuming a steady, uniform volumetric heat generation within the fluid. For the sake of completeness, an overview of the model's development is presented. Physical and numerical modelling Physical domain and geometry The two-dimensional configuration under analysis is presented in Fig. 1. The fluid enters the mold cavity from the left, with a constant velocity, flows through the porous medium, (the shadowed region of Fig. 1) or through the top and right-hand side channels (edge flow), leaving the mold cavity through its right-hand side. In the exit region the outflow conditions are placed further downstream to reduce their influence upon the upstream flow. The flow and energy transfer are assumed to be steady and laminar. The fluid physical properties are taken as constant. Figure 1. Domain and geometry H=25mm; L/H=5; D/H=0.1). (Dimensions: Free fluid region The flow in the open region is governed by the continuity and Navier-Stokes equations. Defining the dimensionless variables: x*i-xi/H, u*i=ui/U, p* = p/(pU2), and Re = pUH//u, the dimensionless version of the governing equations, using the indicial notation (i,j = 1,2), becomes: 520 Sousa A.C.M. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 dx*i (u*i) = 0 dx I u*iu* j I — 9 dp* -------------- Cx* j Ox 1 du* j Re dx (1) (2) Taking the reference temperature difference as ATref = RePrU2/cP , the Prandtl number, Pr = v/a and the dimensionless heat generation the energy conservation equation in the open domain becomes: q , g , ,*=q , g ,,H [paU 9 t T\ 9 dx dx RePr dx 1 dT ) q*g (3) (RePr) Since the ‘artificial’ exit region was considered of small length it is assumed to be a region with no volumetric heat generation. Fluid-saturated porous region The dimensionless mass conservation equation in the porous medium is exactly the same as for the open fluid domain, (Eq. (1)). The fluid flow in the porous domain is governed by the Brinkman-Forchheimer equations, [11,12], that give, in dimensionless form: 9 (1 -------- —u *i u *j dx \e «1 dp* 9 dx* j dx Res dx V + 9u (4) /Da ReDa^ where Ł is the permeability of the porous medium, «1 = ^ /ju is the ratio between the Brinkman viscosity and the fluid viscosity, f = 1.75/v150Ł5 is the Forchheimer coefficient, Da = K/H2 is the Darcy number, and |V*| is the modulus of the dimensionless velocity vector. The thermal conductivity of the porous medium is given by the combination of the porous solid matrix and saturating fluid thermal conductivities. In this work it is used a geometric mean approach [12], given by k = kŁfks 1 ~Ł' =kfRc1 , where Rc is the ratio of thermal conductivities,Rc =ksIkf . The dimensionless energy conservation equation in the porous domain is therefore given by: Rc1~e^ dT 5 I T \- ° ---------1 u*i* I —--------- dx dx RePr dx qg*S (RePrf (5) In the porous domain, only the volume fraction s is filled with liquid, where heat is being generated. Interface At the interface of the fluid-saturated porous and open domains, mass, momentum, and energy balances must be satisfied as discussed in [10]. It is assumed here that there are no phase change phenomena or chemical reactions, and local thermodynamic equilibrium prevails. Mass conservation at each point of the interface implies that, in dimensionless form, (V*.n*)fd = (V*.n*) d (6) where the subscripts fd and pd indicate free domain and porous domain values, respectively, and n* is a unit vector normal to the interface. Taking the dimensionless stress CT*ij = (1Re)(du*i /dx* j + du* j /dx*i I - p*ij , the continuity of the normal and shear stresses is (n*.n*.o* J fd — (n*.n*.o* J d f t*.n*.o* Jfd = (t*.n*.o* J d (7) t* being the dimensionless unit vector tangent to the interface. In addition, it is assumed that the pressure is continuous across the interface, that is, IpA =(p*) (8) The mass conservation equation, Eq. (6), requires equal normal velocities on both open- and porous-domain sides. The essential velocity continuity condition at the interface [13], however, implies: C V*) = (V*) (9) Equation (7) requires that, at the interface, the viscous part of the total normal stress and the shear stress are supported only by the fluid contained in the porous medium. An approach, which has been used in previous work dealing with the Brinkman-Forchheimer model, [13,14]. In what concerns heat transfer, thermodynamic equilibrium at the interface imposes that (T *\ =(t*) (10) V /fd v /pd and the energy balance at the interface gives (n*.q*'')fd=(n*.q*'')pd (11) where the total energy flux is defined as (q ) =u*T*___1 dT* \ * i'fd RePr dx*i (q '' * i) =u*iT*-------* \ i RePr dx*i (12) pd Numerical modelling The physical model is solved using a two-dimensional laminar version of the control-volume finite element method (CVFEM) described in [15]. Similar procedure is followed for the discretization of the Brinkman-Forchheimer equations, with the source terms described by the terms multiplying u*j in Eq. (4). The calculation domain, as fully reported in [10], is first discretized into three-node triangular elements. The nodes of this finite element mesh are the vertices of the triangular elements. The grid is designed so that there is a line of nodes along the entire interface between the Heat transfer distribution for a free/porous system with forced convection and heat generation 521 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 open and the fluid-saturated porous domains. The discretization of the calculation domain yields polygonal control volumes surrounding the nodes in the finite element mesh. Part of the resulting finite element mesh, and a polygonal control volume surrounding a node on the interface (P) are illustrated in Fig. 1. The thermophysical properties are stored at the centroid of each element and assumed to prevail over the respective element. In this work, these properties include mass density, dynamic viscosity, specific heat, and thermal conductivity of the fluid in the open domain, and Brinkman (or effective) dynamic viscosity, effective thermal conductivity, porosity, and permeability in the fluid-saturated porous medium. All dependent variables are stored at the same nodes in the finite element mesh, leading to a co-located formulation [15,16]. The momentum equations in the free and porous domains are integrated over each of the polygonal control volumes. Then, each of these integro-differential equations is approximated by a discretized equation that connects the dependent variables at each node in the finite element mesh to those at its immediate neighboring nodes. In the derivation of algebraic approximations to the various terms in the integro-differential equations, the dependent variables are interpolated over each element by similar interpolation functions, giving rise to an equal-order CVFEM. In the approximation of advection transport terms, the advected dependent variables are interpolated by flow-oriented exponential upwind functions in each element, and linear interpolations of the dependent variables are used to approximate the viscous or diffusion transport terms. The pressure is interpolated linearly in each element. The velocity components that are involved in the mass flux terms are interpolated using the so-called momentum interpolation scheme [17]. Full details of the formulation steps are available in [10]. The temperature field is calculated based on the discretized form of Eqs. 3 and 5, and using the converged velocity field. Figure 2. Control volume associated with a node P located on the interface between adjacent free (1) and fluid-saturated porous (2) regions of the calculation domain. 522 The overall iterative solution procedure can be summarized, at each iteration level, as follows: (i) based on guessed or latest available values of the velocity components, density, and viscosity, the discretized equations for the velocity components are constructed, except for contributions of the pressure-gradient terms, and the appropriate boundary conditions are introduced by manipulation of the coefficients and constants of these equations; (ii) pseudo-velocities and pressure coefficients are evaluated at each node; (iii) discretized equations for the pressure are obtained, and pressure or mass flow boundary conditions are introduced in the coefficients of such equations; (iv) a new pressure field is calculated; (v) the terms corresponding to the pressure gradient are introduced into the discretized momentum equations, based on the newly calculated pressure field, and the new velocity field is evaluated. This sequential procedure is repeated until suitable convergence criteria are satisfied. Then, the energy equations, Eqs. 3 and 5, in discretized form, and combined with the converged velocity field are solved iteratively. In this work, only structured grids with a line-by-line arrangement of the nodes were used. The discretized equations were solved using a simple line-Gauss Seidel iterative procedure, with a block correction algorithm, as a convergence accelerator. Detailed description of the overall iterative solution procedure can be found in [10,15]. This solution procedure was considered to have converged when the maximum normalized residues in the discretized momentum and pressure equations, and discretized energy equations were all less than 10 – 6. Model justification and testing In this section is reported a very small sample of the extensive testing to which the model was subjected, and which is reported in some detail in [7,10]. Transverse flow As already mentioned, models for edge flows based on the Darcy Law, which are formulated using an analytical formulation for the interface [1,2,3], or an equivalent permeability for the edge [1,2,4] do not take into account transverse flow. Close observation of the flow along the gap of width "d" (Fig. 1), reveals the strong influence of the transverse flow, which is depicted in Fig. 3, where the velocity vectors for the flow within the gap are shown. Figure 3 is nearly self-explanatory; Fig. 3(a) clearly denotes a flow that resembles plane Poiseuille flow, however, Figs. 3(b) and (c) present flow complexities, which are well beyond the capability of the Darcy Law – based models already discussed [1,2,3,4]. A.C.M. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 (c) Figure 3. Vectorial form of the velocity within the gap for the two-channel configuration (s = 0.7, Re = 10-1, Da 10- ) in three different regions, namely: a) Halfway the top channel; b) Bend between the top and RHS channel, and c) Outlet. Isothermal testing Flow through a plane channel with a porous plug is selected as an illustrative example. This problem involves flow through a parallel-plate channel with a porous plug, as shown schematically in Fig. 4, under an imposed overall pressure drop. For distances sufficiently far from the porous plug, the flow is fully developed. In the vicinity of the porous plug, the flow is not fully developed. Nevertheless, in this case, the fluid flows in a direction that is essentially normal to the interfaces between the open and fluid-saturated porous domains. As in the problem presented in the preceding subsection, the governing dimensionless parameters are the Reynolds number based on the mean velocity, Re = puH//j , the Darcy number for the porous domain, Da = K/H , the porosity of the porous domain, S , and the viscosity ratio, Vb/ M . Pin H Pout K- 3H H- 2H -h- 3H Figure 4. Geometry and computational domain for the plane channel with a porous plug. It is assumed that at distances of 3H upstream and downstream of the porous plug, the y-direction component of velocity, v, is zero and a cyclic boundary condition [10] can be applied to u, the streamwise component of the velocity. At the center of the porous plug, an essentially plug flow is expected if the permeability is low. The domain dimensions, a fixed value of the streamwise overall pressure drop (Ap =1.0 N/m2), the Darcy number, a fixed value of porosity (e = 0.7 ), and the viscosity ratio were imposed, and the properties p and pi were adjusted in order to obtain the desired Reynolds number, which is Re = 1 for all the results presented for this problem. Preliminary numerical tests of the asymptotic type have shown that a (21 + 21 + 21)x21 mesh, with a geometrically increasing node density towards the interfaces between the free and porous domains interfaces, in the streamwise direction, with a factor of 1.2, and a uniform node distribution in the y direction, gives essentially grid independent results. Results for Da = 10~ and 1 = 1 are presented in Fig. 5a where, on the left-hand side, is presented the centerline u velocity dependence on the streamwise dimensionless coordinate x/H , and on the right-hand side is presented the centerline pressure dependence on x/H . It is observed in Fig. 5a that the velocity field changes markedly due to the presence of the porous plug, this change being examined in this work via the variation of the u velocity along the centerline. On the left-hand and right-hand boundaries of the open domains, the u velocity variation in the y direction is essentially parabolic, and in the central region of the porous plug, this profile adjusts to be more flat (nearly a plug flow profile), with a lower centerline u velocity within the porous domain. The flow field is almost one-dimensional over most of the open and porous domains, but it is two-dimensional in the vicinity of the interfaces between the open and porous domains, where adjustments of the velocity field occur. In this case, as already mentioned, this adjustment occurs in a thin region, the extent of which scales as S/H ~ yjDa = 10~ . From the right-hand side of Fig. 5a it is observed that the centerline pressure presents distinct behaviors in the free and porous regions, with an essentially linearly decreasing profile in each, but with greater slope in the porous region, as expected. These results are in good qualitative agreement with the corresponding results in the literature. For Da = 10-3 and a =1, the results are presented in Fig. 5b. In this case, it is observed that marked changes in the centerline u velocity occur in a narrower region than that for Da = 10-2 and a =1, as now, s/h ~ -JDa = 10 3 . Also, in this problem, a lower Darcy number corresponds to lower permeability, thus the global pressure drop occurs almost fully in the porous domain, as illustrated by the right-hand side of Fig. 5b. Heat transfer distribution for a free/porous system with forced convection and heat generation 523 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 0.090 0.085 0.080 0.075 0.070 0.065 0.060 • ¦ ***N f 9 • \ / t L f v^, ~J 012345678 x/H (a) 0.040 0.036 0.032 0.028 0.024 0.020 "\ r m \ k~< ~J 0123 45678 x/H (b) 0.072 0.070 0.068 0.066 0.064 0.062 0.060 • i • • • ?—^ ^** * * • > • ' 3 V- 012345678 x/H 1.0 0.8 0.6 0.4 0.2 0.0 =\ 01 2345 678 x/H 1.0 0.8 0.6 0.4 0.2 0.0 012345678 x/H 1.0 0.8 0.6 0.4 0.2 0.0 X X (c) 0 1 23 4 5 67 8 x/H Figure 5. Centerline u velocity (left) and pressure (right) for Ap=1.0 N/m2, Ł=0.7 and Re=1 for (a) Da=10-2 and a=1; (b) Da=10-3 and a=1; and (c) Da=10-2 and a=5. Results for the same conditions as those considered previously (Fig. 5a), but with a = 5 , are presented in Fig. 5c. As the effective viscosity is now higher within the porous domain, the centerline u velocity presents a flatter profile than in Fig. 5a, and an increase in the pressure drop over the porous domain is observed. Results and discussion Numerical simulations were conducted using a nonuniform structured grid with 81 columns of nodes in the x direction and 55 rows of nodes in the y direction, with 17 rows of nodes in the channel region (of width d), and 13 rows of nodes in the inlet and outlet zones (of width D). The selection of this grid is determined on the basis of extensive grid convergence testing [18]. The results were obtained for the geometry presented in Fig. 1, considering l/H = 5, D/H =0.1, Ł = 0.54, a = 1, Rc =100, and q, g , ,*=109. The remaining governing dimensionless parameters were subjected to change in order to observe their influence on the resulting flow and temperature distributions. The first analysis that should be undertaken is related with the gap width. From Figs. 6 and 7, it can clearly be 524 Sousa A.C.M. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 observed that a change from d/'H = 0.04 to d/H =0.02, (i.e. a decrease of 0.02 in the dimensionless mold cavity clearance), results in substantial alterations in what concerns the flow through the porous medium. Figure 6 clearly shows a preferential flow through the top and right hand side channels in detriment of the flow through the porous medium, whereas in Fig. 7 it can be observed that a decrease in the clearance width minimizes the edge flow, yielding an increase of the flow through the porous medium. A pressure fall can be observed near the inlet region for d/H = 0.04 whereas, for djH =0.02, the pressure is more uniformly distributed along the bed. In what concerns the global dimensionless pressure difference, it can be observed that this parameter is considerably higher for the smaller gap (i.e. for d/H =0.02 and (Ap*=3.1x105) for d/H = 0.04. In terms of temperature fields, the larger gap leads to a larger region subjected to high temperatures in the vicinity of the mold exit. Furthermore, the maximum global dimensionless temperature difference is considerably higher for the larger gap (i.e. AT* =3.5x104for d/H = 0.04, and AT* =1.9x10 for d/H =0.02). The influence of Darcy number is also relevant, as expected. When comparing the results shown in Figs. 6 and 8 for Da = 10 and Da = 10-8, respectively, representing the permeability variation of the porous medium, it can be noted the emphasis on the preferential flow through the gap towards the exit (i.e. greater contribution to the edge flow effect). In this case, the main pressure decrease occurs closer to the inlet region, with the global dimensionless pressure difference assuming a value of Ap* =2.4x107 . As the interior of the porous medium remains essentially without flow, the hot region increases both in area and temperature level, where the global dimensionless temperature difference becomes AT =4.8x10 . Figure 6. Dimensionless streamlines, isobars and isotherms for d/H = 0.04, Da = 10~6 , Re = 1 and Pr = 1000. The influence of the Reynolds number was studied for Re=1 and Re=0.1 with d/H=0.04, Da=10-6, and Pr=1000. Both figures present similar behavior in what concerns the flow field (isobars and streamlines), however the global dimensionless pressure difference is considerably higher for Re=0.1 than for Re=1, and is equal to Ap* =3.1x106. The isotherms do not exhibit a hot spot and the temperature rise is quite uniform along the bed. However, a lower Reynolds number leads to higher temperatures, with the maximum global dimensionless temperature difference reaching a value of AT =2.2x10 . ?\\ - ¦^ ' / 7 _ Figure 7. Dimensionless streamlines, isobars and isotherms for d/H = 0.02, Da = 10~6, Re = 1 and Pr = 1000. Figure 8. Dimensionless streamlines, isobars and isotherms for d/H=0.04, Da = 10~8, Re = 1 and Pr = 1000. Figure 9. Dimensionless isotherms for d/H =0.04, Da = 10~6, Re = 0.1 and Pr = 1000. In what concerns Prandtl number, it should be noted from Fig. 10, that a rise of this parameter to 10000 from 1000 (Fig. 6) infers a significant fall on the temperature level along the interface i.e AT* = 565 vs. AT* =3.5x10 . The combination of higher Prandtl number with lower Reynolds number (Fig. 11), where the fluid assumes, as compared to the situation depicted in Fig. 6, a greater heat removal capability and lower mass flow rate, respectively, becomes evident from the isotherms, where AT=3.5x104. Heat transfer distribution for a free/porous system with forced convection and heat generation 525 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 519-526 Figure 10. Dimensionless isotherms for d/H =0.04, Da = 10~6, Re = 1 and Pr = 10000 . Figure 11. Dimensionless isotherms for d/H = 0.04, Da = 10 6, Re = 0.1 and Pr = 10000. Conclusion The present work reviews a detailed physical model for situations involving fluid flow and heat transfer in fluid domains partially filled with a porous medium. Special attention was given to the boundary conditions at the interface between the open and porous domains, and to the situations of steady volumetric heat generation. In addition, the physical justification for the model is reported along with an illustrative testing example. The results were obtained for a set of governing dimensionless parameters relevant to the LCM processes. The analysis of the results highlights the influence of the gap width over the flow and thermal fields. In what concerns the thermal field, the temperature distribution is highly affected by the range of the governing parameters. It was found that large regions in the mold are susceptible to become subjected to deficient wetting and poor heat removal capability due to building up of high temperature levels within the component. The predictions identify conditions that may lead to the production of components of inferior quality, which, in extreme cases, have to be rejected. Acknowledgement The author acknowledges the finanical support received for this work from FCT (Portugal) - Project POCTI/EME/36263/2000, and NSERC (Canada) -Discovery Grant 12875 (ACMS). The valuable and numerous contributions to this research project made by Profs. V.A.F. Costa and M.S.A. Oliveira (U. Aveiro), Prof. LA. Oliveira (U. Coimbra), and Prof. B.R. Baliga (McGill U.) are also acknowledged References [1] [2] [3] S. Bickerton, S.G. Advani, Characterization and modeling of race-tracking in liquid composite molding processes, Composites Science and Technology 59 (1999) 2215-2229. A. Hammami, R. Gauvin, F. Trochu, Modeling the edge effect in liquid composites molding, Composites: Part A, 29A (1998) 603-609. G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech., 30 (1967) 197-207. [4] W.-B. Young, C.-L. Lai, Analysis of the edge effect in resin transfer molding, Composites, 28A (1997) 817-822. [5] K.M. Pillai, Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows, Composites, 33A (2002) 1007-1019. [6] S. Bickerton, E.M. Sozer, P. Simácek, S.G. Advani, Fabric structure and mold curvature effects on preform permeability and mold filling in the RTM process. Part II. Predictions and comparisons with experiments, Composites: 31A (2000) 439-458. [7] V.A.F. Costa, M.S.A. Oliveira, A.C.M. Sousa, Numerical simulation of non-Darcian flows through spaces partially filled with a porous medium, in: A.A. Mohamad (Ed.), Proc., The 3rd Intal. Conference on Computational Heat and Mass Transfer, U. of Calgary Press, Calgary, Canada, 2003, pp. 135-143. [8] V.A.F. Costa, A.C.M. Sousa, Modeling of Flow and Thermo-Kinetics During the Cure of Thick Laminated Composites, International Journal of Thermal Sciences, 42 (1) (2003) 15-22. [9] M.K. Alkam, M.A. Al-Nimr, Transient Non-Darcian Forced Convection Flow in a Pipe Partially Filled with a Porous Material, Int. J. Heat Mass Transfer, 41 (1998) 347-356. [10] V.A.F. Costa, L.A. Oliveira, B.R. Baliga, A.C.M. Sousa, Numerical simulation of coupled viscous flows in adjacent porous and open domains using a control volume finite element method (accepted for publication in Numerical Heat Transfer, Nov., 2003). [11] M. Kaviany, Principles of Heat Transfer in Porous Media, Springer Verlag, New York, 1991. [12] D. A. Nield, A. Bejan, Convection in Porous Media, 2nd. Ed., Springer Verlag, New York, 1999. [13] D. K. Gartling, C. E. Hickox, R. C. Givler, Simulation of coupled viscous and porous flow problems, Comp. Fluid Dynamics, 7, (1996) 23-48. [14] A. G. Salinger, R. Aris, J. J. Derby, Finite element formulations for large-scale, coupled flows in adjacent porous and open fluid domains, Int. J. for Num. Methods in Fluids, 18, (1994) 1185-1209. [15] V.A.F. Costa, L.A. Oliveira, A.R. Figueiredo, A Control Volume Based Finite Element Method for Three-Dimensional Incompressible Turbulent Fluid Flow, Heat Transfer, and Related Phenomena, International Journal for Numerical Methods in Fluids, 21 (1995) 591-615. [16] B.R. Baliga, Control-volume finite element methods for fluid flow and heat transfer, in: W.J. Minkowycz, E.M. Sparrow (Eds.), Advances in Numerical Heat Transfer, vol. 1, Chapter 3, 1997, pp. 97-135. [17] C. Prakash, S.V. Patankar, A control-volume-based finite-element method for solving the Navier-Stokes equations using equal-order velocity-pressure interpolation, Num. Heat Transfer, 8 (1985) 259-280. [18] V.A.F. Costa, M.S.A. Oliveira, A.C.M. Sousa, Numerical evaluation of heat dissipation by non-Darcian forced convection in porous/fluid systems, Heat Transfer VIII, Adv. Computational Methods, B. Sundén et al. (Eds.), Southampton, U.K., 2004, 99-108. (Invited Contribution) 526 Sousa A.C.M.