Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 SOLVING TURBULENT FLOW DYNAMICS OF COMPLEX, MULTIPLE BRANCH RADON MITIGATION SYSTEMS L. Moorman, Ph.D. Radon Home Measurement and Mitigation, Inc. Fort Collins, CO ABSTRACT Active radon mitigation systems with a large flexibility in their complexity are simulated via computer calculations assuming incompressible gas flow dynamics. When designing more complex, multiple branch systems, insight into the turbulent flow dynamics inside the systems vent pipes may be of advantage in design and optimization of pipe sizing and radon ventilator characteristics. I will show how their effects can be simulated and assessed during the installation of real systems. For an N-branch, active, non-looped, radon mitigation system I will show how to construct the complete set of non-linear equations with its independent variables and demonstrate how these can be uniquely solved and applied to optimize a system under design. A number of real systems were used to compare with calculations. I will show a number of applications using this program: 1) The effects of sub-slab material and its resistance on the system. 2) The effects of turbulent flow through vent pipes of various diameters. 3) The effects of varying ventilator characteristics. 4) Calculations of multiple pipe branches in real time. INTRODUCTION The equation that governs the pressure distribution for the air in a single branch radon system in its simplest form consists of three components, the pressure loss through the soil and the radon extraction cavity, the pressure loss due to the resistance of the air with the pipe system and the pressure boost by the action of an operating ventilator somewhere in the pipe system. I will first discuss each of these components separately and indicate in each case which theoretical functional behavior approximates the description of the various pressure differences best. After that I will put together the components and describe how the equations for complete radon mitigation systems can be derived for multiple branches. Numerical evaluation of the solutions will result in gaining quantitative insight into the use of the equivalency of pipe sizes in complex radon systems. SUB SLAB MATERIAL AND CAVITY RESISTANCE PARAMETERIZATION. A dimensionless parameter that characterizes the flow for a gas through an opening or capillary is the Knudsen number. This is defined as the ratio of the mean free path of particles in the gas divided by a characteristic dimension such as the diameter of the pipe: ! Kn = (1) d In the movement of atoms or molecules of gas through small capillaries of the sub-slab medium, such as soil, rock or concrete, the flow may be more determined by the 1 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 collisions with the walls of the pores in the porous material than by collisions among the gas particles. In that case, the viscosity of the gas (which is a parameter for the friction inside the gas) does not play a large roll in the dynamics of the system. Therefore, on the one hand, very slow air movement inside a medium when no active soil depressurization exists, direct diffusion of radon gas through small capillaries may be the largest mode of transportation for radon. On the other hand, when enough small cracks and openings exist, or if the pore size in the material is large enough, convective movement of air which drags radon with it through these openings is more likely to dominate the mode of transportation for radon. The sub-slab material plus extraction cavity resistance (friction) for air movement when an active soil depressurization system operates can be measured in situ by placing a fan with relevant power and flow characteristics directly over the final excavated extraction cavity and measuring the pressure across the fan while running the fan at various speeds. When varying pressure differences, ΔPf, and measuring corresponding flow rate s, F, the curve that is generated appears in practice sufficiently linear in the relevant flow range for active mitigation systems. Fig. 1 shows an example of these parameters that were measured in existing homes during mitigations. Fig 1: Flow measured as a function of pressure difference across ventilator in the configuration described in text for various residences and different branches in one residence. 2 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Flows were measured by timing the capture of air from the end of the fan into a well determined size bag of known volume. Accuracy was verified by employing the same method to measure the fan curves as advertised by manufacturers for various fans. Flow measurements were found to be somewhat lower than the manufacturer’s specification, but sufficiently accurate to be able to use this method for this purpose. The resulting slope in Fig 1 is independent of which ventilator is used on a radon extraction cavity for the measurement. The observed linearity allows us to define a (local) linear friction Rg for the combined soil with extraction cavity through the following relationship: !Pf = PFR g (2) This equation can also be understood by interpreting the right hand side as the combined soil pressure drop inside the cavity (ΔPs=PFRg) which must counteract the pressure boost (ΔPf) by the fan in this simple model. The experimental linearity may be related to the fact that air velocities further away from the extraction cavity, representing the majority of the sub-slab material, are small. Deviations from linearity are expected at higher flow rates. By dividing both sides of the Eq. (2) by pressure and flow rate I obtain an equation for the friction: !Pf Rg = (3) PF The unit of resistance is independent of the units used for pressure. If flow is expressed in cubic feet per minute (cfm) this resistance expresses how many minutes it takes for one cubic foot air to move through the extraction cavity (min/cf) [in SI units the flow would be expressed in m3/s and the resistance would be in s/m3]. In Fig. 2 a graph is shown of a number of measured soil resistances graphed against the relative fan depressurization, which is the measured depressurization across the fan compared to the maximum fan depressurization at zero air flow. What can be seen in this “landscape of resistances” is a correlation of low soil resistance with low relative fan depressurization, as in large gravel, and a correlation of high soil resistance with high relative fan depressurization, as is the case in very tight clay. Other sub-slab material types found in between are designated as pea gravel, sand, looser clay, hard sand and clay. It must be emphasized here that the measurement of the resistance is unique to the specific residence and includes not only the sub-slab material type 3 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Fig 2: Resistances measure of excavated radon cavities in existing homes following method of Fig 1. but also any inhomogeneities in the sub-slab material and any leaks due to openings in the slab. The measured resistance defined here is a good quantification of the effective resistance that a radon mitigation system will be subject to via this radon extraction cavity. In this sense, I have found a quantification for the resistance of the extraction cavity with sub-slab material which is also valid when dealing with a multi-branch radon mitigation system. Thus this information can be used in simulations that include the complete system with all branches as will be shown later. VENT PIPE MODELING Laminar, viscous flow is indicated by values of the Knudsen parameter, defined in Eq. (1), smaller than 0.01 indicating that collisions between the particles of the gas are much more abundant and are determining the character of the flow rather than particle-to-wall collisions. This is the situation in air moving through a pipe at small speed and flow rates. At small flow rates, for viscous, laminar flow through a cylindrical tube the velocity of a fluid with viscosity η can be shown to be a parabolic (quadratic) function of the distance from the axis with the largest velocity on the central axis of the tube, and a zero velocity near the wall. The volume rate of fluid crossing any section of the tube with radius r for this flow is given by: #r 4 F= !P (4) 8"L 4 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 This is known as Poiseuille’s formula and is important in medicine since it gives a qualitative understanding of blood flow through our arteries and veins. In addition, by reading Eq. (4) from right to left it can be interpreted that under steady state flow conditions the pressure loss along any section of a cylindrical pipe of length L is a linear function of the flow rate through a cross section of the pipe: 8"L #P = 4 F (5) !r However, when the flow velocity of a fluid becomes sufficiently large, laminar flow breaks down and turbulence sets in. The critical point for the onset of turbulence in a fluid can be characterized by a dimensionless parameter referred to as the Reynolds’ number, which is the ratio of the shear stress in the fluid due to turbulence and the shear stress due to viscosity given by: 2r"v NR = (6) ! In this equation ρ is the density of the fluid and v the velocity. Experiments have shown that the flow will be laminar when the Reynolds’ number is below 2000 and the flow will be turbulent for values above 3000. In between these values the flow can be either laminar or turbulent and can go back and forth. When the air speed increases through a pipe, the same model applies and when the Reynolds’ number increases beyond its critical value the air flow becomes turbulent and the resistance increases. For such situations the linear relationship between pressure loss along the pipe and flow rate is no longer valid and the new relationship is a nonlinear function of flow rate that must go through zero. Thus I can extend the pressure loss behavior to a more general form with two parameters which may be different for each diameter of the pipe: !P( F ) = cF a (7) The variables for the pressure drop along 100 ft of the pipe for different diameters are known from the heating and ventilation industry and values for 2”, 3”, 4” and 6”round tubes are known. The derivative of this pressure loss to flow can be written as: "#P ( F ) = caF a !1 "F The length must be scaled appropriately to the total equivalent length of the pipe where 90 degree elbows and 45 degrees elbows have a certain equivalent length in air resistance to straight section of pipe of the same diameter. This equivalent length approximation is used in the heating and ventilation industry for ducts and assumes that the nonlinearity (power of F) of the components is the same for connections as for straight pipe. If this is not the case a different equivalent length would have to be determined for each component at each air flow rate. In the following I will use the assumption that the equivalent length calculation can be made. VENTILATOR CURVE MODELING 5 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 The modeling of the ventilator curve (fan curve) can be done sufficiently accurate for our goal by three parameters. I assign a new function with three parameters to each ventilator and have verified it can be accurately described for each fan (f) by: (8) g ( f , F ) = q1 + q 2 F b2 This means the derivative to the flow rate of the fan curve is: "g ( f , F ) = q 2 b2 F b2 !1 "F In our model of multiple branch radon mitigation systems I will allow for the possibility that any pipe section between nodes (pipe TEE’s) can have a fan inserted. SINGLE BRANCH RADON MITIGATION SYSTEM For the calculation of a single Branch radon mitigation system, I can now add the individual components keeping in mind that all pressure losses due to sub-slab material and cavity resistance and resistance in pipe system together are counteracted by the pressure boost from the ventilator, resulting in the nonlinear one-dimensional equation: a PF1 R g + c1 F1 ! g ( f , F1 ) = 0 (9) Once the resistance for the sub-slab material (Rg) is measured, it can be shown that a unique solution exists for each fan and pipe diameter. By dividing the equation by pressure and flow I find: a !1 F1 g ( f , F1 ) (10) ! =0 P PF1 This can be employed in a graphic solution even during the mitigation as is shown in the fan resistance response curve in Fig. 3. R g + c1 For any given ventilator, relevant calculations can be made that will give us insight in the flow dynamics of the system by numerically evaluating the solution to Eq. (9) as a function of equivalent length for various cavity resistances. The resulting curves for one commonly used fan are shown in Fig 4. It shows that the relative change of volume rate through the system is very different depending on the resistance. For low cavity resistance the flow rate depends strongly on the equivalent length for short equivalent lengths only whereas for high resistances the flow rate is nearly independent of the equivalent length. 6 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Fig 3: Fan Resistance Response curve: Method of graphically solving Eq. (10). For 4” diameter pipe system and 175 ft equivalent length the red squares indicate the resistance for all flow rates. To solve the equation find the crossing with the fan curve (green and blue squares) indicated by a red circle and shift the circle up by the extraction cavity friction (heavy red solid line) along the fan curve to find the value inside the second circle. The arrow points at a flow rate of approximately 80 cfm which is the steady state flow that this ventilator in this mitigation system will support. For a 3” diameter pipe it is indicated at the second arrow that the flow will go down to approximately 50 cfm., thus loosing 38% of the flow rate through this system. The dashed horizontal line indicates the resistance of a different extraction cavity that was measured with a higher resistance which leads to a flow rate solution even with a 4” pipe of approximately 50 cfm. Similarly I have evaluated this information given in figure 4 for different diameter pipes, which is most relevant for 3” and 2” pipes. I can than evaluate the ratio for flow rate of the 3” pipe with 4” pipe. This ratio is given in Fig. 5 as a function of equivalent length for each cavity resistance. It can be seen that suppression effect of flow for smaller pipe sizes are strongest for low cavity resistances (crawlspace, gravel) and weakest for large cavity resistances (hard sand, or clay). This can be qualitatively reformulated by the statement that high flow rates are easier choked by a smaller diameter pipe than low flow rates, and that when the flow rate is small the effect of pipe diameter is not very important. 7 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Fig4: For a FR150 fan with 4-inch pipe for a single branch system the volume rate is numerically calculated for several extraction cavity resistances as function of equivalent length. The information obtained from these calculations and for similar calculations for 2” diameter vent pipes can be summarized by defining a criterion for flow rate loss. I choose to extract from Fig. 5 a critical equivalent length that corresponds to a flow rate loss of 20% for each resistance. These critical equivalent lengths can be graphed against the sub-slab material and cavity resistance as has been done in figure 6 which in turn allows for the following interpretation. The regions divided by the solid and dashed lines are the regions in which there is a certain equivalence of system. It can be seen that in some cases less than 20% loss is expected by switching to a smaller pipe diameter. This is indicated by the notation: 3” or 4” inside the graph. Similarly, there is a region at large cavity resistances and shorter length systems where the flow dynamics for even 2” pipe would be similar to 3” and 4” pipe diameters. On the other hand it is clear that at small cavity resistances only for very short systems, up to 10 ft equivalent length, the 3” and 4” systems are similar and at 25 ft equivalent lengths one can expect already as much as 30% flow loss when switching to a 3” vent pipe. Nevertheless, the correct interpretation of this graph is that it shows equivalence of systems but does not show whether either of the equivalent systems will work to reduce radon levels. 8 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Fig. 5: For a FR 150 fan for a single branch system the ratio between the volume rate employing 3 inch and 4 inch piping in the system is shown. Fig 6: Equivalent Dynamics Graph (Fan specific, here for FR150). Flow reduction data points for 20% flow loss are connected with straight lines in the graph. The areas in between the lines are marked for which pipe diameters are equivalent within the chosen criterion of the flow loss. The open symbols with dashed lines are for 30% loss. TWO BRANCH RADON MITIGATION SYSTEM 9 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 In deriving the equations for a 2-Branch system I can develop a set of three equations that are similar to the Kirchhoff equations in electronic circuits with the one important difference that the equations here are nonlinear in the flow and that it is not immediately clear that only a single flow solution exists. i) Fan equation: !P2 + !P1 = !Pf ii) iii) Kirchhoff Rule through the two branches: !P2 = !P3 Flow conservation: F1 = F2 + F3 (11) By substituting iii) in i) the set of equations can be written as two nonlinear equations in two independent variables (F2 and F3 ) only: PF2 R g 2 + c 2 F2a 2 + c1 ( F2 + F3 ) a1 ! g ( f , F2 + F3 ) = 0 i’) (12) ii’) PF2 R g 2 + c 2 F2a 2 ! PF3 R g 3 ! c3 F3a3 = 0 The equations can be written in a vector formulation as r r r (13) A( F ) = 0 and solved numerically by finding the root of the two dimensional equation. The solution of Eq. (13) can be found numerically by use of an assumed reasonable r starting value of the flow vector, F0 , and, by iteration, finding the small root deviation vector from the previous approximate solution. For iteration number n I find: "1 ( % r r . /Ai + & # A ( F ) = " L"1 A ( F ) )) d n = "! ,, (14) ! q n pq q n &- /Fk * r # q q &' Ft # $ pq In which I defined the Jacobian matrix: rr & 'A # (15) L = $$ i !! % 'Fk " and the inverse of this matrix is: 1 & L22 ' L12 # $ ! (16) L'1 = det L $% ' L21 L11 !" I than find the next iteration of the solution vector in the two dimensional “flow space” from the previous approximate solution by r addition: r r (17) Fn +1 = Fn + d n r I then choose to continue the iteration until a certain stopping criterion d n " ! , for a sufficiently small value of ε, is reached, at which time the closest approximation to the r solution is Fn for final value of t reached. The method described here to solve Eq. (13) is mathematically known as the Newton–Raphson method of root finding for a nonlinear system of functions with multiple independent variables. N-BRANCH RADON MITIGATION SYSTEM 10 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 More general solutions have been derived for 3- and 4-Branch radon mitigation systems with a possibility of up to 7 ventilators in the various branches of the mitigation system. From this, it is possible to formulate an N-Branch mitigation system with (2N-1) ventilator insertion points. In a non-looped system as is shown schematically in Fig. 7 this would be one ventilator for each of the N branches, one ventilator in the common upper stack of the system (Branch 1) and (N-2) ventilators for the interconnecting pieces between the nodes of the N pipe branches. As an example, when deriving the equation for 4 branches I write down the fan equation, three Kirchhoff equations and the continuity equation. This leads to four non-linear equations with four independent unknown variables. The four independent variables are the flow values in each of the four branches that connect to the ground: & F2 # $ ! r $ F3 ! (18) F =$ ! F4 $ ! $F ! % 5" Fig. 7: Schematics of a 4-branch radon mitigation system with the parameter definitions used in the derivations. and the 4-dimensional mitigation system equation thus is: r r r (19) A( F ) = 0 11 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Similar to Eq. (15) I can derive the 4-dimensional Jacobian matrix: In this: &/ + . rr & *A # $ . $ L = $$ i !! = $ % *Fk " $ 0 $ 0 % / / / # ! '- ' ) ') ' )! ', ' ( ' (! ! 0 , ' + !" (20) ( = a ! f1 ' = b ! f2 & = c ! f3 % = d ! f4 $ = e ! f5 # = A ! FA " = B ! FB in which the flow resistance components of the upper stack are:: a = PF1(2345) = c1 a1 ( F2 + F3 + F4 + F5 ) a1 !1 A = PA(3B) = c A a A ( F3 + F4 + F5 ) a A !1 B = PB(45) = c B a B ( F4 + F5 ) aB !1 and in which the flow resistance components of the branches with sub-slab material connections are: b = RP 2 = PRg 2 + c2 a2 F2a2 !1 c = RP3 = PRg 3 + c3 a3 F3a3 !1 d = RP 4 = PRg 4 + c4 a4 F4a4 !1 e = RP5 = PRg 5 + c5 a 5 F5a5 !1 In addition the active ventilator boost derivatives in the common upper part of the vent stack is given by: f1 = f1 (2345) = q1b1 ( F2 + F3 + F4 + F5 ) b1 !1 and for the individual branches i = 2, 3, 4, or 5 these ventilator boost derivatives are given by: f i = f i (i ) = qi bi Fi bi !1 For the interconnecting sections A and B between the various node points the ventilator boost derivatives are given by: f A = q A b A ( F3 + F4 + F5 ) bA !1 f B = q B bB ( F4 + F5 ) bB !1 Solving the set of equations defined in Eq. (19) is similar to the approach of solving Eq. (13), as was done in Eq. (17) with use of Eq. (14), following the Newton-Raphson 12 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 method for a set of non-linear equations with multiple independent variables. In this case the 4x4 Jacobian matrix In Eq. (20) must be inverted. Although there is an analytical expression for this inversion, the numerical evaluation of this is cumbersome. Instead, a numerical evaluation of the inverted 4x4-matrices was applied directly via the Quattro Pro matrix inversion function and solutions of the set of equations were calculated with an iterative approach by using three embedded macros. Verification of the validity of the inversion process was done by multiplication with the original matrix and by testing the convergence of the solution. With fairly generic starting values of the flows, the four dimensional calculations converged most of the time in less than 8 steps to a result of the final flow that is close to the accuracy of the computer. The equations presented contain all of the physics that was discussed before. Additional physical effects that can be included are gravitational effects which can be responsible for a fan to act differently when moving air horizontally over 25 ft compared to moving air over a vertical lift of 25 ft, as well as the lift due to temperature differences of air inside the pipe compared to outside air (the stack effect for a passive system). Also certain aspects of the Bernoulli effect can be included in these equations, and have not yet been included in the above formulation. Finally I must mention that these equations are for incompressible flow, which in its defense is a sufficiently accurate approximation for the air parameter regime in which radon systems generally operate. In Fig 8 I show an example of a calculation for a house where we installed a three branch system with the HP220 ventilator in the attic of the garage 2 ft under the discharge point. Slab areas were 600, 1050, and 100 sf to be serviced by the three branches 2, 3 and 4, respectively. The lower half of this figure shows the vacuum (negative) pressures that the computer program calculated to solve Eq. (19) for the cavity parameters measured and equivalent lengths of pipes calculated. The horizontal axis gives relevant observation points along the system with the Sky-limit pressure on the right and the Sub-slab limit pressure on the left, which will always be the same and at zero pressure. At the mark “P-hole” the pressures in each of the three extraction cavities is given and “P1-below fan” is the pressure in the top stack of the system immediately below the ventilator in the attic of the garage. The program allows for placing the system at further distances away from the discharge point for instance when installing an outside system.. At the marker P_45B the pressure is given at the node where branches 4 and 5 come together and the air flows into pipe B. The same is the case for the other node points. The top half of the figure shows the flow rates indicated on the right vertical axis of the numerical solution that the computer calculated for each marker along the horizontal axis. In this case an HP220 fan was simulated, the top stack was 45 equivalent feet long and branch 2, 3, and 4 each 65 equivalent feet long with 2, 3, and 4 inch diameter pipe, respectively and branch 5 was 40 equivalent feet long with 3 inch diameter pipe. A and B were 10 and 20 equivalent feet of 4 inch diameter pipe. 13 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Cavity resistances were for cavity 2: 2×10-5 min/cf, for cavity 3: 8×10-5 min/cf and for cavity 4: 1×10-3 min/cf and cavity 5 was essentially deactivated in this calculation by giving it a virtual resistance of 100 min/cf. It is interesting to see that the flows indicated at P_hole through branch 2, even though the vacuum pressure is small, is almost as large as through branch 3, despite the fact that branch 2 only uses 2 inch piping and branch 3 has 3 inch piping. This is due to the fact that the cavity resistance of branch 2 was a factor four lower than the cavity resistance for branch 3. Fig. 8: The numerical program shows the resolution in a graph that can be conveniently read of. The vacuum pressures in the three radon extraction cavities under the slab are given at P_hole, as are the three airflows in the individual branches. When calculating Reynolds numbers the program showed that branch 5 and B were operating in the laminar regime, rather than the turbulent regime as is the case for the other branches as shown in Fig. 9. 14 Proceedings of the American Association of Radon Scientists and Technologists 2008 International Symposium Las Vegas NV, September 14-17, 2008. AARST © 2008 Fig. 9: Reynolds numbers in the various pipe branches indicating the flow dynamics applicable. CONCLUSION A theory for single-branch radon mitigation systems was developed in which experimental data were used to support the theory that the sub-slab material and extraction cavity resistance can be sufficiently well described by a single resistance parameter. The differences between laminar and turbulent flow through round pipes were discussed and the roll of the Reynolds’ number as a criterion for the onset of turbulence in the air flow. A method utilizing the Fan Resistance Response graph was developed to aid in graphically solving the 1-branch mitigation system equation. By numerically solving the single branch radon mitigation equation for many different situations I have derived criteria for the equivalence of different size piping under various circumstances, and developed a quantitative insight in those circumstances when pipe sizes cannot be considered equivalent. These insights were graphically summarized in the Equivalent Dynamics Graph which can be helpful during installation of radon mitigation systems possibly leading to more appropriate pipe size choices in the future. As an example of an N-Branch system, for a four-branch, active, non-looped, radon mitigation system I derived the complete set of non-linear equations with its independent variables and demonstrated how these can be solved numerically and applied to optimize a system under design. A three-branch example was discussed in detail showing that this system has simultaneously turbulent and laminar flow in different branches. 15