本案例演示利用OpenFOAM编程实现颗粒跟踪并绘制颗粒轨迹。

1 文件结构

利用下面的命令创建文件结构。

run
foamNewApp demo16

文件结构如下图所示。

2 程序源码

本案例只需要一个源文件demo16.C，其内容如下。

 
#include "fvCFD.H"
 
/*
定义一个函数，用于写入包含由点列表定义的单条多段线的vtk文件。
有关vtk文件格式规范，请参阅：
https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
*/
void saveTrajectoryToVtk(DynamicList &path, fvMesh &mesh)
{
    Info << endl
         << "将颗粒路径写入到VTK文件" << endl;
 
    // 创建一个文件夹用于存放VTK文件
    fileName VTK_dir = mesh.time().path() / "VTK";
    mkDir(VTK_dir);
 
    // 创建文件指针
    autoPtr vtkFilePtr;
    vtkFilePtr.reset(new OFstream(VTK_dir / "particle_path.vtk"));
 
    // 写入文件头
    vtkFilePtr() << "# vtk DataFile Version 2.0" << endl
                 << "particle_path" << endl
                 << "ASCII" << endl
                 << "DATASET POLYDATA" << endl;
    vtkFilePtr() << nl;
 
    // 写入点头部
    vtkFilePtr() << "POINTS " << path.size() << " DOUBLE" << endl;
    //  写入点坐标
    forAll(path, ipt)
    {
        vtkFilePtr() << path[ipt].x() << "  " << path[ipt].y() << "  " << path[ipt].z() << endl;
    }
    vtkFilePtr() << nl;
 
    // 写入线的头部与列表
    vtkFilePtr() << "LINES 1 " << path.size() + 1 << endl;
    vtkFilePtr() << path.size() << endl;
    forAll(path, ipt)
    {
        vtkFilePtr() << ipt << endl;
    }
 
    Info << tab << "vtk文件写入完成！" << endl;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
 
    instantList times = runTime.times(); //得到所有的时间
    runTime.setTime(times.last(), 0);    //设置当前时刻为最后一个时间点
    Info << nl << "使用的数据时刻为：" << runTime.timeName() << endl;
 
    // 读取速度场U
    Info << nl << "读取速度场" << endl;
    volVectorField U(
        IOobject(
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE),
        mesh);
 
    // 定义颗粒的初始位置，这里简单起见定义了一个固定坐标
    point particle(0.01, 0.11875, 0);
 
    //准备粒子位置列表和几个标量变量来存储其总行程时间和距离
    DynamicList particlePositions;
    scalar timeTaken(0);
    scalar distanceTravelled(0);
 
    // 添加初始化位置到颗粒位置列表中
    particlePositions.append(particle);
    // 找到无质量颗粒所在的网格，利用网格流速来确定颗粒的下一个位置
    label cellID = mesh.findCell(particle);
    Info << nl << "初始颗粒位于网格 " << cellID << " 中" << endl;
 
    // 描述其他的变量
    vector curPos = particle;
    vector newPos(0, 0, 0);
    label iterCount(1);
    const label maxIters(100);
 
    Info << nl << "开始颗粒跟踪" << endl;
    while (cellID != -1)
    {
        vector velocity = U[cellID];                   //得到网格的速度
        scalar charLen = Foam::cbrt(mesh.V()[cellID]); //得到网格的特征尺寸
        scalar dt = charLen / mag(velocity);           //得到特征时间步长
 
        newPos = curPos + velocity * dt; //计算得到新的颗粒位置
 
        scalar dist = mag(newPos - curPos); //得到颗粒在当前时间步的位移
 
        // 输出信息
        Info << nl << "Lagrangian time step: " << iterCount << nl
             << tab << "current position = " << curPos << nl
             << tab << "new position = " << newPos << nl
             << tab << "local distance travelled = " << dist << nl
             << tab << "local time taken = " << dt << nl
             << tab << "currently reciding cell no. = " << cellID << endl;
 
        // 将数据添加到列表中
        distanceTravelled += dist;
        timeTaken += dt;
        particlePositions.append(newPos);
        curPos = newPos;
        iterCount++;
 
        cellID = mesh.findCell(curPos);
        if (iterCount > maxIters)
        {
            FatalErrorInFunction << "达到最大的跟踪时间步数！" << abort(FatalError);
        }
    }
 
    if (cellID == -1) //颗粒离开计算域
    {
        Info << "颗粒离开计算区域！" << endl
             << "颗粒运行的总距离：" << distanceTravelled << nl
             << "颗粒停留时间：" << timeTaken << endl;
    }
 
    // 写入文件
    saveTrajectoryToVtk(particlePositions, mesh);
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
    Info << nl;
    runTime.printExecutionTime(Info);
 
    Info << "Endn"
         << endl;
 
    return 0;
}
 

3 测试程序

本程序是一个后处理程序，在案例求解完毕后使用。

测试案例计算得到的速度场。

在案例根路径下执行程序demo16，结果如下图所示。

此时在根目录下创建了一个VTK文件夹，其中包含文件particle_path.vtk。

利用paraview读取vtk文件查看轨迹，如下图所示。

（完毕）