Palabos中的单位转换

当我们每次运行palabos算例的时候，tmp文件夹都会附上一个log文件，在这里我们可以检查算例的尺寸，雷诺数等信息。

其源码对应位置在src/core/units.h中，如下：

template<typename T>
void writeLogFile(IncomprFlowParam<T> const& parameters,
                  std::string const& title)
{
    std::string fullName = global::directories().getLogOutDir() + "plbLog.dat";
    plb_ofstream ofile(fullName.c_str());
    ofile << title << "\n\n";
    ofile << "Velocity in lattice units: u=" << parameters.getLatticeU() << "\n";
    ofile << "Reynolds number:           Re=" << parameters.getRe() << "\n";
    ofile << "Lattice resolution:        N=" << parameters.getResolution() << "\n";
    ofile << "Relaxation frequency:      omega=" << parameters.getOmega() << "\n";
    ofile << "Extent of the system:      lx=" << parameters.getLx() << "\n";
    ofile << "Extent of the system:      ly=" << parameters.getLy() << "\n";
    ofile << "Extent of the system:      lz=" << parameters.getLz() << "\n";
    ofile << "Grid spacing deltaX:       dx=" << parameters.getDeltaX() << "\n";
    ofile << "Time step deltaT:          dt=" << parameters.getDeltaT() << "\n";
}

可见，输出的内容是根据parameters提取出来的，那么这个parameters是否还有更多可以调用的功能呢？

首先parameters是从算例cpp文件中定义的。

IncomprFlowParam<T> parameters(
            (T) 1e-2,  // uMax
            (T) 360.,  // Re
            100,       // N
            4.,        // lx
            1.         // ly 
    );

在其类的解释里，第一个latticeU为格子单位的特征速度，第二个为雷诺数，第三个为分辨率（尺寸为1的格子有这么多的cells，这里palabos的cells就是储存分布函数的小格点），后面3个分别为x,y,z方向的无量纲单位长度。

它在不可压流体部分的功能如下，符合网上随意可搜到的单位转换公式：

  IncomprFlowParam(T latticeU_, T Re_, plint resolution_, T lx_, T ly_, T lz_=T() )
    : latticeU(latticeU_), Re(Re_), resolution(resolution_), lx(lx_), ly(ly_), lz(lz_)
    { 
        physicalU      = (T)1;
        physicalLength = (T)1; 
    }
    /// velocity in lattice units (proportional to Mach number)
    T getLatticeU() const { return latticeU; }
    /// velocity in physical units
    T getPhysicalU() const { return physicalU; }
    /// Reynolds number
    T getRe() const      { return Re; }
    /// physical resolution
    T getPhysicalLength() const { return physicalLength; }
    /// resolution
    plint getResolution() const { return resolution; }
    /// x-length in dimensionless units
    T getLx() const      { return getPhysicalLength()*lx; }
    /// y-length in dimensionless units
    T getLy() const      { return getPhysicalLength()*ly; }
    /// z-length in dimensionless units
    T getLz() const      { return getPhysicalLength()*lz; }
    /// lattice spacing in dimensionless units
    T getDeltaX() const  { return (T)getPhysicalLength()/(T)getResolution(); }
    /// time step in dimensionless units
    T getDeltaT() const  { return getDeltaX()*getLatticeU()/getPhysicalU(); }
    /// conversion from dimensionless to lattice units for space coordinate
    plint nCell(T l) const { return (int)(l/getDeltaX()+(T)0.5); }
    /// conversion from dimensionless to lattice units for time coordinate
    plint nStep(T t) const { return (int)(t/getDeltaT()+(T)0.5); }
    /// number of lattice cells in x-direction
    plint getNx(bool offLattice=false) const { return nCell(getLx())+1+(int)offLattice; }
    /// number of lattice cells in y-direction
    plint getNy(bool offLattice=false) const { return nCell(getLy())+1+(int)offLattice; }
    /// number of lattice cells in z-direction
    plint getNz(bool offLattice=false) const { return nCell(getLz())+1+(int)offLattice; }
    /// viscosity in lattice units
    T getLatticeNu() const { return getLatticeU()*(T)getResolution()/Re; }
    /// relaxation time
    T getTau() const       { return (T)3*getLatticeNu()+(T)0.5; }
    /// relaxation frequency
    T getOmega() const     { return (T)1 / getTau(); }

由此可见，你可以自定义输出内容包括某坐标cells数量，松弛时间等参数。

此外，这里还有physicalU和physicalLength两个量，都是1，这是从物理单位转化成无量纲单位再到格子单位用的的参考速度和长度尺度，dt的计算理论上为时间尺度/达到该时长的迭代次数。

在Jonas Latt(2008) Choice of units in lattice Boltzmann simulations第四页上半部分中提到过dt的选择不是很直接，考虑到稳定性因素会与dx相关，这里不过多叙述。

在同文件夹下的units.cpp中，则有更多的单位转换关系，节选如下，其意义一目了然：

double Units3D::physTime(plint lbIter) const {
    return lbIter*dt_;
}

plint Units3D::numCells(double physLength) const {
    return util::roundToInt(physLength/dx_);
}

double Units3D::physLength(plint numCells) const {
    return numCells*dx_;
}

double Units3D::lbVel(double physVel) const {
    return physVel * dt_/dx_;
}

double Units3D::physVel(double lbVel) const {
    return lbVel * dx_/dt_;
}

以下为damBreak3D的重力处理。

// Gravity in lattice units.
T gLB = 9.8 * delta_t * delta_t/delta_x;
externalForce = Array<T,3>(0., 0., -gLB);

下面是unit.cpp里面的关于加速度的描述。

Array<double,3> Units3D::lbAcceleration(Array<double,3> const& physAcceleration) const {
    return physAcceleration * dt_*dt_/dx_;
}

最后我稍微整理了一下Palabos里面的单位转换，并上传至我的github了。
https://github.com/Yulan-Fang/PalabosCodeExplanation
找到Palabos的单位.pdf即可。

–0917/2021更新了第二版。