Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 33
Текст из файла (страница 33)
Function cfdComputeGradientGauss0, displayed inListing 7.34, shows the details of the Gauss Gradient implementation.7.1 uFVMfunction phiGrad = cfdComputeGradientGauss0(phi,theMesh)%===================================================% written by the CFD Group @ AUB, Fall 2014%===================================================%if(nargin<2)theMesh = cfdGetMesh;end%----------------------------------------------------% Initialize phiGrad Array%----------------------------------------------------phiGrad = zeros(theMesh.numberOfElements+theMesh.numberOfBFaces,3);%----------------------------------------------------% INTERIOR FACES contribution to gradient%----------------------------------------------------% get the list of interior faces indices and the list of boundary facesindicesiFaces = 1:theMesh.numberOfInteriorFaces;iBFaces = theMesh.numberOfInteriorFaces+1:theMesh.numberOfFaces;% get the list of element indices and the list of boundary element indicesiElements = 1:theMesh.numberOfElements;iBElements = theMesh.numberOfElements+1:theMesh.numberOfElements+theMesh.numberOfBFaces;% get the list of owners and neighbours for all interior facesiOwners = [theMesh.faces(iFaces).iOwner]';iNeighbours = [theMesh.faces(iFaces).iNeighbour]';% get the surface vector of all interior facesSf = [theMesh.faces(ifaces).Sf]’;% get the geometric factor for all interior facesgf = [theMesh.faces(iFaces).gf]';% compute the linear interpolation of phi into all interior facesphi_f = gf.*phi(iNeighbours,iComponent) + (1-gf).*phi(iOwners,iComponent);% loop over faces and add contribution of phi face flux to the owner andneighbour elements of the facesfor iFace=iFacesphiGrad(iOwners(iFace),:) = phiGrad(iOwners(iFace),:) +phi_f(iFace)*Sf(iFace,:);phiGrad(iNeighbours(iFace),:) = phiGrad(iNeighbours(iFace),:)phi_f(iFace)*Sf(iFace,:);end%----------------------------------------------------% BOUNDARY FACES contribution to gradient%-----------------------------------------------------% get the list ofelements owning a boundary faceiOwners_b = [theMesh.faces(ibFaces).iOwner]’;% get the boundary values of phiphi_b = phi(iBElements,iComponent);% get the surface vector of all boundary facesSb = [theMesh.faces(iBFaces).Sf]';% loop over all boundary faces and add phi flux contribution to gradientfor k=1:theMesh.numberOfBFacesphiGrad(iOwners_b(k),:) = phiGrad(iOwners_b(k),:) + phi_b(k)*Sb(k,:);endListing 7.34 Computing the Gauss gradient1897 The Finite Volume Mesh in OpenFOAM® and uFVM190%----------------------------------------------------% Get Average Gradient by dividing with element volume%% get volume of all elements in meshvolumes = [theMesh.elements(iElements).volume]’;% compute the average gradient by dividing the sum of all phi fluxes for anelement by the volume of the elementfor iElement =1:theMesh.numberOfElementsphiGrad(iElement,:) = phiGrad(iElement,:)/volumes(iElement);end%----------------------------------------------------% Set boundary Gradient equal to associated element GradientphiGrad(iBElements,:) = phiGrad(iOwners_b,:);Listing 7.34 (continued)The loop over the boundary faces in Listing 7.34 did not take into account thepatches to which the boundary faces belong, rather it just looped over all boundaryfaces.
A version that would loop over the various patches and then over theirrespective boundary faces would be written as shown in Listing 7.35.%----------------------------------------------------% BOUNDARY FACES contribution to gradient%----------------------------------------------------% get the list number of boundary Patches%theNumberOfPatches = theMesh.numberOfBoundaries;theEquation = cfdGetModel(theEquationName);%for iPatch=1:theNumberOfPatches%theBoundary = theMesh.boundaries(iPatch);numberOfBFaces = theBoundary.numberOfBFaces;%iFaceStart = theBoundary.startFace;iFaceEnd = iFaceStart+numberOfBFaces-1;iBFaces = iFaceStart:iFaceEnd;%iElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces;iElementEnd = iElementStart+numberOfBFaces-1;iBElements = iElementStart:iElementEnd;iBOwners = [theMesh.faces(iBFaces).iOwner]?;phi_b = phi(iBElements,iComponent);Sb = [theMesh.faces(iBFaces).Sf]';% loop over all boundary faces and add phi flux contribution to gradientfor k= 1: numberOfBFacesphiGrad(iBOwners(k),:) = phiGrad(iBOwners(k),:) + phi_b(k)*Sb(k,:);endendListing 7.35 Boundary faces contribution to gradient implemented by looping over the variouspatches and then over their respective boundary faces7.2 OpenFOAM®1917.2 OpenFOAM®OpenFOAM® [1] uses a finite volume cell-centered discretization of the domainand handles unstructured mesh data format based on the so called face-addressingstorage.
The aim of this data structure is to provide maximum flexibility in thedefinition of unstructured grids in order to allow for the use of a general polyhedronshape.As shown in Fig. 7.7, a polyhedron is a three-dimensional solid consisting of acollection of plane polygons, joined at their edges.
Tetrahedra (polyhedra with fourtriangular faces) and Hexahedra (polyhedra with six faces), for instance, arepolyhedra where the faces are made by polygons of three and four edges, respectively. The possibility to describe any three-dimensional shape and use it as a finitevolume element for the discretization of the equations gives multiple advantagesand flexibility during mesh generation.An efficient way to describe a mesh constructed using polyhedral elements is theface addressing. In face addressing the element shape is of no consequence to theequation discretization process, requiring the use of global addressing whenforming the coefficients.
(Details on the formation of coefficients from the local andglobal perspective will be detailed later on).In OpenFOAM® the data for points, faces, and elements are stored in a numberof lists (arrays), as shown in Fig. 7.8. The list of points contains the threedimensional spatial coordinates, defined as vectors, which correspond to the vertices of the actual mesh. The dimensional unit is the meter.
Moreover each vertexhas a label defined by the position in the list, and because of the use of the C++computer language the label counting starts from zero.The faces are defined by a list of vertex labels referring to points in which theordering is such that each two adjacent points are connected by an edge. The facelist is organized in a way that all internal faces appear first in the list followed byfaces related to the first boundary, then faces related to the second boundary, and soon. Lists of boundary faces are also named patches. It is important to remember thatinternal faces belong to two cells while boundary faces belong to one cell only.Fig.
7.7 A polyhedron element7 The Finite Volume Mesh in OpenFOAM® and uFVM192...points listxyz(coordinates)0123faces list01.n(points indices)012elements list01.k(faces indices)012npnfneFig. 7.8 Points, faces, and elements lists in OpenFOAM®Finally the element or cell list is defined by a list of indices, where the position inthat list is the cell index, the first index at any position is the number of faces for thatelement, and the face indices at a given position represent the faces for that cell.OpenFOAM® can read computational grids from any mesh generator softwarecapable of writing the mandatory files needed by OpenFOAM®.
As explainedearlier, the mesh is named polyMesh and has to be defined with a set of proper files.These files are placed in the directory “constant/polyMesh”. Based on the previousdescription the names of the files whose contents are self-explanatory are given by•••••pointsfacesownerneighbourboundaryIn OpenFOAM® the entire boundary of the domain is described as a list ofpatches in the file boundary.
The file boundary is a list of all defined patchescontaining a proper entry named dictionary with each patch declared using thepatch name, its type, its number of faces (nFaces), and the starting face (startFace).The general syntax used is as written in Listing 7.36.PatchName{type patch;nFaces 'number of faces';startFace 'starting face';}Listing 7.36 Script used to create a patch7.2 OpenFOAM®193An example of a boundary file with two boundary types is shown inListing 7.37.2(inlet{typephysicalTypenFacesstartFace}outlet{typephysicalTypenFacesstartFace}patch;automatic;100;9850;patch;automatic;100;9950;)Listing 7.37 A boundary file with two boundary typesFrom a programming point of view it is worth giving a brief introduction of theC++ classes that handle the mesh and allow access to specific data.The base class that handles the “low-level” structure is called primitiveMesh. Itis a generic class that wraps the geometric information of the mesh withoutassuming any particular form of discretization.
It is the basic class for low levelinformation about the mesh.Examples of some of the functions that are members of the primitiveMesh classare shown in Listing 7.38.constconstconstconstconstconstconstlabelListList &labelListList &cellList &vectorField &vectorField &scalarField &vectorField &cellCells() constpointCells() constcells() constcellCentres() constfaceCentres() constcellVolumes() constfaceAreas() constListing 7.38 Some functions of the primitiveMesh classAs can be seen, this class allows obtaining specific data of the mesh (e.g., cellvolume, cell centers, face centers, etc.).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
