上一节概述了OpenFOAM中边界条件的实现。 在这一节中,描述了两个新的边界条件的植入。

OpenFOAM提供了大量不同的边界条件可供选择,并且在代码的这一部分也进行了社区开发,其中最突出的是swak4Foam贡献的groovyBC边界条件。在编写符合自身需求的新边界条件之前,明智的做法是查看代码库中是否已经提供了此功能,或者是否可以通过groovyBC边界条件进行建模。

关于如何在OpenFOAM中编写自己的边界条件,互联网上已经有相当多的信息。 本章中的例子是独立于现有材料编写的。 第一个例子展示了如何扩展OpenFOAM中任何边界条件的功能而不修改它,具体目的是减少边界处的再循环。第二个示例演示如何开发一个新的pointPatchField来将预定义的运动应用到patch。 这个预定义的运动是从必须由用户提供的表格数据中计算出来的。 这两个例子都强调了OpenFOAM中边界条件的根抽象基类所固定的类界面的正确使用,即fvPatchField和pointPatchField类。

10.3.1再循环控制边界条件

在本例中,介绍了如何在运行期间通过附加计算(功能)扩展现有边界条件的方法。假设有一个边界条件,在仿真执行时以某种方式更新边界值。在某个点,基于模拟结果,模拟的条件(例如,在另一边界处的压力)发信号通知附加的边界操作是必要的。发生这种情况时,新的延伸边界条件会在执行时间启用额外的计算,并相应地修改边界值。

一个很好的技术例子是充满理想气体的加热密闭容器。当热通量进入容器时,容器内的压力上升。扩展边界条件将测量容器盖处的压力,并且当压力达到特定值时打开盖。从数值上讲,该边界条件将在模拟过程中根据盖处的压力改变其类型,从封闭壁面改变为出口边界条件。

在本章中给出的例子中,再循环是在边界处测量的。 当再循环发生时,附加计算将边界条件的性质修改为流入条件。

注:再循环边界条件示例在并行执行中可能不起作用,这取决于修改的边界是否是处理器域的一部分。本例的目的是展示网格、对象注册表和物理场如何相互连接,而不是如何并行化边界条件更新。

假设只使用继承就可以实现这样的扩展。 然而,这将需要在OpenFOAM中扩展每个边界条件,以说明需要执行的启用RTS的额外计算。 通过使用多重继承来扩展每个边界条件将导致对现有边界条件的修改,只是为了考虑可能的扩展(即循环控制),这取决于用户的选择,可以在运行时使用,也可以不使用。 显然,这绝不是一个通用的方法,它使得在运行时不修改现有代码就很难实现扩展。 当一个新功能需要在运行时添加到一个已经存在的类层次结构中而不修改该层次结构的模型时,可以使用面向对象的设计模式装饰器模式。 关于decorator模式和其他OOD模式的细节在书中由[1]给出。 图10.4阐明了零次边界条件的工作原理,该边界条件是用任意函数修饰的。

图10.4 装饰边界条件的工作原理。 边界条件以标准的方式运行,直到满足这样的条件,即需要扩展计算并将其委托给装饰器。在这一点上,边界条件的行为就像它切换到了装饰器类型一样。

10.3.1.1 使用decorator模式向BC添加功能

边界条件修饰器本身就是一个边界条件,因为它也修改边界场。 因此,装饰器继承自OpenFOAM中所有边界条件的fvPatchField抽象基类。 除了从fvPatchField继承之外,装饰器还组成它自己的基类(fvPatchField)的对象。

为了澄清这一点,图10.5中给出了应用于边界条件层次结构的装饰器模式的统一建模语言(UML)类协作图。 如图所示,装饰器和其他边界条件一样是类层次结构的一部分。 因此,对于OpenFOAM中的边界条件,使用抽象基类fvPatchField施加这样的IS-A关系可以使边界条件修饰器充当OpenFOAM客户端代码其余部分的边界条件。 这与普通继承的原理完全相同,只是扩展了存储继承的对象的实例。

图10.5 OpenFOAM边界条件的修饰器

装饰器必须实现抽象类FvPatchField规定的所有纯虚方法。 由于它也合成了一个普通的边界条件,装饰器将根据程序员规定的条件将函数调用委托给装饰的边界条件。 decorator提供的扩展可以根据程序员的需要进行设计。 如图10.5所示,通过结合继承和组合,修饰的边界条件可以在运行时切换自己的类型,因为它将计算委托给修饰器。

10.3.1.2 在边界条件中加入再循环控制

作为在OpenFOAM中为任何边界条件添加运行时功能的一个例子,我们选择流动再循环作为控制参数,来流速度作为边界条件的施加作用。 流动再循环由边界处的体积流量交替符号来识别。 区域边界上体积通量中的交替符号表示涡旋(涡)越过边界。 因此,流体通过边界的一部分流入区域,并通过边界的另一部分流出区域。

在算例模拟中,边界条件检查扩展边界的流动是否存在回流,并试图控制修饰边界条件以减少回流。 这种类型的边界条件修饰器最有可能被设置为一些预期流出的边界条件。 本例中的控件是以一种相当简单的方式完成的:通过修改修饰的边界条件,使其场被增加的流入值覆盖。

当然,这个示例是特定于流入/流出情况的,但是这个示例的目标不是处理流控制。 再循环控制边界条件修饰器不同于OpenFOAM标准边界条件 它直接修改由另一个边界条件计算的场值,而不考虑被修饰的边界条件的类型。 示例代码存储库中已经提供了完成的再循环控制边界条件。 为了便于理解所给出的示例,在描述之后应该使用在文本编辑器中打开的示例代码库中的再循环控制边界条件的代码。 在下面的讨论中,文件和类的名称与示例代码存储库中的名称相同。

注:边界条件必须编译为共享库。 它们永远不应该直接实现到应用程序代码中,因为这极大地限制了它们的可用性以及与其他OpenFOAM程序员的共享。

边界条件编译成的库在运行时与求解器应用程序动态链接。 在本教程开始时,需要创建库目录:

?>  mkdir -p primerBoundaryConditions/recirculationControl
?>  cd !$

并且需要复制一个现有边界条件的类文件,我们将其用作再循环控制边界条件的骨架文件:

?>  cp \
$FOAM_SRC/finiteVolume/fields/fvPatchFields/basic/zeroGradient/* .

字符串“Zerogradient”的所有实例都将在文件名和类名中重命名为“RecirculationControl”。 修改名称后,退出RecirculationControl目录,并创建编译配置文件夹make:

?>  cd ..
?>  $FOAM_SRC/../wmake/scripts/wmakeFilesAndOptions

并修改make/files以考虑到要编译的是一个库,而不是一个应用程序,因此行

EXE = $(FOAM_APPBIN)/primerBoundaryConditions

需要修改为:

LIB = $(FOAM_USER_LIBBIN)/libofPrimerBoundaryConditions

请注意,脚本wmakeFilesAndOptions将把所有*.C文件插入make/files文件中,这意味着行:

recirculationControl/recirculationControlFvPatchField.C

必须在编译前移除,因为它限制了循环控制边界条件的类模板定义。 这将在编译期间导致类重新定义错误,因为实际的边界条件类是使用宏从张量属性的循环控制类模板中实例化的。

makePatchFields(recirculationControl);

一旦从make/files中删除了文件recirculationControlFvPatchField.C,库就可以进行第一次编译测试了。 编译可以通过命令

?>  wmake

从primerBoundaryConditions目录中。 清单49与清单50分别为make/files及make/options文件的最终版本。

# Make/files文件内容
recirculationControl/recirculationControlFvPatchFields.C
LIB = $(FOAM_USER_LIBBIN)/libofPrimerBoundaryConditions
# Make/options文件内容
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

编译成功后,创建了一个自包含的共享边界条件库的框架实现。 在Decorator功能的进一步开发开始之前,要用实际的模拟运行来测试边界条件。 在这一点上,边界条件与它所基于的zeroGradientFvPatchField具有相同的功能,但名称不同。 它的功能可以通过运行任何模拟案例来测试,只要您在system/controlDict文件中定义所需的链接库:

libs ("libofPrimerBoundaryConditions.so")

这将加载边界条件共享库,并将其链接到OpenFOAM可执行文件。

注:在编写自定义库的过程中,强烈建议执行此类集成测试。使用版本控制系统(第6章)也可以改进工作流。例如,尝试使用icoFoam求解器将再循环控制边界条件应用于空腔教程案例中速度场的fixedWalls边界。

此时,库与OpenFOAM的集成已经过测试,并被认为是有效的。图10.5所示的设计的实施可以开始了。实现Decorator模式的严格的面向对象方法是从fvPatchField的抽象Decorator类开始。在这种情况下,recirculationControl必须作为具体的装饰器模型从它派生出来。相反,在这个例子中使用了一个具体的装饰器作为起点,抽象的实现留给读者作为本节末尾的练习。通过实施该练习,可以在OpenFOAM中为任何边界条件添加不同的功能,而无需修改它们的实现。除了增进对OpenFOAM中边界条件的了解之外,这是进行练习的另一个好处。

对再循环控制边界条件的第一个修改应用于类模板声明文件recirculationControlFvPatchField. H。它继承自fvPatchField,fvPatchField已经由zeroGradientFvPatchField完成。recirculationControlFvPatchField类定义的重要部分显示在代码清单中。

在清单51中,可以看到装饰器模式:存在从fvPatchField的继承,但fvPatchField的私有属性也由再循环控制边界条件类合成。在baseTypeTmp_中实例化本地副本是使用OpenFOAM的智能指针对象tmp完成的,因为它提供了额外的功能,如垃圾收集。这种智能指针依赖于RAII C++习惯用法,极大地简化了指针的处理。声明为常量的其余类属性将从边界条件的字典中读取。

// 程序清单51
template<class Type> 
class recirculationControlFvPatchField
:
    public fvPatchField<Type> 
    {
        protected:
            // Base boundary condition.
            tmp<fvPatchField<Type> >  baseTypeTmp_;
            const word applyControl_;
            const word baseTypeName_;
            const word fluxFieldName_;
            const word controlledPatchName_;
            const Type maxValue_;
            scalar recirculationRate_;

为了实例化应修饰的边界条件,必须通过边界条件的字典传递特定名称,并将其存储在baseTypeName_ attribute中。为了提供应用或不应用再循环控制的选项,实现了开关applyControl_。如果它为false,则不会在边界上施加任何内容,因此关闭它并简单地使用修饰的边界条件是很容易的。此外,它还将报告扩展边界条件下的再循环量。延伸边界条件的类型由baseTypeName_变数决定,然后使用该变数根据类别名称建立具体的延伸边界条件。此边界条件存储在baseTypeTmp_中。

为了计算再循环,边界条件需要知道体积通量场的名称,该名称是fluxFieldName_属性定义的。受控修补程序字段的名称由controlledPatchName_定义。它可能采用的最大值由成员maxValue_定义,recirculationRate_存储扩展边界条件上负体积通量的百分比。

方法是扩展和修改构造函数初始化列表,以考虑新的私有类属性。为了保持描述简洁,以下文本中仅显示字典构造函数。源代码库中提供了完整的实现。

再循环控制边界条件的构造函数如清单52所示。在初始化列表中,baseTypeTmp_未被赋值,因此取任意值。在建构函式本身中,会使用New选取器建立fvPatchField并指派给baseTypeTmp_。这种构造对象的方式称为工厂方法(新选择器),并在抽象fvPatchField类中定义。它在运行时根据字典中提供的名称选择recirculationControlFvPatchField的基类,构造函数使用字典来初始化修饰边界条件。在执行此构造函数后,将初始化控制函数所需的受保护属性以及修饰边界条件。

template<class Type> 
Foam::recirculationControlFvPatchField<Type> ::
recirculationControlFvPatchField
(
    const recirculationControlFvPatchField<Type> & ptf,
    const fvPatch& p,
    const DimensionedField<Type, volMesh> & iF,
    const fvPatchFieldMapper& mapper
)
:
    fvPatchField<Type> (ptf, p, iF, mapper),
    baseTypeTmp_(),
    applyControl_(ptf.applyControl_),
    baseTypeName_(ptf.baseTypeName_),
    fluxFieldName_(ptf.fluxFieldName_),
    controlledPatchName_(ptf.controlledPatchName_),
    maxValue_(ptf.maxValue_),
    recirculationRate_(ptf.recirculationRate_)
{
    // Instantiate the baseType based on the dictionary entries.
    baseTypeTmp_ = fvPatchField<Type> ::New
    (
        ptf.baseTypeTmp_,
        p,
        iF,
        mapper
    );
}

recirculationControlFvPatchField的其余构造函数必须以类似方式进行修改,以说明新的类属性。一旦修改了构造函数,图10.5中列出的成员函数也需要修改。它们需要将工作委托给存储在baseTypeTmp_中的修饰边界条件。所有四个成员函数的实现都是相同的,所以清单53中只给出了其中一个。

// 程序清单53
template<class Type> 
Foam::tmp<Foam::Field<Type> > 
Foam::recirculationControlFvPatchField<Type> ::valueInternalCoeffs
(
    const Foam::tmp<Foam::Field<scalar> >  & f
) const
{
    return baseTypeTmp_-> valueInternalCoeffs(f);
}

这种委托需要出现在装饰器实现的所有部分中,在这些部分中不执行流的控制。在这种情况下,边界条件应作为装饰边界条件。图10.5中的具体边界条件模型对此进行了说明。

作为实现再循环控制边界条件的最后一步,需要实现负责边界条件运算的成员函数(updateCoeffs)。实现分为两个主要部分。第一部分检查边界条件是否已更新,这是OpenFOAM边界条件中的常见做法。如果不是这样,则使用通量场fluxFieldName_的用户定义名称查找通量场。

正、负体积通量由扩展边界条件计算,如清单54所示。一旦计算出正和负体积流量,再循环率(负流量与总流量之比)将使用清单55中所示的代码进行计算。

// 程序清单54
if (this-> updated())
{
    return;
}
typedef GeometricField <Type, fvPatchField, volMesh>  VolumetricField;
// Get the flux field
const Field<scalar> & phip =
this-> patch().template lookupPatchField
<
    surfaceScalarField,
    scalar
> (fluxFieldName_);
// Compute the total and the negative volumetric flux.
scalar totalFlux = 0;
scalar negativeFlux = 0;
forAll (phip, I)
{
    totalFlux += mag(phip[I]);
    if (phip[I] < 0)
    {
        negativeFlux += mag(phip[I]);
    }
}
// 程序清单55
// Compute recirculation rate.
scalar newRecirculationRate = min
(
    1,
    negativeFlux / (totalFlux + SMALL)
);

Info << "Total flux " << totalFlux << endl;
Info << "Recirculation flux " << negativeFlux << endl;
Info << "Recirculation ratio " << newRecirculationRate << endl;
// If there is no recirculation.
if (negativeFlux < SMALL)
{
    // Update the decorated boundary condition.
    baseTypeTmp_-> updateCoeffs();
    // Mark the BC updated.
    fvPatchField<Type> ::updateCoeffs();
    return;
}

在不发生再循环的条件下,再循环控制和装饰边界条件的行为相同。因此,场的修改被委托给封装的混凝土边界条件模型。注意,因为边界条件装饰器本身是边界条件,所以每次发生更新时,边界条件状态必须被设置为最新。通常,必须为此调用fvPatchField抽象类的成员函数updateCoeffs。updateCoeffs成员函数中负责控制再循环的部分如清单56所示。

// 程序清单56
if (
    (applyControl_ == "yes") &&
    (newRecirculationRate >  recirculationRate_)
)
{
    Info << "Executing control..." << endl;
    // Get the name of the internal field.
    const word volFieldName = this-> internalField().name();
    // Get access to the regitstry.
    const objectRegistry& db = this-> db();
    // Find the GeometricField in the registry using
    // the internal field name.
    const VolumetricField& vfConst =
    db.lookupObject<VolumetricField> (volFieldName);
    // Cast away constness to be able to control
    // other boundary patch fields.
    VolumetricField& vf = const_cast<VolumetricField&> (vfConst);
    // Get the non-const reference to the boundary
    // field of the GeometricField.
    typename VolumetricField::Boundary& bf =
    vf.boundaryFieldRef();
    // Find the controlled boundary patch field using
    // the name defined by the user.
    forAll (bf, patchI)
    {
        // Control the boundary patch field using the recirculation rate.
        const fvPatch& p = bf[patchI].patch();
        if (p.name() == controlledPatchName_)
        {
            if (! bf[patchI].updated())
            {
                // Envoke a standard update first to avoid the field
                // being later overwritten.
                bf[patchI].updateCoeffs();
            }
            // Compute new boundary field values.
            Field<Type>  newValues (bf[patchI]);
            scalar maxNewValue = mag(max(newValues));
            if (maxNewValue < SMALL)
            {
                bf[patchI] == 0.1 * maxValue_;
            } else if (maxNewValue < mag(maxValue_))
            {
                // Impose control on the controlled inlet patch field.
                bf[patchI] == newValues * 1.01;
            }
        }
    }
}

该边界条件将需要对另一边界条件的场进行非常量的访问,因为它将修改边界场值。因此,它通过强制丢弃注册表提供的VolumetricField的常量来破坏objectRegistry类的封装。清单57显示了对另一个边界字段的非常量访问。

//程序清单57
// Cast away constness to be able to control other boundary patch fields.
VolumetricField& vf = const_cast<VolumetricField&> (vfConst);
// Get the non-const reference to the boundary field of the GeometricField.
typename VolumetricField::GeometricBoundaryField& bf = vf.boundaryField();

注:在任何可能的情况下,都应该避免以这种方式抛弃constness:它使封装点无效--只能由类成员函数修改的对象状态。此示例使用常量转换仅指向几何字段和对象注册表之间可能的协作。

如果objectRegistry的接口只提供对已注册对象的常量访问,那么抛弃常量并修改对象状态就是在欺骗程序员。他/她将不期望经由objectRegistry的接口能够改变VolumetricField对象状态。算法1阐明了伪代码中的再循环控制算法。

该边界条件主要是作为将装饰器模式应用于OpenFOAM中的边界条件类层次结构的示例。然而,这里应用的设计可以很好地用于存在入口/出口边界条件的情况。在某些情况下,可以通过增加例如入口压力或速度来实现再循环的减少。不过,请注意,这只是说明两个主要问题的示例:首先,它显示了如何使用VolumetricField的设计来耦合不同边界条件之间的功能。其次,说明了如何在OpenFOAM中编程一个新的边界条件,该边界条件是从fvPatchField导出的。

10.3.1.3 测试再循环控制边界条件

实现updateCoeffs方法后,即可使用边界条件。作为一个测试用例,使用了一个带有后向台阶的简单后向通道。在通道出口采用回流控制,控制边界条件为后向台阶壁。

图10.6 再循环控制测试用例的几何设置和初始条件

图10.6描述了流动的初始形态,其中向后台阶的速度设置为零,因为它是一个不可渗透的壁面。当出口出现再循环时,再循环控制边界条件将超越后向台阶的零速度。进一步,它将把应用于左壁面边界的边界条件转换为入流,以便将再循环排出。

图10.7给出了有无再循环控制边界条件下速度场的比较。

图10.7 再循环控制测试用例的速度场,无和有再循环控制,分别有零梯度和再循环出口

注:此测试用例的设置可在示例用例库的ofbookcases/chapter 10/recirculationControlChannel文件夹中获得,并使用icoFoam求解器对瞬态层流不可压缩单相流进行模拟。

练习:

此示例中未实现抽象装饰器:修改recirculationControlFvPatchField,以便将抽象装饰器添加到类层次结构中。该装饰器概括了边界条件的装饰,从而能够在运行时向OpenFOAM中的任何边界条件添加任何功能。

10.3.2 网格运动边界条件

本小节说明了用于网格运动的新边界条件的构造。网格的运动依赖于以边界条件的形式为网格的点定义的位移或速度,该边界条件与属于网格边界的边界点相关联,在OpenFOAM中称为pointPatchField。与基于fvPatchField的边界条件不同,pointPatchField类型的边界条件不会将边界值存储在边界字段中,而是用于修改内部字段的值。更改内部字段值和任何其他操作由求解器或使用这些边界条件的其他类处理。网格运动边界条件所使用的矢量量将定义特定网格点的速度或位移,这取决于网格运动求解器的选择。

在本节中,网格运动解算器将使用dynamicFvMesh库,在模拟示例中,它将使用新边界条件定义的边界处的位移解算点位移的拉普拉斯方程。由于拉普拉斯方程用于模拟扩散输运,在网格边界处规定的位移将平滑地扩散到周围的网格,这确保了变形单元的更高质量。对于此类应用程序,存储点变形的字段称为pointDisplacement。

本节中介绍的边界条件从输入文件中读取面片重心的位置和方向,并将位移应用于网格边界(相对于其先前位置),并且必须与动态网格结合使用,否则将无法读取pointDisplacement字段。边界条件功能由OpenFOAM的两个现有组件组成:

  1. 计算面片重心(COG)的位置,该位置是基于字典中包含的指定运动计算的。这种规定的运动必须以表格的形式出现,并在每个数据点之间进行线性插值。所有这些都已在tabulated6DoFMotion中实现,它是一个dynamicFvMesh类,可根据指定的运动移动整个网格。
  2. 将向量值赋给pointPatchField。 这类边界条件的一个例子是oscillatingDisplacementPointPatchVectorField

对于从fvPatchField导出的边界条件,仅更改域边界上的值。任何边界条件都不会对模拟或现场进行额外更改,并且功能以逻辑方式封装。对于fvPatchField边界条件,流求解器使用边界值计算场变量。相同的原则适用于从pointPatchField导出的边界条件。边界条件仅规定速度或位移,而实际网格更改由专用网格运动解算器(dynamicFvMesh库的一部分)执行。在下文中,将简要介绍这两个组件,并突出显示与新边界条件相关的部分。

10.3.2.1 读取运动数据

在OpenFOAM的现有代码库中进行简单搜索,发现有一个网格运动解算器,它从表格文件中读取运动数据,与此边界条件的计划类似。但是,在这种情况下,会将相同的位移套用至所有网网格节点,从而产生将网面做为刚体移动的网面运动:不会发生网格变形,网格点的相对位置不会改变。出于本教程的目的,可以使用运动计算,然后将其应用于面片点,使网格边界作为刚体移动。当物体相对于网格(流域)的相对运动小时,这种网格运动可能是有益的。在这种情况下,如果运动以类似扩散的方式传播到流域中,则远离具有规定位移的边界的网格的运动将接近于零。位移从物体消失的速度由位移扩散系数的大小决定。

基于表格数据实现刚体网格运动的代码包含在tabulated6DoFMotionFvMesh类中,该类派生自solidBodyMotionFunction,可在以下位置找到:

$FOAM_SRC/dynamicMesh/motionSolvers/displacement/solidBody/\
solidBodyMotionFunctions/tabulated6DoFMotion/

在正式版本中有一个示例案例,它使用tabulated6DoFMotion来规定封闭罐的运动。本教程可在此处找到:

$FOAM_TUTORIALS/multiphase/interFoam/laminar/sloshingTank3D6DoF

包含运动数据(时间点)的字典被建立为列表(使用列表数据结构),该列表由时间t中的平移和旋转向量组成。此数据的示例可在上述教程中的constant/6DoF.dat中找到。下面的代码片段说明了字典建立的原理。

(
    (t1 ((Translation_Vector_1) (Rot_Vector_1))
    (t2 ((Translation_Vector_2) (Rot_Vector_2))
    ...
)

执行基于样条的插值以获得字典的数据点之间的位置和取向数据。平移和旋转向量是相对于原始坐标系定义的,如图10.8所示。

图10.8:相对于全局坐标系从t0位置移动到tn位置的示例面片的图示

由于第13章仅涉及动态网格,因此为了更清楚起见,仅提供了solidBodyMotionFvMesh类型的动态网格的工作原理的简要总结:

  • 从solidBodyMotionFvMesh派生的动态网格仅处理实体的运动
  • 不能执行拓扑更改,实体面片也不能更改其形状
  • 运动本身不是由solidBodyMotionFvMesh定义的,这是solidBodyMotionFunction和派生类的作用
  • 这简化了运动计算与实际网格运动算法的分离
  • tabulatedSixDoF类继承自solidBodyMotionFunction,后者是OpenFOAM中所有实体运动函数的基类。

如果在dynamicMeshDict中选择了类型为solidBodyMotionFvMesh的动态网格,则在使用运行时选择构建dynamicFvMesh期间将实例化solidBodyMotionFvMesh。从常量/dynamicMeshDict的子字典中读取特定参数,常量/dynamicMeshDict是称为SBMFCoeffs_的私有类属性。下文描述了该边界条件的程序代码的相关部分。

清单58显示了从输入字典中阅读列表数据。数据和文件名从dynamicMeshDict的子字典读入类型为Tuple2的列表。Tuple2类是一种数据结构,它存储两个可以是不同类型的对象。在上面的示例中,标量与translationRotationVectors的实例存储在Tuple2中,第一个实例是时间,第二个实例是嵌套向量,一个用于平移,一个用于旋转。因此,运动文件的内容直接存储在该数据结构中。为了提供对该数据的轻松访问,创建了单独的列表,如上面代码片段的后半部分所示。

fileName newTimeDataFileName
(
    fileName(SBMFCoeffs_.lookup("timeDataFileName")).expand()
);
IFstream dataStream(timeDataFileName_);
List<Tuple2<scalar, translationRotationVectors>  >  timeValues
(
    dataStream
);
times_.setSize(timeValues.size());
values_.setSize(timeValues.size());
forAll(timeValues, i)
{
    times_[i] = timeValues[i].first();
    values_[i] = timeValues[i].second();
}

从输入数据计算转换是由公共成员函数transformation()执行的,相关内容如清单59所示。第二分配内插当前时间t的位置和取向数据。所有角度都必须转换成弧度,最后用quaternions和septernions.组合了变换的表示。

scalar t = time_.value();
// -- Some lines were spared --
translationRotationVectors TRV = interpolateSplineXY
(
    t,
    times_,
    values_
);
// Convert the rotational motion from deg to rad
TRV[1] *= pi/180.0;
quaternion R(TRV[1].x(), TRV[1].y(), TRV[1].z());
septernion TR(septernion(CofG_ + TRV[0])*R*septernion(-CofG_));

构造对象当然是通过构造函数来完成的,构造函数有两个参数,如下面的清单所示。第一个是对包含tabulated6DoFMotion所需数据的字典的引用,它是数据文件的路径。将引用传递给“时间”可简化数据点之间的插值:

tabulated6DoFMotion
(
    const dictionary& SBMFCoeffs,
    const Time& runTime
);

找到负责点运动的代码后,下一步是找到从pointPatchField派生的边界条件,该边界条件执行与我们计划实现的任务类似的任务:由此我们计划建立我们自己边界条件。

10.3.2.2 调整现有边界条件

现有的oscillatingDisplacementPointPatchVectorField边界条件是导出网格运动边界条件的良好起点。它根据时间相关正弦曲线将位移值应用于存储在边界上的点中的每个值。可以在fvMotionSolver的子目录中找到oscillatingDisplacementPointPatchVectorField的源代码:

?>  $FOAM_SRC/fvMotionSolver/pointPatchFields/derived/oscillatingDisplacement/

如本章开头关于fvPatchField类型边界条件的讨论,边界条件的实际功能在evaluate()或updateCoeffs()成员函数中实现。在oscillatingDisplacementPointPatchVectorField边界条件的情况下,成员函数updateCoeffs()计算边界网格的每个点的位移向量。位移向量定义如下:

amplitude_*sin(omega_*t.value())

其中amplitude_和omega_都是从字典中读取的标量值(omega是角旋转)。这个方法的实现可以如清单60中所示。

//程序清单60
void oscillatingDisplacementPointPatchVectorField::updateCoeffs()
{
    if (this-> updated())
    {
        const polyMesh& mesh =
        this-> dimensionedInternalField().mesh()();
        const Time& t = mesh.time();
        Field<vector> ::operator=(amplitude_*sin(omega_*t.value()));
        fixedValuePointPatchField<vector> ::updateCoeffs();
    }
    fixedValuePointPatchField<vector> ::updateCoeffs();
}

updateCoeffs()方法中Field<vector> ::operator=(amplitude_*sin(omega_*t.value()));的使用Field的指定运算符将适当的置换值指定给网格的边界点。此指令之后是对父类的updateCoeffs()的调用。通过将面片的置换指定给面片点,动态网格解算器可以处理实际的网格运动。

10.3.2.3 组合边界条件

这个边界条件的工作版本可以在随书分发的示例代码存储库中找到,该存储库与前一章的再循环控制边界条件一起捆绑到库primerBoundaryConditions.中。 为了更容易地隐藏,您可能希望在文本编辑器中打开即用网格运动边界条件的源代码,同时遵循此处描述的步骤。 与其他编程示例类似,本示例的边界概念应该使用wmake libso编译到动态库中。 像往常一样,第一步是创建一个新目录来存储边界条件:

?>  mkdir -p $WM_PROJECT_USER_DIR/applications/tabulatedRigidBodyDisplacement
?>  cd $WM_PROJECT_USER_DIR/applications/tabulatedRigidBodyDisplacement

下一步是将oscillatingDisplacementPointPatchVectorField复制到新目录:

?>  cp $FOAM_SRC/fvMotionSolver/pointPatchFields/derived/oscillatingDisplacement/* .

若要保持tabulatedRigidBodyDisplacement边界条件的名称正确,所有出现的oscillatingDisplacement都必须由tabulatedRigidBodyDisplacement取代。这既适用于文件名,也适用于源文件本身内部的匹配。删除 *.dep文件后,必须相应地重命名剩余的C和H文件。

?>  rm *.dep
?>  mv oscillatingDisplacementPointPatchVectorField.H tabulatedRigidBodyDisplacement.H
?>  mv oscillatingDisplacementPointPatchVectorField.C tabulatedRigidBodyDisplacement.C
?>  sed -i "s/oscillating/tabulatedRigidBody/g" *.[HC]

首先要检查的是重命名的oscillatingDisplacementPointPatchVectorField是否仍能正确编译。要检查这一点,需要创建OpenFOAM典型的Make/文件和Make/选项文件:

?>  mkdir Make
?>  touch Make/files
?>  touch Make/options

文件的内容很短,因为边界条件仅包含一个源文件:

tabulatedRigidBodyDisplacementPointPatchVectorField.C
LIB = $(FOAM_USER_LIBBIN)/libtabulatedRigidBodyDisplacement

为简单起见,代码库中提供的所有示例边界条件都被编译到一个库中。此库的名称与上述代码段中定义的名称不同。

不过,这个源文件有很多依赖项。其他各种库及其头文件必须链接到此库:

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/fileFormats/lnInclude

LIB_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lfileFormats

通过在包含Make文件夹的目录中发出wmake libso命令,测试边界条件是否仍能编译。如果编译时没有出现任何错误和警告,请在您选择的编辑器中打开头文件和源文件,并应用以下更改。首先,头文件必须包含在tabulatedRigidBodyDisplacementPointPatchVectorField.H中:

#include "fixedValuePointPatchField.H"
#include "solidBodyMotionFunction.H"

源文件必须包含更多的头文件:

#include "tabulatedRigidBodyDisplacementPointPatchVectorField.H"
#include "pointPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "Time.H"
#include "fvMesh.H"
#include "IFstream.H"
#include "transformField.H"

所需功能的实际实现是在updatecoeffs()成员函数中完成的。 尽管这里使用了一个私有属性,但它没有实现:一个constant字典,它包含在0/目录中边界条件字典中定义的所有数据。 这个字典被添加到头文件中的私有属性中:

//- Store the contents of the boundary condition's dicitonary
const dictionary dict_;

边界条件的每个构造函数都必须初始化新的私有属性,这通常是通过调用dicitonary的null-constructor来完成的。 一个例子是从PointPatch和DimensionedField构造边界条件的构造函数,如清单61所示。 如果边界条件是通过从0目录读取特定文件来构造的,则调用以下构造函数。 事实上,边界条件的字典被传递给构造函数,并且必须存储在边界条件中以便以后处理:

// 清单61
tabulatedRigidBodyDisplacementPointPatchVectorField::
tabulatedRigidBodyDisplacementPointPatchVectorField
(
    const pointPatch& p,
    const DimensionedField<vector, pointMesh> & iF
)
:
    fixedValuePointPatchField<vector> (p, iF),
    dict_()
{}


tabulatedRigidBodyDisplacementPointPatchVectorField::
tabulatedRigidBodyDisplacementPointPatchVectorField
(
    const pointPatch& p,
    const DimensionedField<vector, pointMesh> & iF,
    const dictionary& dict
)
:
    fixedValuePointPatchField<vector> (p, iF, dict),
    dict_(dict)
{
    updateCoeffs();
}

此构造函数是唯一在构造期间调用updateCoeffs()的构造函数。updateCoeffs成员函数如清单62所示。最重要的行是定义SBMFPtr的行。它基于字典和Time对象为solidBodyMotionFunction构造autoPtr。由于此代码源自dynamicFvMesh库,因此传递给构造函数的字典是dynamicMeshDict的子字典。此子字典包含用户在dynamicMeshDict中选择的solidBodyMotionFunction所需的所有参数。由于此边界条件的运动参数定义应基于每个边界而非全局执行,因此传递给构造函数的字典应从0目录中的边界条件读取,而不是从dynamicMeshDict读取。这就是私有成员dict_的作用:它只读取一次,并在每次调用updateCoeffs()时传递给solidBodyMotionFunction。

// 清单62
void tabulatedRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
{
    if (this-> updated())
    {
        return;
    }
    const polyMesh& mesh = this-> dimensionedInternalField().mesh()();
    const Time& t = mesh.time();
    const pointPatch& ptPatch = this-> patch();
    autoPtr<solidBodyMotionFunction>  SBMFPtr
    (
        solidBodyMotionFunction::New(dict_, t)
    );
    pointField vectorIO(mesh.points().size(),vector::zero);
    vectorIO = transform
    (
        SBMFPtr().transformation(),
        ptPatch.localPoints()
    );
    Field<vector> ::operator=
    (
        vectorIO-ptPatch.localPoints()
    );
    fixedValuePointPatchField<vector> ::updateCoeffs();
}

下一行与tabulated6DofMotion中的行类似。 应用转换后的补丁点的绝对位置存储在vectorIO.中。 由于实际运动相对于先前的点位置,必须计算差值,这是直接在调用赋值算子中完成的。转换是使用间septernions进行的,在许多方面优于转换矩阵。

既然源代码已经准备好了,请再次编译库。

?>  wclean
?>  wmake libso

10.3.2.4 在示例案例中执行模拟

首先,应使用示例代码库中准备好并经过测试的代码执行示例案例中的模拟。一旦运行了分布式代码,字典中的必要输入参数以及应用于字段的修改应该是清楚的,并且可以测试所实现的边界条件。范例案例tabulatedMotionObject位于范例案例储存库中。该三维情况是新实现的平板刚体运动边界条件的功能的简单演示。图10.9显示了案例设置的示意图,包括一个立方体域和一个从中间切除的体积。由该切割过程创建的边界由movingObject面片表示。

图10.9:TabulatedMotionObject示例案例的域和移动边界说明

movingObject将根据6DoF.dat文件中常量/中包含的数据由新的边界条件进行位移。与基本的OpenFOAM情况相比,执行新的边界条件需要新的配置文件:0/pointDisplacement和常量/dynamicMeshDict文件。使用新的边界条件为pointDisplacement字段的movingObject面片设置初始边界值,如清单63所示。为了使用边界条件,需要将在运行时动态加载新库的指令添加到system/controlDict配置文件中。请注意,如果使用本节中提供的代码编译边界条件,则库的名称会有所不同(“libofPrimerBoundaryConditions.so“)。

//清单63
movingObject
{
    type    tabulatedRigidBodyDisplacement;
    value    uniform (0 0 0);
    solidBodyMotionFunction tabulated6DoFMotion;
    tabulated6DoFMotionCoeffs
    {
        CofG        ( 0 0 0 );
        timeDataFileName  "constant/6DoF.dat";
    }
}

需要dynamicMeshDict来控制网格运动解算器。在这种情况下,使用了displacementLaplacian smoother,它基于拉普拉斯方程并平滑点位移。如果周围的点不随膜片移动,近端网格将迅速变形和破碎。清单64显示了dynamicMeshDict的内容。

    // 清单64
    dynamicFvMesh dynamicMotionSolverFvMesh;
    motionSolverLibs ("libfvMotionSolvers.so");
    solver
    displacementLaplacian;
    displacementLaplacianCoeffs
    {
        diffusivity
        inverseDistance (floatingObject);
    }
}

可以通过在Simulation case目录中执行allrun脚本来运行模拟。 此脚本执行所有必要的步骤,如网格生成,以轻松完成此示例。 求解器MovedynamicMesh工具只调用动态网格例程,没有任何与流相关的计算,这使得它相对于具有动态网格功能的通常流求解器来说相对较快。 除了执行网格运动操作,它还执行大量与网格相关的质量检查。

案例的后处理非常简单,可通过目视进行:可以在Paraview中检查该案例。为案例设置动画时,会显示内部面片在域内翻滚和摇摆,类似于输入文件中提供的数据。使用带有“Crickle Clip2”选项的“Clip”过滤器将域切成两半,可以阐明与动态网格有关的一些有趣的事情:由边界条件施加的点位移由网格运动求解器耗散的方式。拉普拉斯方程将位移耗散到内部网格点,相应地移动这些点,并在面片平移和旋转时保持良好的网格质量。图10.10显示了变换前后最终面片位置的图像。

图10.10 边界转换前后

总结

本章介绍了OpenFOAM中FV法和网格运动边界条件的设计与实现。 这两个边界条件族分别被建模为两个分离的类层次结构,FVPatchField和PointPatchField。 动态多态性允许现有实现的不同扩展和组合。 将边界条件与存储在内部网格点中的物理场一起实现为类层次结构,允许OpenFOAM库在运行时确定边界条件的类型,而无需在每次为场设置新的边界条件时重新编译代码。 动态库的加载机制允许将边界条件的新实现编译到单独的库中(例如FinitEvolume库)。 这反过来又产生了一种直接添加新边界条件的方法,而不必重新编译数值库。 这些优点允许程序员通过开发(和共享)自我维持的库来扩展边界条件和OpenFOAM的其他部分。

开发新的边界条件涉及到在类层次结构中找到一个入口点,从这个入口可以扩展层次结构。 选择类作为边界条件进一步发展的起点取决于所需边界条件实现的功能,以及是否已经存在类似的边界条件。 计算部分边界条件位于Updatecoeffs成员函数中。 除了适当地重命名源文件和修改新类以进行继承之外,程序员很可能会完成大部分工作来实现这个函数的主体。