Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab.pdf), страница 96

Startingwith the mass flow rate field (i.e., the mDot flux), it is updated by executing thebelow statement (Listing 15.14).15.10Computational Pointers645mDot += ppEqn.flux();Listing 15.14 Statement used to update the mass flow rate field to satisfy continuityThe flux() function in Listing 15.14 provides the matrix multiplication, the extradiagonal matrix coefficients, and the corresponding solution.

Recalling the finitevolume discretization of the diffusion term, each extra diagonal of coefficientsrepresents a face of the mesh. Thus the update of the fluxes can be performed in amore consistent way using directly the matrix coefficients and cell values.A simplified version of the flux() function is shown in Listing 15.15.for (label face=0; face<lowerAddr.size(); face++){mDotPrime[face] =upperCoeffs[face]*pp[upperAddr[face]]- lowerCoeffs[face]*pp[lowerAddr[face]];}return mDotPrime;Listing 15.15 A simplified version of the flux() function where the flux correction mDotPrime iscomputedIn Listing 15.15 the correction flux mDotPrime (Eq.

15.101) is basically evaluated by performing a loop over the faces using the upper and lower coefficients ofthe matrix and multiplying these coefficients with the corresponding cell values.Finally the velocity and pressure at cell centroids are updated using the scriptshown in Listing 15.16,scalar URF = mesh.equationRelaxationFactor("pp");p += URF*pp;p.correctBoundaryConditions();U -= fvc::grad(pp)*DU;U.correctBoundaryConditions();Listing 15.16 Update of the velocity and pressure fields at cell centroidswhere the variable URF is the explicit relaxation factor for pressure update λp,necessary for a stable SIMPLE solver.64615 Fluid Flow Computation: Incompressible FlowssimpleFoamImprovedIn simpleFoamImproved the Rhie-Chow interpolation is extended to account forthe relaxation of the velocity field.

This is translated into the syntax presented belowin Listing 15.17 that expands the generic Rhie-Chow interpolation to include theadditional term in Eq. (15.196).// Rhie-Chow interplationmdotf =(U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f)) &mesh.Sf() )+(scalar(1) - URFU)*(mdotf.prevIter() - (U_avg_prevIter_f &mesh.Sf()));Listing 15.17 Improved Rhie-Chow interpolation accounting for under relaxationThus the fields mdotf.prevIter() and U.prevIter() need to be defined.simpleFoamTransientIn simpleFoamTransient the Rhie-Chow interpolation is extended to account forthe transient term.

Thus the expression for the mass flow rate becomes (Listing15.18).// Rhie-Chow interplationmdotf =(U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f)) &mesh.Sf() )+(scalar(1) - URFU)*(mdotf.prevIter() - (U_avg_prevIter_f &mesh.Sf()))+ DTf*( mdotf_old - (U_old_f& mesh.Sf()));Listing 15.18 Rhie-Chow interpolation accounting for the effects of under-relaxation and theunsteady termFurthermore the main loop is modified to add a transient loop, with the maincode becoming as shown in Listing 15.19.15.10Computational Pointers647pimpleControl pimple(mesh);// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** * * //Info<< "\nStarting time loop\n" << endl;while (runTime.run()){#include "readTimeControls.H"#include "CourantNo.H"#include "setDeltaT.H"runTime++;Info<< "Time = " << runTime.timeName() << nl << endl;scalar iter=0;while (pimple.loop()){iter++;Info<< "Iteration = " << iter << nl << endl;//U.storePrevIter();mdotf.storePrevIter();// Pressure-velocity SIMPLE corrector#include "UEqn.H"#include "ppEqn.H"}runTime.write();Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"<< " ClockTime = " << runTime.elapsedClockTime() << " s"<< nl << endl;}Listing 15.19 The main loop used for solving unsteady flow problemssimpleFoamBuoyancyThe simpleFoamBuoyancy solver adds to simpleFoamTransient the followingcapabilities: (i) the solution of the energy equation, (ii) the inclusion of a bodysource term in the momentum equation, and (iii) an account of the effects of thebody force term redistribution in the Rhie-Chow interpolation.The codes used to introduce these modifications are shown in Listings (15.20),(15.21), and (15.22).64815 Fluid Flow Computation: Incompressible FlowsThe code needed to solve the energy equation is given in Listing (15.20).// Solve the Energy equationfvScalarMatrix TEqn(fvm::ddt(T)+ fvm::div(phi, T)- fvm::laplacian(K, T));TEqn.relax();TEqn.solve();Listing 15.20 The script used to solve the energy equationThe source term in the momentum equation is implemented as (Listing 15.21),surfaceVectorField B_f = linearInterpolate(- beta*(T-To)*g);volVectorField B_reconstructed = fvc::average(B_f);solve(UEqn == -fvc::grad(p) + B_reconstructed);Listing 15.21 Accounting for the body force term source in the momentum equationwhile the calculation of the mass fluxes using the Rhie-Chow interpolation are asdisplayed in Listing (15.22).surfaceVectorField B_reconstructed_f =linearInterpolate(B_reconstructed);// Rhie-Chow interplationphi = (U_avg_f & mesh.Sf())- DUf*( (gradp_f*mesh.magSf())-(gradp_avg_f&mesh.Sf()))+ (scalar(1) - URFU)*(phi.prevIter() - (U_avg_prevIter_f &mesh.Sf()))+ DTf*(phi_old - (U_old_f& mesh.Sf()))+ DUf* ( (B_f& mesh.Sf()) - (B_reconstructed_f& mesh.Sf()) ) ;Listing 15.22 The modified Rhie-Chow interpolation accounting for body forces15.11Closure64915.11 ClosureThis chapter presented the segregated pressure-based approach for solvingincompressible flow problems on collocated grids.

It also demonstrated that thesuccess of the Rhie-Chow interpolation on collocated grids is due to its formation ofa pseudo-momentum equation at the cell face that has a tight pressure gradientstencil similar to the one resulting from a staggered grid formulation.

In addition,the details of implementing the most commonly encountered boundary conditionsin the momentum and pressure-correction equations were discussed. The nextchapter will extend the pressure based method to predict compressible fluid flow atall speeds.15.12 ExercisesExercise 1A portion of a water-supply system is shown in Fig.

15.31. The flow rate m_ in apipe section is given bym_ ¼ CDpwhere Δp is the pressure drop over the length of the pipe section, and C is thehydraulic conductance. The following data is known:p1 ¼ 400; p2 ¼ 350m_ F ¼ 25CA ¼ 0:4; CB ¼ 0:2; CC ¼ 0:1; CD ¼ 0:3; CE ¼ 0:2Find p3, p4, p5, m_ A ; m_ B ; m_ C ; m_ D and m_ E using the following procedure••••Start with a guess for p3, p4, and p5.Compute m_ values based on the guessed pressures.Construct the pressure-correction equations and solve for p03 , p04 and p05 .Update the pressures and the m_ values65015 Fluid Flow Computation: Incompressible FlowsDo you need to iterate? Why?p2 = 350mBp3mEmF = 25mCp5mDp1 = 400p4mAFig.

15.31 A portion of a water-supply systemExercise 2A one dimensional flow through a porous material is governed bycjuju þdp¼0dxwhere c is a constant. The continuity equation isdu¼0dxx2 x1 ¼ 1SA ¼ 3x3 x2 ¼ 2SB ¼ 2Use the SIMPLE procedure for the grid shown in Fig.

15.32 to compute p2, uA,and uB from the following data:cA ¼ 0:3p1 ¼ 150cB ¼ 0:15p3 ¼ 18with the size and area at the center of each control volume given byDxA ¼ 1; AA ¼ 3DxB ¼ 3; AB ¼ 21uA2Fig. 15.32 One dimensional flow in a porous materialuB315.12Exercises651Exercise 3In the Steady, one dimensional, constant density situation shown in Fig.

15.33, thevelocity u is calculated at locations B and C, while the pressure is calculated atlocations 1 and 2. The velocity correction formulae are written asuB ¼ DB ðp1 p2 ÞanduC ¼ DC ðp2 p3 Þwhere the values of DB and DC are 3 and 4, respectively. The boundary conditionsare uA = 5 and p3 = 70.(a) If at a given stage in the iteration process, the momentum equations giveuB ¼ 4 and uC ¼ 6, calculate the values of p1 and p2.(b) Explain how you could obtain the values of p1 and p3 if the right handboundary condition is given as uC = 5 instead of p3 = 70.uA1uBuC23Fig. 15.33 Incompressible flow in a one dimensional domainExercise 4Consider the main control volume shown in Fig.

15.34. A staggered mesh is usedwith the u and v velocity components stored as shown. The following quantities aregiven: uw = 7, vs = 3, pN = 0 and pE = 50. The flow is steady and the density isconstant. The momentum equations for ue and vn are given by:ue ¼ De ðpE pC Þvn ¼ Dn ðpN pC ÞAlso given De = 2, Dn = 1.6, and the control volume has Δ x = Δy = 1.ðnÞ(a) Starting with a guessed value of pC ¼ 50, use the SIMPLE algorithm to findue and vn.(b) Is an iteration loop needed for this problem? Explain.pNvnpEpCuwuevsFig.

15.34 A main control volume in a two-dimensional staggered grid arrangement65215 Fluid Flow Computation: Incompressible FlowsExercise 5Consider the simplified one-dimensional Forchheimer model for flow in porousmedia given bybu2 ¼ kdpdxwith the continuity equation given bydðeuÞ¼0dxIn the above equations b is a constant and ε is the porosity coefficient thataccounts for the effective porous area.Devise a SIMPLE-like procedure to compute pC, ue, and uw for the following data:Dx ¼ 0:1;Dy ¼ 1bW ¼ 5;ee ¼ 0:9;bC ¼ 4; bE ¼ 3ew ¼ 0:6pW ¼ 40;pE ¼ 200Start with the following initial values for velocity and pressure (Fig. 15.35):ue ¼ uw ¼ 3 and pC ¼ 100:xuwpWWueypCpWpEEpEFig.