在OpenFOAM中获取每个网格的Skewness

常用的网格划分软件中关于网格偏度（Skewness）只有统计结果，通常只能看到平均值、最大值及粗略的分布情况，而无法看到每一个网格单元对应的Skewness。因此，这里我希望借助OpenFOAM建立一个生成每个网格单元对应Skewness的小工具。

在OpenFOAM中，可以利用checkMesh工具检查当前网格的质量。checkMesh检测结果包含了网格的Skewness，其输出举例如下：

# ......
Checking geometry...Overall domain bounding box (-0.50095 -0.499013 -1) (0.55 0.499013 1.1)Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)Mesh has 3 solution (non-empty) directions (1 1 1)Boundary openness (4.31927e-16 3.08134e-16 1.46623e-16) OK.Max cell openness = 2.90386e-16 OK.Max aspect ratio = 3.69768 OK.Minimum face area = 0.000292504. Maximum face area = 0.00215511.  Face area magnitudes OK.Min volume = 8.27995e-06. Max volume = 7.0365e-05.  Total volume = 1.14596.  Cell volumes OK.Mesh non-orthogonality Max: 39.9086 average: 4.756Non-orthogonality check OK.Face pyramids OK.Max skewness = 3.7248 OK.Coupled point location match (average 0) OK.Mesh OK.

可以看到，checkMesh工具仅仅能够给出最大的Skewness，而无法给出所有网格单元Skewness的详细信息。但无论如何，我们也可以借助checkMesh工具来研究与Skewness计算相关的代码。

checkMesh工具位于$FOAM_APP/utilities/mesh/manipulation/checkMesh文件夹中，该文件夹下包含诸多文件，但通过内容搜索可以知道，与计算Skewness相关的代码仅在checkGeometry.C文件中。相关的内容如下：

faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
if (mesh.checkFaceSkewness(true, &faces))
{noFailedChecks++;label nFaces = returnReduce(faces.size(), sumOp<label>());if (nFaces > 0){Info << "  <<Writing " << nFaces<< " skew faces to set " << faces.name() << endl;faces.instance() = mesh.pointsInstance();faces.write();}
}

这段代码的含义暂且不细研究。联想到checkMesh工具仅能够给出网格整体的最大Skewness，因此我查了一下网格整体Skewness的计算方法。Fluent用户文档中给出：

Skewness is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume. Highly skewed cells can decrease accuracy and destabilize the solution. For example, optimal quadrilateral meshes will have vertex angles close to 90 degrees, while triangular meshes should preferably have angles of close to 60 degrees and have all angles less than 90 degrees. A general rule is that the maximum skewness for a triangular/tetrahedral mesh in most flows should be kept below 0.95, with an average value that is less than 0.33. A maximum value above 0.95 may lead to convergence difficulties and may require changing the solver controls, such as reducing under-relaxation factors and/or switching to the pressure-based coupled solver.

即Skewness是一个面元与规则面元之间的变形差异情况。此外，Fluent用户文档中也给出了面元的Skewness与网格单元、网格整体的Skewness的计算方式：

Normalized Equiangular Skewness
In the normalized angle deviation method, skewness is defined (in general) as
max[θmax−θe180−θe,θe−θminθe]max\left[\frac{\theta_{max}-\theta_e}{180-\theta_e}, \frac{\theta_e-\theta_{min}}{\theta_e} \right] max[180−θeθmax−θe,θeθe−θmin]
where
θmax\theta_{max}θmax = largest angle in the face or cell
θmin\theta_{min}θmin = smallest angle in the face or cell
θe\theta_{e}θe = angle for an equiangular face/cell (such as 60 for a triangle, 90 for a quad, and so on)
The cell skewness will be the maximum skewness computed for any face. For example, an ideal pyramid (skewness=0) is one in which the 4 triangular faces are equilateral (and equiangular) and the quadrilateral base face is a square.

因此，一个网格单元的Skewness值即为它所包含的面元的最大Skewness值。而网格整体的Skewness值即为所有面元的最大Skewness值。显而易见，面元Skewness较大会影响插值精度，因此在绘制网格时，需要注意生成网格的Skewness范围。

接下来我们就寻找checkFaceSkewness()函数所在的位置。由于前面提到，此函数在checkMesh工具中用于对网格的Skewness进行检查，因此其一定包含计算Skewness的过程。在OpenFOAM的源码文件夹$FOAM_SRC/OpenFOAM中，我们通过Ag/Ack工具搜索checkFaceSkewness字串，可以看到包含此函数的文件主要有两个：

meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

通过比较这两个源文件中checkFaceSkewness()函数所包含的输出信息我们可以知道，checkMesh工具使用的是primitiveMeshCheck.C中所定义的checkFaceSkewness()。到这里，令人惊喜的事情发生了，在checkFaceSkewness()函数中，可以发现它使用了：

tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness(...)

这个函数来计算某个面的Skewness，因此我们继续寻找这个函数，终于在meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C这个文件中看到了faceSkewness()这个函数的定义：

Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
(const primitiveMesh& mesh,const pointField& p,const vectorField& fCtrs,const vectorField& fAreas,const vectorField& cellCtrs
)
{const labelList& own = mesh.faceOwner();const labelList& nei = mesh.faceNeighbour();tmp<scalarField> tskew(new scalarField(mesh.nFaces()));scalarField& skew = tskew.ref();forAll(nei, facei){skew[facei] = faceSkewness(mesh,p,fCtrs,fAreas,facei,cellCtrs[own[facei]],cellCtrs[nei[facei]]);}// Boundary faces: consider them to have only skewness error.// (i.e. treat as if mirror cell on other side)for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++){skew[facei] = boundaryFaceSkewness(mesh,p,fCtrs,fAreas,facei,cellCtrs[own[facei]]);}return tskew;
}

该函数在计算一个非壁面面元的Skewness时，将其与包含此面元的两个网格单元一起分析。输入的数据包括该面元对应顶点的坐标、两个网格单元中心连线向量、面元中心与网格单元中心连线向量、面元的面法相向量。而壁面面元也有自己的特别计算方式。通过该函数，可以输出当前面元的Skewness。

根据以上信息，我所写的计算每个网格单元Skewness的小工具步骤包括：

读取网格信息
使用faceSkewness()函数计算所有面元的Skewness
根据面元与网格单元的对应关系，将每个面元的Skewness映射到对应的网格单元上，并只保留最大的Skewness
输出网格单元上保存的Skewness

源码如下：

#include "IFstream.H"
#include "OFstream.H"
#include "fvCFD.H"
#include "primitiveMesh.H"
#include "primitiveMeshTools.H"
#include "meshSearch.H"int main(int argc, char *argv[])
{#include "setRootCase.H"#include "createTime.H"#include "createMesh.H"const pointField &p = mesh.points();    // All pointsconst vectorField &C = mesh.C();        // All cell centersconst vectorField &fC = mesh.Cf();      // All face centersconst vectorField &fA = mesh.Sf();      // All face area vectors/* const faceList &faces = mesh.faces(); */const labelList &faceOwn = mesh.faceOwner();        // face owner cellconst labelList &faceNei = mesh.faceNeighbour();    // face neighbour cell/*  Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness(const primitiveMesh& mesh,const pointField& p,const vectorField& fCtrs,const vectorField& fAreas,const vectorField& cellCtrs)*/// Calculate the skewness of all faces.scalarField fskewness =primitiveMeshTools::faceSkewness(mesh, p, fC, fA, C);// Initialize the skewness of all cells.scalarField cskewness = mesh.V();forAll(cskewness, cellId) {cskewness[cellId] = 0.0;}// Calculate the skewness of all cells.label faceNeiId = 0;forAll(fskewness, faceId) {cskewness[faceOwn[faceId]]= max(cskewness[faceOwn[faceId]], fskewness[faceId]);if (mesh.isInternalFace(faceId)) {cskewness[faceNei[faceNeiId]]= max(cskewness[faceNei[faceNeiId]], fskewness[faceId]);faceNeiId++;}}// Print.double avg = 0;int num = 0;forAll(cskewness, cellId) {/* Info << cellId << "\t" << cskewness[cellId] << endl; */avg += cskewness[cellId];num++;}avg /= num;Info << avg << endl;return 0;
}

在OpenFOAM中获取每个网格的Skewness相关推荐

openfoam中网格自适应功能dynamicRefineFvMesh
最近看到有些视频介绍openfoam中的网格自适应功能,照着网上的设置,并没有发现网格会被自动加密,原因是使用的openfoam 版本不同,有些操作和案例位置已经发生变化.这个在初学者中很容易引起很多 ...
html中获取modelandview中的json数据_从Bitmap中获取YUV数据的两种方式
从Bitmap中我们能获取到的是RGB颜色分量,当需要获取YUV数据的时候,则需要先提取R,G,B分量的值,然后将RGB转化为YUV(根据具体的YUV的排列格式做相应的Y,U,V分量的排列) 所以这篇 ...
OpenFOAM中slip和noslip介绍（滑移条件无滑移条件）【翻译】
OpenFOAM中slip和noslip介绍(滑移条件&无滑移条件)[翻译] 翻译自:CFD-online 帖子地址:http://www.cfd-online.com/Forums/open ...
5895. 获取单值网格的最小操作数
5895. 获取单值网格的最小操作数给你一支股票价格的数据流.数据流中每一条记录包含一个时间戳和该时间点股票对应的价格 . 不巧的是,由于股票市场内在的波动性,股票价格记录可能不是按时间顺序到 ...
在OpenFOAM中标记某些区域自适应加密
OpenFOAM中自带的自适应加密只支持三维的,这里分两种情况,一种是三维的自适应加密,一种是二维的自适应加密一.三维情况 1.1 更改求解器算例不过这样的加密是三维的,如果是二维的加密看下边 ...
openfoam计算旋转体滑移网格方法和MRF方法（附案例代码）
1 旋转体数值模拟 openfoam中处理旋转体的方法主要包括两种方法:滑移网格法和多重参考系方法,在滑移网格方法中主要采用的为任意网格界面法(AMI),多重参考系方法则是通过设定多重参考系区域,该区 ...
从OVF矢量场文件中获取磁斯格明子的位置和半径的粗略方法（trace skyrmion）
文章目录前言一.使用oommf的avf2odt命令行程序获取斯格明子中心位置的示例二.当磁体系的单个xy平面层仅有一个斯格明子的情况 1.读取所有磁化文件中的指定磁化分量 2.筛选出每一个xy平 ...
知识图：从图和数据库中获取知识
知识图:从图和数据库中获取知识知识图到底是什么,以及关于它们的所有炒作是什么?如果想成为世界各地的Airbnbs,Amazon,Google和LinkedIn,那么学会区分真实的炒作,定义不同类型的 ...
JAVA中获取当前系统时间
JAVA中获取当前系统时间转自:http://www.cnblogs.com/Matrix54/archive/2012/05/01/2478158.html 一. 获取当前系统时间和日期并格式化输 ...

在OpenFOAM中获取每个网格的Skewness

在OpenFOAM中获取每个网格的Skewness相关推荐

最新文章

热门文章