Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 45
Текст из файла (страница 45)
Combining Eqs. (9.10) and (9.11) yieldsrf 0 ¼rf neenð9:12Þwith f 0 located, gC is computed asrF rf 0 dFf 0gC ¼¼krF rC k dFCð9:13Þ2789Gradient ComputationThen, the calculation procedure involves the following steps:During the first iteration, calculate the gradient field over the entire domain asfollows:1. Calculate /f 0 using/f 0 ¼ gC /C þ ð1 gC Þ/F1 X2.
Calculate r/C using r/C ¼/ 0 SfVC f nbðCÞ fFrom the second iteration onward, correct the gradient field according to thefollowing procedure:/f ¼ /f 0 þ gC ðr/ÞC rf rC þ ð1 gC Þðr/ÞF rf rF3. Update /f using4. Update r/C using r/C ¼1 X/ SfVC f nbðCÞ f5. Go back to step 3 and repeat.Option 2: For f 0 chosen to be at the centre of segment [CF] [Fig. 9.4a in twodimensions and Fig. 9.4b in three dimensions] the equations become simpler andthe calculation of the gradient field over the domain proceeds as follows:During the first iteration, calculate the gradient field over the entire domain asfollows:/C þ /F21 X2.
Calculate r/C using r/C ¼/ 0 SfVC f nbðCÞ f1. Calculate /f 0 using/f 0 ¼From the second iteration onward, correct the gradient field according to thefollowing procedure: 3. Update /f using/f ¼ /f 0 þ 0:5 ðr/ÞC þ ðr/ÞF rf 0:5 ðrC þ rF Þ1 X4. Update r/C using r/C ¼/ SfVC f nbðCÞ f5.
Go back to step 3 and repeat(b)(a)SfCffFfFfeCSfFig. 9.4 Correction to Non-Conjunctionality using the midpoint approach: a two dimensionalconfiguration; b three dimensional configuration9.2 Green-Gauss Gradient279(b)(a)SfCffFfFfeCSfFig. 9.5 Correction to Non-Conjunctionality using the minimum distance approach: a twodimensional configuration; b three dimensional configurationOption 3: The position of f 0 can be chosen such that the distance ff 0 is the shortestpossible [Fig. 9.5a in two dimensions and Fig. 9.5b in three dimensions], i.e., ½ff 0 perpendicular to [CF].
This leads to a more accurate computation of the gradientduring the first iteration. In this case f 0 is computed by minimizing the distancebetween f and f 0 . In general rf 0 is described byr f 0 ¼ r C þ qð r C r F Þð9:14Þwhere 0 < q < 1.Denoting the distance f 0 f by d, then its square is obtained as d 2 ¼ rf rf 0 rf rf 0 ð9:15Þ¼ rf rC qðrC rF Þ rf rC qðrC rF Þ 2¼ rf rC rf rC 2q rf rC ðrC rF Þ þ q ðrC rF Þ ðrC rF ÞMinimizing the function d2 with respect to q, yields@ ðd 2 Þ¼ 0 ) 2 rf rC ðrC rF Þ þ 2qðrC rF Þ ðrC rF Þ ¼ 0@qð9:16ÞSolving, q is obtained asq¼rCf rCFrCF rCFð9:17ÞKnowing the q values, the gradient calculation over the domain proceeds asfollows:During the first iteration, calculate the gradient field over the entire domain asfollows:1. Calculate rf 0 using2.
Calculate gC usingrCf rCFrf 0 ¼ rC ð rC rF Þr rCFCF gC ¼ rF rf 0 =krF rC k28093. Calculate /f 0 using4. Calculate r/C usingGradient Computation/f 0 ¼ gC /C þ ð1 gC Þ/F1 Xr/C ¼/ 0 SfVC f nbðCÞ fFrom the second iteration onward, correct the gradient field according to thefollowing procedure:5. Calculate r/f 0 using6. Update /f using7.
Update r/C usingr/f 0 ¼ gC r/C þ ð1 gC Þr/F/f ¼ /f 0 þ r/f 0 rf rf 01 Xr/C ¼/ SfVC f nbðCÞ f8. Go back to step 5 and repeatExample 1For the mesh shown in Fig. 9.6, the coordinates of the grid point C and itsneighbors F1 through F6 are given byCð13; 11ÞF1 ð4:5; 9:5ÞF2 ð8; 3ÞF3 ð17; 3:5ÞF4 ð22; 10Þ F5 ð16; 20ÞFig.
9.6 Grid layout forExamples 1 and 2F6 ð7; 18ÞyF5F6n6n1f6f1Sf1F1f5n2f4Cf1′n5f3f2n4F4n3F2F3xwhile the nodes n1 through n6 are located atn1 ð9; 14Þn4 ð17; 9Þn2 ð8; 8Þn3 ð12; 5Þn5 ð17:5; 14Þ n6 ð12; 17Þ9.2 Green-Gauss Gradient281If the values of the dependent variable / at the centroids are known to be/C ¼ 167/F1 ¼ 56:75/F4 ¼ 252/F2 ¼ 35/F3 ¼ 80/F5 ¼ 356 /F6 ¼ 151and the values of the gradient of /, ðr/Þ, at all neighboring elements to Care given byr/F1 ¼ 10:5i þ 5:5j r/F2 ¼ 4i þ 9j r/F3 ¼ 4:5i þ 18jr/F4 ¼ 11i þ 23jr/F5 ¼ 21i þ 17jr/F6 ¼ 19i þ 8jIf the volume of cell C is VC = 76, find the gradient r/C usinga.
The Green-Gauss method with no correctionb. The Green-Gauss method alongside correction to skewness with f 0 chosento be at the centre of segment [CF].Solution The Green-Gauss method with no correctiona. The interpolation factors are needed to perform the calculations. This, inturn, necessitates computing the coordinates of the face centroids. Thecoordinates of the centroid f1 are found asxf1 ¼ 0:5 ðxn1 þ xn2 Þ ¼ 0:5 ð9 þ 8Þ ¼ 8:5yf1 ¼ 0:5 ðyn1 þ yn2 Þ ¼ 0:5 ð14 þ 8Þ ¼ 11) f1 ð8:5; 11ÞIn a similar way, the coordinates of other face centroids are found to bef2 ð10; 6:5Þf3 ð14:5; 7Þf4 ð17:25; 11:5Þf5 ð14:75; 15:5Þf6 ð10:5; 15:5ÞThe surface vectors are calculated asSf1 ¼ 6i þ jSf4 ¼ 5i 0:5jSf2 ¼ 3i 4jSf5 ¼ 3i þ 5:5jSf3 ¼ 4i 5jSf6 ¼ 3i þ 3j2829Gradient ComputationThe interpolation factor (gC)1 is computed using9F1 f1>>>>F1 f1 þ Cf1>=qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi22F1 f1 ¼ ð4:5 8:5Þ þð9:5 11Þ ¼ 4:272 > ) ðgC Þ1 ¼ 0:487>>qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi>>;22Cf1 ¼ ð13 8:5Þ þð11 11Þ ¼ 4:5ð gC Þ 1 ¼In a similar way, the other interpolation factors are found to beðgC Þ2 ¼ 0:427 ðgC Þ3 ¼ 0:502 ðgC Þ4 ¼ 0:538 ðgC Þ5 ¼ 0:492 ðgC Þ6 ¼ 0:455Using Eq.
(9.5) the /f values are computed as/f1 ¼ 110:442/f2 ¼ 91:364/f4 ¼ 206:27/f5 ¼ 263:012 /f6 ¼ 158:28/f3 ¼ 123:674Using Eq. (9.4), r/C is calculated as8 9110:442 ð6Þ þ 91:364 ð3Þ þ 123:674 4 þ>>>i þ>>>=206:27 5 þ 263:012 3 þ 158:28 ð3Þ1 <r/C ¼76 >110:442 ð1Þ þ 91:364 ð4Þ þ 123:674 ð5Þ þ >>>>j>;:206:27 ð0:5Þ þ 263:012 5:5 þ 158:28 ð3Þ¼ 11:889i þ 12:433jb. The Green-Gauss method alongside correction to skewness with f 0 chosento be at the centre of segment [CF].The values at the f 0 locations are computed as half the sum of the values atthe nodes straddling the face and are given by/f1 ¼ 111:875/f2 ¼ 101/f3 ¼ 123:5/f4 ¼ 209:5/f5 ¼ 261:5/f6 ¼ 158:5Using these values, the first estimate for the gradient is obtained usingEq. (9.4) asr/C ¼ 11:53i þ 11:826j9.2 Green-Gauss Gradient283Defining df = rf − 0.5 × (rC + rF), the various values are found to bedf1 ¼ 0:25i þ 0:75j df2 ¼ 0:5i 0:5j df3 ¼ 0:5i 0:25jdf4 ¼ 0:25i þ jdf5 ¼ 0:25idf6 ¼ 0:5i þ jUsing /f ¼ /f 0 þ 0:5 ½ðr/ÞC þ ðr/ÞF rf 0:5 ðrC þ rF Þ theupdated values at the faces are obtained as/f1 ¼ 115:619/f2 ¼ 91:911/f4 ¼ 224:097/f5 ¼ 265:566 /f6 ¼ 176:046/f3 ¼ 115:764These values are used in Eq.
(9.4) yielding the updated value of thegradient asr/C ¼ 11:614i þ 13:761jMethod 2: Extended Stencil [2]The value of /f at the surface centroid f can be computed as the mean of thevalues at the vertices defining the surface. This necessitates the estimation ofthe properties at the vertices. The properties at a vertex node are calculated usingthe weighted average of the properties within the cells surrounding that node.Figure 9.7 shows the cells that are considered for the weighted average of theproperties at the vertex nodes n1 and n2. The weight is taken as the inverse ofthe distance of the vertex from the cell centre [3]. The resulting equation for theproperties at the vertices are written as,Fig.
9.7 Cells Contributingto node n for the WeightedAverage2849NBðnÞP/n ¼k¼1NBðnÞPk¼1Gradient Computation/Fkkr n r F k k1kr n r F k kð9:18Þwhere n refers to the vertex node, Fk to the neighboring cell node, NB(n) the totalnumber of cell nodes surrounding the vertex node n, and krn rFk k the distancefrom the vertex node to the centroid of the neighboring cells.Once the values /n at the vertices are found, the values /f at the surfacecentroids are calculated followed by the gradients at the control volume centroids.In two dimensional situations, /f is computed as/f ¼/n1 þ /n22ð9:19ÞThen the gradient at C is found usingr/C ¼/n1 þ /n21 X1 X/f Sf ¼SfVC f ¼nbðCÞVC f nbðCÞ2fð9:20ÞIn three dimensional situations, the calculations are a little more involved as thenumber of a face vertices depends on the element type.
The value of /f is foundfrom the values at the vertices usingnbðfPÞ/f ¼k¼1nbðfPÞk¼1/nkrn rf k1r n r f kð9:21Þwhere n represents the number of vertices of face f. Once the values /f are calculated, the gradient at C is computed using Eq. (9.4).One of the disadvantages of this approach is that information from the wrongside of the cell face also contributes to the weighted average values of the conservedvariables. This can be overcome by using upwind biased gradients as discussed byCabello [4].
The higher order calculations based on the upwind biased gradients,however, have both higher memory overheads required to store the informationabout the cells used for the upwind biased gradient calculation and increased codingcomplexity.9.2 Green-Gauss Gradient285Example 2Using the data of example I, calculate /f1 using the extended stencilapproach via Eq. (9.16).Solution First the distances have to be calculated and are given byqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffið9 7Þ2 þð14 18Þ2 ¼ 4:472qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffikrn1 rF1 k ¼ ð9 4:5Þ2 þð14 9:5Þ2 ¼ 6:364qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffikrn1 rC k ¼ ð9 13Þ2 þð14 11Þ2 ¼ 5qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffikrn2 rF1 k ¼ ð8 4:5Þ2 þð8 9:5Þ2 ¼ 3:808qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffikrn2 rF2 k ¼ ð8 8Þ2 þð8 3Þ2 ¼ 5qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffikrn2 rC k ¼ ð8 13Þ2 þð8 11Þ2 ¼ 5:831krn1 rF6 k ¼The values at nodes n1 and n2 are computed using Eq.
(9.15) as15156:75 16756:75 35167þþþ þ5 ¼ 131:009 / ¼ 3:808 5 5:831 ¼ 79:708/n1 ¼ 4:472 6:364n2111111þþþ þ4:472 6:364 53:808 5 5:831The value of /f1 is found to be/f1 ¼ 0:5ð131:009 þ 79:708Þ ¼ 105:35859.3Least-Square GradientUsing least-square methods to compute the gradients [5] offers more flexibility withregard to the order of accuracy achieved [6] and the stencil used [7]. In theleast-square method, the divergence-based gradient can be recovered as a specialcase. This flexibility comes at a cost, as proper weighting is needed for the stencilterms, and computation of the weights adds to the computational cost.
The methodis described next.Considering a control volume and its immediate neighbors (Fig. 9.8), the changein centroid values between C and F is given by ð/F /C Þ, if the cell gradient(r/C ) is exact, then the difference can also be computed as2869Gradient ComputationFig. 9.8 A control volume with its immediate neighbors/F ¼ /C þ ðr/ÞC ðrF rC Þ|fflfflfflfflfflffl{zfflfflfflfflfflffl}ð9:22ÞrCFHowever unless the solution field is linear the cell gradient cannot be exactbecause C has more neighbors than the gradient vector has components.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
