Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 86
Текст из файла (страница 86)
Start with a guessed pressure and velocity fields p(n) and u(n), respectively.2. Solve the momentum equation given by Eq. (15.27) to obtain a new velocityfield uf .3. Update the mass flow rates using the momentum satisfying velocity field toobtain the m_ f field.4. Using the new mass flow rates solve the pressure correction equation to obtain apressure correction field p0 .5.
Update the pressure and velocity fields to obtain continuity-satisfying fieldsusing the following equations:15.2A Preliminary Derivationuf¼573ufþu0fu0f¼Duf 0@p@x fðnÞpC ¼ pC þ p0C_ f þ m_ 0fm_ f ¼m 0@p_m0f ¼ qf Duf Dyf@x fð15:45Þ6. set uðnÞ ¼ u and pðnÞ ¼ p7. Go back to step 2 and repeat until convergence.The SIMPLE algorithm is best illustrated via the example presented next.Example 1Flow in a Pipe NetworkA portion of a water pipe system is shown in Fig.
15.4. The momentumequation for the flow in the pipes can be written asm_ ¼ quA ¼ DDPwhere DA = 0.5, DB = DF = 0.4, DC = DE = 0.3, DD = 0.19, DG = 0.1875,and DH = 0.35. Using the SIMPLE algorithm, calculate the unknown massflow rates and pressures in the system.p2 = 350p5 = 300mBmEmDmAp1 = 400p6mCp4 = 50mI = 50mGp3mFp7 = 80p8mHp9 = 200Fig. 15.4 A portion of a water pipe systemSolutionIn this system, the mass flow rate field is used as a variable instead of thevelocity field. This is not problematic since the momentum equation has beensimplified by dropping the convection and diffusion terms as their values arenegligible compared to the pressure head.57415 Fluid Flow Computation: Incompressible FlowsThe solution using the SIMPLE algorithm starts by first computing aninitial velocity field using the assumed pressure field, and then predicting apressure field that enforces continuity on the just computed velocity field.This procedure is summarized as1.
Start with a guessed pressure field.2. Compute the mass flow rates using the given momentum equation.3. Construct a pressure correction equation that enforces continuity (massconservation) and use it to correct the pressure and velocity fields.No iterations will be needed since there are no non-linear effects inducedby a convection term.step 1Start by assigning guessed values to the pressure at the locations wheresolutions are to be found. Thus assume the following:ðnÞðnÞðnÞp3 ¼ 300 p6 ¼ 200 p8 ¼ 120 (other values could have been used)step 2Based on the assumed pressure values calculate the various mass flow ratesusing the momentum equations according toðnÞm_ A ¼ DA ðp1 p3 Þ ¼ 0:5 ð400 300Þ ¼ 50ðnÞm_ B ¼ DB ðp3 p2 Þ ¼ 0:4 ð300 350Þ ¼ 20ðnÞm_ C ¼ DC ðp4 p3 Þ ¼ 0:3 ð50 300Þ ¼ 75ðnÞðnÞm_ D ¼ DD ðp3 p6 Þ ¼ 0:19 ð300 200Þ ¼ 19ðnÞm_ E ¼ DE ðp6 p5 Þ ¼ 0:3 ð200 300Þ ¼ 30ðnÞm_ F ¼ DF ðp7 p6 Þ ¼ 0:4 ð80 200Þ ¼ 48ðnÞðnÞm_ G ¼ DG ðp6 p8 Þ ¼ 0:1875 ð200 120Þ ¼ 15ðnÞm_ H ¼ DH ðp9 p8 Þ ¼ 0:35 ð200 120Þ ¼ 28step 3P Check whether the mass flow rates satisfy continuity by computingm_ k atkall interior points, i.e.,X Node 3 :m_ k ¼ m_ B þ m_ D m_ A m_ C ¼ 20 þ 19 50 þ 75 ¼ 24kX m_ k ¼ m_ G þ m_ E m_ D m_ F ¼ 15 30 19 þ 48 ¼ 14Node 6 :kNode 8 :Xkm_ k ¼ m_ I m_ G m_ H ¼ 50 15 28 ¼ 715.2A Preliminary Derivation575Since mass conservation is not satisfied, correction fields are needed andpressure correction equations are derived as follows:XkX X m_ k þ m_ 0k ¼ 0 )m_ 0k ¼ m_ kkkIn term of pressure corrections, mass flow rate corrections can beexpressed asm_ 0A ¼ DA p03m_ 0B ¼ DB p03m_ 0C ¼ DC p03m_ 0D ¼ DD p03 p06 m_ 0E ¼ DE p06m_ 0F ¼ DF p06m_ 0G ¼ DG p06 p08m_ 0H ¼ DH p08Note that p01 , p02 , p04 , p05 and p07 are set to zero since the correspondingpressure values are known and hence represent the exact values.The flow field at nodes 3, 6, and 8 are respectively given bym_ 0B þ m_ 0D m_ 0A m_ 0C ¼ m_ B þ m_ D m_ A m_ Cm_ 0G þ m_ 0E m_ 0D m_ 0F ¼ m_ G þ m_ E m_ D m_ F m_ 0G m_ 0H ¼ m_ I m_ G m_ Hand using pressure corrections as DB p03 þ DD p03 p06 DA p03 DC p03 ¼ m_ B þ m_ D m_ A m_ C DG p06 p08 þ DE p06 DD p03 p06 DF p06 ¼ m_ G þ m_ E m_ D m_ F 0 DG p6 p08 DH p08 ¼ m_ I m_ G m_ HAfter simplification the above equations become 1:39 p03 0:19 p06 ¼ m_ B þ m_ D m_ A m_ C 1:0775 p06 0:1875 p08 0:19 p03 ¼ m_ G þ m_ E m_ D m_ F 0:1875 p06 þ 0:5375 p08 ¼ m_ I m_ G m_ H57615 Fluid Flow Computation: Incompressible FlowsSubstituting the tentative mass flow rates, the various correction fields satisfy 1:39 p03 0:19 p06 ¼ 24 0 1:0775 p6 0:1875 p08 0:19 p03 ¼ 14 0:1875 p06 þ 0:5375 p08 ¼ 7Solving the system of pressure correction equations yields8 0< p3 ¼ 20p0 ¼ 20: 60p8 ¼ 20With the pressure correction computed, the velocity and pressure fields cannow be updated to produce a mass conserving velocity field.
The mass flowrates are computed as_ A þ m_ 0A ¼ m_ A 0:5p03 ¼ 50 0:5ð20Þ ¼ 60m_ A ¼mm_ B ¼ m_ B þ m_ 0B ¼ m_ B þ 0:4p03 ¼ 20 þ 0:4ð20Þ ¼ 28_ C þ m_ 0C ¼ m_ C 0:3p03 ¼ 75 0:3ð20Þ ¼ 69m_ C ¼m_ D þ m_ 0D ¼ m_ D þ 0:19 p03 p06 ¼ 19 þ 0:19ð20 þ 20Þ ¼ 19m_ D ¼m_ E þ m_ 0E ¼ m_ E þ 0:3p06 ¼ 30 þ 0:3ð20Þ ¼ 36m_ E ¼mm_ F ¼ m_ F þ m_ 0F ¼ m_ F 0:4p06 ¼ 48 0:4ð20Þ ¼ 40_ G þ m_ 0G ¼ m_ G þ 0:1875 p06 p08 ¼ 15 þ 0:1875ð20 þ 20Þ ¼ 15m_ G ¼m_ H þ m_ 0H ¼ m_ H 0:35p08 ¼ 28 0:35ð20Þ ¼ 35m_ H ¼mwhile the pressure is updated usingðnÞp3 ¼ p3 þ p03 ¼ 300 20 ¼ 280ðnÞp6 ¼ p6 þ p06 ¼ 200 20 ¼ 180ðnÞp8 ¼ p8 þ p08 ¼ 120 20 ¼ 100Treat the corrected values as a new guess and repeat. Better estimate for themass flow rates are computed using the momentum equations as15.2A Preliminary Derivation577ðnÞm_ A ¼ DA p1 p3 ¼ 0:5 ð400 280Þ ¼ 60ðnÞm_ B ¼ DB p3 p2 ¼ 0:4 ð280 350Þ ¼ 28ðnÞm_ C ¼ DC p4 p3 ¼ 0:3 ð50 280Þ ¼ 69ðnÞðnÞm_ D ¼ DD p3 p6 ¼ 0:19 ð280 180Þ ¼ 19ðnÞm_ E ¼ DE p6 p5 ¼ 0:3 ð180 300Þ ¼ 36ðnÞm_ F ¼ DF p7 p6 ¼ 0:4 ð80 180Þ ¼ 40ðnÞðnÞm_ G ¼ DG p6 p8 ¼ 0:1875 ð180 100Þ ¼ 15ðnÞm_ H ¼ DH p9 p8 ¼ 0:35 ð200 100Þ ¼ 35The imbalance in the mass flow rate at nodes 3, 6, and 8 are computed asXkXkXkm_ k ¼ m_ B þ m_ D m_ A m_ C ¼ 28 þ 19 60 þ 69 ¼ 0m_ k ¼ m_ G þ m_ E m_ D m_ F ¼ 15 36 19 þ 40 ¼ 0m_ k ¼ m_ I m_ G m_ H ¼ 50 15 35 ¼ 0The pressure correction equations become 1:39 p03 0:19 p06 ¼ 0 1:0775 p06 0:1875 p08 0:19 p03 ¼ 0 0 0 0:1875 p6 þ 0:5375 p8 ¼ 0The solution to the pressure correction field is found to be8 0< p3 ¼ 0p0 ¼ 0: 60p8 ¼ 0Thus the solution is obtained in one iteration.57815 Fluid Flow Computation: Incompressible Flows15.2.7 Pressure Correction Equation in Two DimensionalStaggered Cartesian GridsIn a two dimensional Cartesian grid, three grid systems are used.
One for the uvelocity component, a second one for the v-velocity component, and a third gridsystem for the pressure and other variables as illustrated in Fig. 15.5.The derivations presented above for the pressure correction equation in onedimensional domains can be easily extended into multi dimensional situations. Forelement C shown in Fig. 15.6, the pressure correction equation is obtained as000000apC p0C þ apE p0E þ apW p0W þ apN p0N þ apS p0S ¼ bpCð15:46Þ0qe Due DyCq Du DyCapW ¼ w wdxedxwvv00qDDxqDDxCCapN ¼ n napS ¼ s sdydy 0 n 0 sp0ppp0p0aC ¼ aE þ aW þ a N þ aS0bpC ¼ m_ e þ m_ w þ m_ n þ m_ sð15:47Þwhere0apE ¼ Fig. 15.5 u, v, andp elements in a twodimensional Cartesianstaggered gridmain CVpuvv velocity CVu velocity CV15.2A Preliminary Derivation579Fig. 15.6 A two dimensionalCartesian element for thederivation of the pressurecorrection equation( x )wNW( x )eNvn( y )nNEnueuwWw( y )svsSWeCE( y )CsSSE( x )CExample 2In the two dimensional problem shown in Fig. 15.7, the following quantitiesare given uw = 50, vs = 20, pN = 0 and pE = 10.The flow is steady and the density is uniform and equal to 1.
Themomentum equations for ue and vn are given byue ¼ de ðpE pC Þvn ¼ dn ðpN pC Þwhere the constants de = 1 and dn = 0.25. The element shown hasΔx = Δy = 1. Use the SIMPLE algorithm to compute the values of ue, vn, andpC .pN = 0uw = 50pE = 10vs = 20Fig. 15.7 The two dimensional domain used in Example 258015 Fluid Flow Computation: Incompressible FlowsSolutionStart by assigning a guessed value to the pressure at the C location wheresolution is to be found.
Thus assume the following:ðnÞpC ¼ 100 (other values could have been guessed)Based on the assumed pressure value calculate the various velocities usingthe momentum equation according toðnÞue ¼ de pE pC ¼ 1 ð10 100Þ ¼ 90ðnÞvn ¼ dn pN pC ¼ 0:25 ð0 100Þ ¼ 25Since the density is uniform and equal to 1 and Δx = Δy = 1, thenm_ e ¼ ue ¼ 90m_ n ¼ vn ¼ 25Check whetherP the mass flow rates satisfy continuity over element C bym_ f at the element faces.computingfXf m_ f ¼ m_ e m_ w þ m_ n m_ s ¼ 90 50 þ 25 20 ¼ 45In the above mass conservation equation the negative sign for m_ w and m_ sis explicitly used.
Since mass conservation is not satisfied, correction fieldsare needed and a pressure correction equation is derived as follows:Xfm_ f þ m_ 0f¼0)X X m_ 0f ¼ m_ fffIn term of pressure corrections, mass flow rate corrections can beexpressed asm_ 0e ¼ u0e ¼ p0Cm_ 0n ¼ v0n ¼ 0:25 p0CThe correction equation becomesXm_ 0e þ m_ 0n ¼ m_ f ) 1:25p0C ¼ 45 ) p0C ¼ 36fApplying the correction to the mass flow rate and pressure fields, continuity satisfying fields are obtained as15.2A Preliminary Derivation581_ e þ m_ 0e ¼ m_ e þ p0C ¼ 90 36 ¼ 54m_ e ¼ mm_ n ¼ m_ n þ m_ 0n ¼ m_ n þ 0:25p0C ¼ 25 0:25ð36Þ ¼ 160pC ¼ pC þ pC ¼ 100 36 ¼ 64Treat the corrected values as a new guess and repeat.m_ e ¼ de pE pC ¼ 1 ð10 64Þ ¼ 54m_ n ¼ dn pN pC ¼ 0:25 ð0 64Þ ¼ 16The imbalance in the mass flow rate is computed asXfm_ f ¼ m_ e m_ w þ m_ n m_ s ¼ 54 50 þ 16 20 ¼ 0The pressure correction equations become 1:25 p0c ¼ 0The solution to the pressure correction field is found to bep0C ¼ 0Thus the solution is obtained in one iteration asue ¼ 54vn ¼ 16pC ¼ 6415.2.8 Pressure Correction Equation in Three DimensionalStaggered Cartesian GridWithout going into details and for completeness of presentation, the pressure correction equation over the three dimensional staggered Cartesian grid shown inFig.
15.8, where the u, v, and w velocity components are stored at the (e, w), (n, s),and (t, b) element faces, respectively, is given by0000000apC p0C þ apE p0E þ apW p0W þ apN p0N þ apS p0S þ apT p0T þ apB p0B ¼ bCp0ð15:48Þ58215 Fluid Flow Computation: Incompressible FlowsTNynztxzxwWtnyweCbxeEszbysBSFig. 15.8 A three dimensional Cartesian element for the derivation of the pressure correctionequationwhere00qe Due Dye Dzeq Du Dyw DzwapW ¼ w wdxedxwqn Dvn Dxn Dznqs Dvs Dxs Dzsp0¼aS ¼ dyndysww0q D Dxt Dytq D Dxb Dyb¼ t tapB ¼ b bdztdzb 00pp0p0p0p0¼ aE þ aW þ aN þ aS þ aT þ apB¼ m_ e þ m_ w þ m_ n þ m_ s þ m_ t þ m_ bapE ¼ 0apN0apT0apC0bpCð15:49Þ15.3 Disadvantages of the Staggered GridThe use of staggered grids was critical to the development of the SIMPLE algorithm.
Nevertheless adopting a staggered grid arrangement has its disadvantages. Asmentioned above, in two and three dimensions, three and four staggered gridsystems, respectively, are required with the velocity components integrated overdifferent elements, as shown in Fig. 15.5 for a two dimensional situation.Besides the memory requirement to store a grid system for every velocitycomponent and a grid system for pressure and other variables, the staggeringprocedure itself becomes an issue for non-Cartesian grids and more so forunstructured grids.15.3Disadvantages of the Staggered Grid583Fig.