Code of the Week #22: SRFModel rpm

rheoworks · 2019 年2 月 24 日 16:00

#201909

适用版本：foam-extend-3.1

本周要讨论的这个语句位于 tutorials/incompressible/simpleSRFFoam/mixer，是一个算例文件里面的，文件是 constant/SRFProperties

SRFModel  rpm;    // 选择 SRF 子模型

axis (0 0 1);     // 旋转轴向量

rpmCoeffs
{
    rpm  5000.0;  // 转速, r/min
}

在这个算例文件中，完成了旋转机械的稳态计算。计算的结果如下：

计算的是叶片在搅拌器中旋转时的流场。
这里面涉及到比较特殊的两个力是：科里奥力 (Coriolis)，离心力 (Centrifugal force)，处理的方法就是把这两个力加入到动量方程中，即

\nabla\cdot(\boldsymbol{U}_r \boldsymbol{U}_r) + \overbrace{2\boldsymbol{\Omega}\times \boldsymbol{U}_r}^{\text{Coriolis}} + \overbrace{\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times \boldsymbol{r})}^{\text{Centrifugal}} =\nabla\frac{p}{\rho} + \nu\nabla\cdot\nabla \boldsymbol{U}_r\\ \nabla\cdot\boldsymbol{U}_r = 0

\boldsymbol{U}_r=\boldsymbol{U}_i-\boldsymbol{\Omega}\times\boldsymbol{r} 表示相对速度， \boldsymbol{U}_i 表示绝对速度， \boldsymbol{r} 表示半径 (向量)。

从上式可以看到， Coriolis 这项含有未知量 \boldsymbol{U}_r ，因此，这一项可以加入到 U 方程中；而 Centrifugal 这一项不涉及求解量，因此，放入 U 方程的等号右边 (RHS, Right-Hand-Side)，源代码位于 applications/solvers/incompressible/simpleSRFFoam，其中， U 方程的部分如下：

        {
            // Momentum predictor
            tmp<fvVectorMatrix> UrelEqn
            (
                fvm::div(phi, Urel)           // U 方程 LHS 第一项
              + turbulence->divDevReff(Urel)  // RHS 第二项
              + SRF->Su()                     // 科里奥力 + 离心力
            );
...
        }
...
        // Recalculate Uabs
        Uabs = Urel + SRF->U();               // 根据 Urel 计算绝对速度

这里，我们来看一下 SRFModel 中的情况，源代码位于 src/finiteVolume/cfdTools/general/SRF/SRFModel，看一下 SRFModel.C

Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
Foam::SRF::SRFModel::Fcoriolis() const
{
    return tmp<DimensionedField<vector, volMesh> >
    (
        new DimensionedField<vector, volMesh>
        (
            IOobject
            (
                "Fcoriolis",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            2.0*omega_ ^ Urel_ // 计算科里奥力
        )
    );
}


Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
Foam::SRF::SRFModel::Fcentrifugal() const
{
    return tmp<DimensionedField<vector, volMesh> >
    (
        new DimensionedField<vector, volMesh>
        (
            IOobject
            (
                "Fcentrifugal",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            omega_ ^ (omega_ ^ mesh_.C()) // 计算离心力
        )
    );
}


Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
Foam::SRF::SRFModel::Su() const
{
    return Fcoriolis() + Fcentrifugal(); // 返回源项参数
}

SRF->U() 的代码如下：

Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::U() const
{
    return tmp<volVectorField>
    (
        new volVectorField
        (
            IOobject
            (
                "Usrf",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            omega_ ^ (mesh_.C() - axis_*(axis_ & mesh_.C())) // Omega * r
        )
    );
}

该模型的计算需要单独的 Properties 文件，其定义由 SRFModel 的构造函数完成：

Foam::SRF::SRFModel::SRFModel
(
    const word& type,
    const volVectorField& Urel
)
:
    IOdictionary
    (
        IOobject
        (
            "SRFProperties",
            Urel.time().constant(),
            Urel.db(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),
    Urel_(Urel),
    mesh_(Urel_.mesh()),
    axis_(lookup("axis")),
    SRFModelCoeffs_(subDict(type + "Coeffs")),
    omega_(dimensionedVector("omega", dimless/dimTime, vector::zero))
{
    // Normalise the axis
    if (mag(axis_) < SMALL)
    {
        FatalErrorIn
        (
            "SRF::SRFModel::SRFModel\n"
            "(\n"
            "    const word& type,\n"
            "    const volVectorField& Urel\n"
            ")"
        )   << "Zero length axis: " << axis_ << ".  This is not allowed."
            << abort(FatalError);
    }

    axis_ /= mag(axis_);
}

在本算例文件中，选用了 rpm 模型，是 SRFModel 的一个子模型，只是对角速度作了新的转换，按每分钟多少转来设置，更加符合工程实际。其中 rpm 部分的代码 (src/finiteVolume/cfdTools/general/SRF/SRFModel/rpm) 如下：

Foam::SRF::rpm::rpm
(
    const volVectorField& U
)
:
    SRFModel(typeName, U),
    rpm_(readScalar(SRFModelCoeffs_.lookup("rpm")))
{
    // Initialise the angular velocity
    omega_.value() = axis_*rpm_*2.0*mathematicalConstant::pi/60.0;  // 转速 -> 角速度转换
}

在 foam-extend-3.1 版本中， constant/SRFProperties 中的 SRFModel 只能选择 rpm 子模型。

参考文献

Rotating machinery training at OFW11. retrieved on 2018-04-06