Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 32
Текст из файла (страница 32)
In addition, the distance vector Tjoining the centroids of the owner and neighbor elements, the geometric factor gf,the distance vector CN from the owner element centroid to the surface centroid, andthe normal distance to the wall walldist of the owner element centroid are alsostored. Other components will be described later as needed. Moreover, it should benoted that the neighbor for a boundary face is set to −1.As displayed in Listing 7.15 for the element with index 20, which is a boundaryelement, the struct elements contains three lists of indices for neighboring elements(iNeighbours), faces (iFaces), and nodes (iNodes).m.elements(20)ans =index:iNeighbours:iFaces:iNodes:volume:faceSign:numberOfNeighbours:centroid:20[100 103][33 34 1317 1493 1494][168 79 616 705 617 80]3.2484[1 1 1 1 1]2[3x1 double]Listing 7.15 An example of information related to an element7.1 uFVM181The order of the element indices and faces indices are synchronized such thatelements are related to faces in the same order, and boundary faces are defined atthe end of the list of faces.
The struct elements also stores information about theelement centroid and its volume. The faceSign list indicates whether the element isan owner (+1) or neighbor (−1) for the respective faces.Elements can be identified on the mesh as in Fig. 7.3 using the followingcommand (Listing 7.16):cfdPlotElements([20 300])Listing 7.16 Script needed for uFVM to display selected elements [here elements (20) and (300)]highlighted on the meshFig. 7.3 A two dimensionalview of the mesh over anelbow highlighting elements(20) and (300)Information about element (20) were displayed above while the attributes ofelement (300), which has two boundary faces, are as shown in Listing 7.17.7 The Finite Volume Mesh in OpenFOAM® and uFVM182m.elements(300)ans =index:iNeighbours:iFaces:iNodes:volume:faceSign:numberOfNeighbours:centroid:300[278 302 590][407 435 436 2053 2054][283 820 679 142 290 827]1.9083[-1 1 1 1 1]3[3x1 double]Listing 7.17 An example of information related to element (300) with two boundary facesFinally information about the various boundary patches, i.e., the list of theboundary faces is stored in struct boundaries.
Each boundary array containsinformation about the starting index of the first boundary face, the number ofboundary faces belonging to the patch, in addition to the physical type of the patchand its name. For the example considered, information about boundary patch {1} isshown in Listing 7.18.>> m. boundaries(1)ans =userName:index:type:numberOfBFaces:startFace:'wall-4'1'wall'1001301Listing 7.18 An example of information stored for a wall boundaryMoreover, information about the first boundary face, for example, is obtained asshown in Listing 7.19.>> m.faces(1301)ans =iNodes:index:iOwner:iNeighbour:centroid:Sf:area:CN:geoDiff:T:gf:walldist:iOwnerNeighbourCoef:iNeighbourOwnerCoef:[38 53 590 575]13011-1[3x1 double][3x1 double]3.7510[3x1 double]5.6264[3x1 double]10.6667[][]Listing 7.19 Information about a boundary face7.1 uFVM183To be noted is the index of the first boundary face, which is 1300+1 since inMatlab® arrays start at index 1 while in the C computer language arrays start at index 0.Thus to loop over the faces of a certain patch, the index of the starting face andthe number of faces from the struct boundary associated with the said patch areneeded.
For example to loop over the faces of patch 2, the following script(Listing 7.20) can be used:theMesh = cfdGetMesh;iPatch = 2;iBFaces = cfdGetFaceIndicesForBoundaryIndex(iPatch)for iBFace=iBFacestheBFace = theMesh.faces(iBFace);disp(theBFace) %display theBFace internal fieldsendListing 7.20 Looping over boundary patch facescfdGetFaceIndicesForBoundaryIndex is defined in Listing 7.21 as follows:theIndices = cfdGetFaceIndicesForBoundaryIndex(theBoundaryIndex)%theBoundary = cfdGetBoundary(theBoundaryIndex);theNumberOfBFaces = theBoundary.numberOfBFaces;theStartFace = theBoundary.startFace;theIndices = [theStartFace:theStartFace+theNumberOfBFaces-1];%endListing 7.21 Indices of boundary faces for a specific patch7.1.4 uFVM Geometric FieldsIn addition to all data stored in the struct mesh, information about the model to besolved and the values of the fields of interest should be stored to be accessed whenneeded.
Three types of locale can be identified and fields of different types can bedefined on them. Namely for locale (Elements, Faces, and Nodes) and for types(scalars, vectors, and tensors). The various fields are first defined based on theirlocale.7.1.4.1 The Element FieldsAn element field is constructed using the following script (Listing 7.22):cfdSetupMeshField(theUserName,theLocale,theType,theTimeStep)Listing 7.22 Script needed to construct an element field7 The Finite Volume Mesh in OpenFOAM® and uFVM184where theUserName is the name of the field, theLocale is the geometric entity overwhich it is defined (Elements, Faces, Nodes), theType defines the type of the arrayelements (Scalar or Vector), and finally the TimeStep indicates the relative time ofthe field (Step0, Step1, etc.).
For example, the following script (Listing 7.23):>> UField = cfdSetupMeshField('U:water','Elements','Vector','Step0')UField =userName: 'U:water'name: 'U_fluid01'type: 'Vector'locale: 'Elements'phi: [2908x3 double]Listing 7.23 Example of setting up a vector fieldsets up a vector field defined over Elements at time step 0, i.e., at the current timestep.As shown in Fig. 7.4, the array has the size of the NumberOfElements+theNumberOfBoundaryFaces since the element array will include the valuesdefined for each element in the mesh in addition to each boundary face. Theboundary face values represent boundary conditions for that field.
These boundaryvalues are grouped in terms of boundary patches.elements array01 2boundary faces...NE...patch#1(np1 faces)...patch#2(np2 faces)NFpatch#nFig. 7.4 Size of the field arrayGenerally they can be accessed as follows: For example to initialize theboundary values of the UField at patch 1 to a value [1 0 0], the following should bewritten (Listing 7.24):7.1 uFVM185% get the meshtheMesh = cfdGetMesh;% get information about the boundary patchtheBoundary = theMesh.boundaries(iPatch);numberOfElements = theMesh.numberOfElements;numberOfInteriorFaces = theMesh.numberOfInteriorFaces;numberOfBFaces = theBoundary.numberOfBFaces;% Starting faceiFaceStart = theBoundary.startFace;% get information about starting and ending elementsiElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces;iElementEnd = iElementStart+numberOfBFaces-1;% define the indices as an index arrayiBElements = iElementStart:iElementEnd;>> UField.phi(iBElements,:) =cfdComputeFormulaAtLocale('[1;0;0]','BPatch1','Vector')ans =100100100Listing 7.24 The script needed to initialize the boundary value of a field in a patchwhere the statement in Listing 7.25 given bycfdComputeFormulaAtLocale(theFormula,theLocale,theType)Listing 7.25 Statement used to compute values over a localeevaluates the values in theFormula over theLocale and returns and array of theappropriate length of type theType (Scalar or Vector).7.1.4.2 The Face FieldsA face field is constructed with theLocale set to ‘Faces’ as (Listing 7.26)cfdSetupMeshField(theUserName,theLocale,theType,theTimeStep)Listing 7.26 Statement used to construct a face fieldAgain, as shown in Fig.
7.5, the length of the array is equal to numberOfFaces,which combines the numberOfInteriorFaces plus the sum of all boundary faces.7 The Finite Volume Mesh in OpenFOAM® and uFVM186faces arrayinterior faces01 2boundary faces...NIF......patch#1(np1 faces)startFacefor patch#1patch#2(np2 faces)startFacefor patch#2NFpatch#n...Fig. 7.5 Size of faces arrayThe boundary faces for any boundary patch can be accessed as (Listing 7.27)theMesh = cfdGetMesh;numberOfElements = theMesh.numberOfElements;numberOfInteriorFaces = theMesh.numberOfInteriorFaces;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;thBFaces = theMesh.faces(iBFaces)Listing 7.27 Script to access the various arrays of struct facesand the boundary elements for a scalar field can then be retrieved using Listing 7.28 asphi_b = phi[iBElements]Listing 7.28 Accessing the boundary elements of a scalar phi7.1.4.3 The Node FieldThe node field is a list where the indices of vertices of faces are stored.
Each face isreferred to by the indices of its vertices in the points list (Fig. 7.6). The position of7.1 uFVM187the face in the list refers to the face index. Therefore face 0 is the first entry in thelist, face 1 is the second entry in the list, and so on.points array01 2...NPFig. 7.6 A points array7.1.5 Working with the uFVM MeshLooping over elements, interior faces, boundary faces, boundary elements, or evenboundary patches are common operations that are performed during the discretization and solution cycles.
It is worth reviewing the mechanisms that allow performing these operations readily.7.1.5.1 Looping Over ElementsThis is a simple loop to implement as the number of elements is already known and theelements are indexed from 1 to the numberOfElements (0 to numberOfElements-1for OpenFOAM®). Also element fields are indexed in a similar fashion, so accessingthem is as simple as accessing elements. Thus the loop is simply as shown inListing 7.29.for iElement=1:numberOfElementstheElement = theMesh.elements(iElement)phi(iElement) % this is field phi at centroid of element iElement....endListing 7.29 Script used to loop over elementsTo access the boundary elements, the script shown in Listing 7.30 is used.for boundary = 1:numberOfBoundariestheBoundary = theMesh.boundaries{1};endListing 7.30 Script used to access boundary elements1887 The Finite Volume Mesh in OpenFOAM® and uFVM7.1.5.2 Looping Over FacesThe faces array is constructed so that all interior faces run through indices 1 tonumberOfInteriorFaces followed by the boundary faces that are also arrangedaccording to the boundary patch to which they belong.
Thus looping over theinterior faces is straight forward and is written as in Listing 7.31.for iFace=1:numberOfInteriorFacestheFace = theMesh.faces(iFace)endListing 7.31 Script used to loop over interior facesLooping over the boundary faces can be done as shown in Listing 7.32.for iBFace= numberOfInteriorFaces+1:numberOfFacestheBFace = theMesh.faces(iBFace)endListing 7.32 Script used to loop over boundary facesIf the interest is to loop over the boundary faces of a certain patch, then the loopruns from a start face defined in the boundary condition to nFaces also defined inthe boundary condition. Thus a loop over the faces of boundary patch n would bewritten as depicted in Listing 7.33.startFace = theMesh.boundaries(n).startFacenFaces = theMesh.boundaries(n).nFacesfor iBFace = startFace: startFace+nFaces-1Listing 7.33 Script used to loop over boundary faces of a certain patchThe subroutine cfdGetFaceIndicesForBoundaryIndex is used to directly getthe vector startFace: startFace+nFaces-1.7.1.6 Computing the Gauss GradientComputing the gradient of an element field in uFVM, necessitates the use of manyof the above routines.