Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 46
Текст из файла (страница 46)
In the leastsquare methods, a gradient is computed by an optimization procedure that finds theminimum of the function GC defined asGC ¼NBðCÞX n2 owk /Fk ð/C þ r/C rCFk Þk¼1¼NBðCÞXk¼1( 2 ) ð9:23Þ @/@/@/wk D/k Dxkþ DykþDzk@x C@y C@z Cwhere wk is some weighting factor. The various terms in the above equationrepresentD/k ¼ /Fk /CDxk ¼ rCFk iDyk ¼ rCFk jð9:24ÞDzk ¼ rCFk kThe function GC is minimized by enforcing the conditions@G@G@G C ¼ C ¼ C ¼ 0@/@/@/@@@@x@y@zð9:25Þ9.3 Least-Square Gradient287to yield the following set of three equations in three unknowns: @/@/@/2wk Dxk D/k þ DxkþDykþDzk¼0@x@y@z CCCk¼1 NBðCÞX @/@/@/2wk Dyk D/k þ DxkþDykþDzk¼ 0 ð9:26Þ@x C@y C@z Ck¼1 NBðCÞX @/@/@/2wk Dzk D/k þ DxkþDykþDzk¼0@x@y@z CCCk¼1NBðCÞX which can be written in matrix form as2NBðCÞPNBðCÞPwk Dxk Dxkwk Dxk Dyk66 k¼1k¼16 NBðCÞNBðCÞ6 PP6wk Dyk Dxkwk Dyk Dyk66 k¼1k¼16 NBðCÞNBðCÞP4 Pwk Dzk Dxkwk Dzk Dykk¼1k¼132NBðCÞPwk Dxk D/k 7676 k¼176 NBðCÞ76 P7¼6wDyD/kkk7676 k¼176 NBðCÞ54 Pwk Dzk D/k32 3@/wk Dxk Dzk 7676 @x C 7k¼176 7NBðCÞ76 @/ 7P77wk Dyk Dzk 7766@y76 C 7k¼177NBðCÞP54 @/ 5wk Dzk Dzk@z Ck¼1NBðCÞPð9:27Þk¼1The solution to the above set of equations yields the gradient ðr/ÞC .
A solutionexists provided that the matrix on the left hand side is not singular. Moreover, thechoice of the wk is important in determining the properties of the gradient. Forexample if wk is chosen to be 1 for all neighbors of C, then all neighboring pointswill have the same weight in the computation of the gradient irrespective of whetherthey are near or far from point C. Actually points that are farther from C will have amore important influence as the error function will be more affected by their error.Another choice for wk, which was used earlier with the extended stencil method,is the inverse distance between C and F given bywk ¼11¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffijrFk rC jDx2 þ Dy2 þ Dz2FkFkFkð9:28Þ2889Fig.
9.9 A three dimensionalCartesian control volumeGradient ComputationTNynztxzWxwtnyweCbxeEszbysBSOther options that can be pursued include the inverse distance raised to anypower n such thatwk ¼1jrFk rC jnð9:29Þwhere n can be set to 1, 2, 3, etc.As mentioned above, the divergence based gradient is a special case of theleast-square formulation. This can be shown for a Cartesian grid (see Fig.
9.9)where substituting the geometric quantities into Eq. (9.27) would yield the following set of three equations.2xE xW4000yN yS0323 23ð@/=@xÞC0/E /W0 54 ð@/=@yÞC 5 ¼ 4 /N /S 5ð@/=@zÞCzT zB/T /Bð9:30ÞSolving the above equation, the derivatives in the various directions are found to be @//E /W@//N /S@// /B¼¼¼ Tð9:31Þ@x C xE xW@y C yN yS@z C zT zBThe obtained values are exactly the ones given in Eq. (9.3), demonstrating that thedivergence-based gradient is a special case of the least-square method.
Finally it caneasily be demonstrated that the accuracy of the resulting gradient is at least first order.Indeed, the Taylor series expansion of the / value around node C can be written as /ðrÞ /ðrC Þ ¼ ðr/ÞC ðr rC Þ þ O r2which when solved for ðr/ÞC results in O(r).ð9:32Þ9.4 Interpolating Gradients to Faces9.4289Interpolating Gradients to FacesIt was shown in Chap.
8 that the discretization of the diffusion term innon-orthogonal grids requires the use of correction terms involving gradients atcontrol volume faces. Thus in this situation the gradients need to be interpolated fromthe control volume centroids where they were computed to the control volume faceswhere they will be used. Figure 9.10a shows the stencil used in the Gauss gradientcomputation for two neighboring control volumes. The gradient at the face will haveexactly the same stencil, while ideally it should be similar to that of Fig. 9.10b.A better insight is gained by considering the configuration in Fig. 9.11a, whichshows the gradients r/C and r/F of the variable / at the two nodes C and F,respectively.
The interpolated gradient at the face, r/f , is obtained by averagingthe values at nodes C and F, as shown in Fig. 9.11b. It is important for the stencil ofthe gradient at the face to be heavily based on the nodes straddling the face, whichis not guaranteed by this simple averaging practice. As schematically displayed inFig. 9.11c, this can be accomplished by forcing the face gradient along the CFdirection to be equal to the local gradient defined by the values of / at C andF. Mathematically this can be written asr/f ¼ r/f þ/F /C r/f eCF eCFdCF|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}ð9:33ÞCorrection interpolated face gradientwherer/f ¼ gC r/C þ gF r/F ;(a)eCF ¼dCF;dCFdCF ¼ rF rC(b)CffFeCF stencilCffFeCF stencilFig.
9.10 a Interpolated and b corrected gradient at a control volume faceð9:34Þ2909Fig. 9.11 a Schematics of thegradients at the two nodesC and F straddling face f;b computing r/f as a simpleaverage of r/C and r/Fusing Eq. (9.31); c thegradient at the face with itsvalue heavily based on thenodes straddling the faceGradient Computation(a)()F()fSf()CCf(()C)f()F((c))fd CFeCF(b)(Fn(()f))feCF eCFFeCFdCFCeCFwhere rF and rC are position vectors, as displayed in Fig.
9.3. This approach isapplicable to both structured and unstructured grids. For unstructured grids, whilethe stencil used for the face gradient might not be decreased, the gradient across aface will still be based on the nodes straddling that face.9.59.5.1Computational PointersuFVMIn uFVM, the functions cfdComputeGradientGauss0 and cfdComputeGradientNodal are used to compute the element gradients. The Gauss gradient routine(Listing 9.1) takes as input the phi element array and returns the element gradientarray. The values are interpolated to the faces using a weighted factor with nocorrection for non-conjunctionality, hence the 0 digit in the function name.9.5 Computational Pointers291function phiGrad = cfdComputeGradientGauss0(phi,theMesh)%===================================================% written by the CFD Group @ AUB, Fall 2006%===================================================%if(nargin<2)theMesh = cfdGetMesh;endtheSize = size(phi);theNumberOfComponents = theSize(2);if(theNumberOfComponents > 3)echo('********* ERROR **********');exit;end%----------------------------------------------------% INTERIOR FACES contribution to gradient%----------------------------------------------------iFaces = 1:theMesh.numberOfInteriorFaces;iBFaces = theMesh.numberOfInteriorFaces+1:theMesh.numberOfFaces;iElements = 1:theMesh.numberOfElements;iBElements=theMesh.numberOfElements+1:theMesh.numberOfElements+theMesh.numberOfBFaces;iOwners = [theMesh.faces(iFaces).iOwner]';iNeighbours = [theMesh.faces(iFaces).iNeighbour]';Sf = [theMesh.faces(iFaces).Sf]';gf = [theMesh.faces(iFaces).gf]';%----------------------------------------------------% Initialize phiGrad Array%----------------------------------------------------phiGrad=zeros(theMesh.numberOfElements+theMesh.numberOfBElements,3,theNumberOfComponents);for iComponent=1:theNumberOfComponentsphi_f = gf.*phi(iNeighbours,iComponent) + (1-gf).*phi(iOwners,iComponent);%for iFace=iFaces=phiGrad(iOwners(iFace),:,iComponent)phiGrad(iOwners(iFace),:,iComponent)+ phi_f(iFace)*Sf(iFace,:);phiGrad(iNeighbours(iFace),:,iComponent)=phiGrad(iNeighbours(iFace),:,iComponent) - phi_f(iFace)*Sf(iFace,:);endend%----------------------------------------------------% BOUNDARY FACES contribution to gradient%----------------------------------------------------iBOwners = [theMesh.faces(iBFaces).iOwner]';phi_b = phi(iBElements,iComponent);Sb = [theMesh.faces(iBFaces).Sf]';for iComponent=1:theNumberOfComponents%for k=1:theMesh.numberOfBFacesphiGrad(iBOwners(k),:,iComponent)phiGrad(iBOwners(k),:,iComponent)+ phi_b(k)*Sb(k,:);endend=Listing 9.1 Function used in uFVM to compute the gradient field at the centroids of elementsfollowing the Green-Gauss approach with phi values at the face computed using simple a weightedaverage interpolation technique with no correction to non-conjunctionality2929Gradient Computation%----------------------------------------------------% Get Average Gradient by dividing with element volume%----------------------------------------------------volumes = [theMesh.elements(iElements).volume]';for iComponent=1:theNumberOfComponentsfor iElement =1:theMesh.numberOfElementsphiGrad(iElement,:,iComponent) = phiGrad(iElement,:,iComponent)/volumes(iElement);endend%----------------------------------------------------% Set boundary Gradient equal to associated element% Gradient%----------------------------------------------------phiGrad(iBElements,:,:) = phiGrad(iBOwners,:,:);endListing 9.1 (continued)The nodal gradient routine, listed below (Listing 9.2), is very similar to theGauss gradient routine except that it uses the functions cfdInterpolateFromElementsToNodes and cfdInterpolateFromNodesToFaces to interpolate phi valuesto the face that are subsequently used in the Gauss algorithm.%%=====================================================% INTERIOR Gradients%=====================================================theNumberOfElements = theMesh.numberOfElements;theNumberOfBElements = theMesh.numberOfBElements;theNumberOfInteriorFaces = theMesh.numberOfInteriorFaces;%phiNodes = cfdInterpolateFromElementsToNodes(phi);phi_f = cfdInterpolateFromNodesToFaces(phiNodes);%phiGrad = zeros(3,theNumberOfElements+theNumberOfBElements);%----------------------------------------------------% INTERIOR FACES contribution to gradient%----------------------------------------------------fvmFaces = theMesh.faces;% interpolate phi to faces%for iFace=1:theNumberOfInteriorFaces%theFace = fvmFaces(iFace);%iElement1 = theFace.iOwner;iElement2 = theFace.iNeighbour;%Listing 9.2 Function used in uFVM to compute the gradient field at the centroids of elementsfollowing the Green-Gauss approach with phi values at the face computed using nodal values, i.e.,the extended stencil approach9.5 Computational Pointers293Sf = theFace.Sf;%%phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_f(iFace)*Sf;phiGrad(:,iElement2) = phiGrad(:,iElement2) - phi_f(iFace)*Sf;end%=====================================================% BOUNDARY FACES contribution to gradient%=====================================================for iBPatch=1:theNumberOfBElements%iBFace = theNumberOfInteriorFaces+iBPatch;iBElement = theNumberOfElements+iBPatch;theFace = fvmFaces(iBFace);%iElement1 = theFace.iOwner;%Sb = theFace.Sf;phi_b = phi(iBElement);%phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_b*Sf;end%----------------------------------------------------% Get Average Gradient by dividing with element volume%----------------------------------------------------for iElement =1:theNumberOfElementstheElement = fvmElements(iElement);phiGrad(:,iElement) = phiGrad(:,iElement)/theElement.volume;end%----------------------------------------------------% Set boundary Gradient equal to associated element% Gradient%----------------------------------------------------for iBPatch = 1:theNumberOfBElementsiBElement = iBPatch+theNumberOfElements;iBFace = iBPatch+theNumberOfInteriorFaces;theBFace = fvmFaces(iBFace);iOwner = theBFace.iOwner;phiGrad(:,iBElement) = phiGrad(:,iOwner);endListing 9.2 (continued)For face interpolation of gradients, a variety of interpolation options are allowedin the function cfdInterpolateGradientsFromElementsToInteriorFaces.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
