OpenFOAM无反射边界条件源码学习

2022-06-24 00:48:57

求解流声分解法的shen方程时，需要对变量的边界作无反射处理。OpenFOAM提供的无反射边界条件有advective和waveTransmissive这两种，但这两种似乎都不能满足笔者的需要，可能要进行某些修改。在这之前，先来学习一下这两种无反射边界条件的代码以及它俩的区别，这对自己的运用起到帮助。这里会参考日本博主的博文。他用的是英文，而且只起了个开头，后面重要的没写（等他好久了都没补上）。我写点中国人看得懂的，顺便把后面的内容也补上。

1. 无反射边界条件

基于有限体积法，OpenFOAM的边界条件就是对边界上的面元通量进行计算（拟合）。在有限体积法里，面上的值通常用相邻两单元的体元值（cell value）插值近似求得。但是在边界上，边界面元只有一边和边界单元体接壤，另一边是计算域之外，所以需要根据不同的需求来对这块边界面元上的值进行设置或者计算。常用边界条件的有fixedValue（固定边界元上的值），zeroGradient（将接壤单元的体元值赋值到面元上）等。在高精度的有限差分方法中，边界可以通过边界外几个点进行插值拟合。而在有限体积法里，受到方法本身的限制，计算复杂的高精度的边界条件就没那么容易了。一些有特别需要的边界条件则要求解对应的偏微分方程才能得到。本文所讲到的无反射边界条件，就属于这一类。当我们需要波（声波、涡波、熵波）流出边界而不会反射回来，对计算域内造成影响时，就需要应用无反射边界条件。

2 advective

2.1 控制方程解析

advective在流体力学里是对流的意思。在边界条件上，主要体现在流体的流出。通过解方程：

来得到面上的通量。其中是wave velocity (see advectiveFvPatchField.H)，由函数advectiveSpeed()函数计算得到。是面元法向量。就是需要无反射处理的流场变量。形式上相当于物质导数，又称全导数、随体导数等，本质上是拉格朗日法的概念之一。它是运动的流体微团的物理量随时间的变化率，它等于该物理量由当地时间变化所引起的变化率与由流体对流引起的变化率的和。然而它和一般我们说的物质导数又有一些区别，因为它只关注边界面上的情况，即仅关注流体在边界这一面上的对流。在计算上，它其实是fixedValue和zeroGradient的结合，当：

时，为fixedValue。而当：

时，相当于zeroGradient。这两项其中一个为零，另一个也为零。可以理解为流体微团的物理量随时间的变化率，等于该物理量由当地时间变化所引起的变化率与由流体通过边界流出引起的变化率的和。拉格朗日法关注单个流体微团。通俗来讲，就是流体微团通过边界面时，其本身对应的物理量不发生改变（全导数为0）。信息“完整地”通过了边界，也就没有造成反射，从而影响到计算域内的情况。想要更直观地了解，可以看看博文里面提供的视频。

的计算由以下代码进行：

template<class Type>
Foam::tmp<Foam::scalarField>
Foam::advectiveFvPatchField<Type>::advectionSpeed() const
{
    const surfaceScalarField& phi =
        this->db().objectRegistry::template lookupObject<surfaceScalarField>
        (phiName_);

    fvsPatchField<scalar> phip =
        this->patch().template lookupPatchField<surfaceScalarField, scalar>
        (
            phiName_
        );

    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
    {
        const fvPatchScalarField& rhop =
            this->patch().template lookupPatchField<volScalarField, scalar>
            (
                rhoName_
            );

        return phip/(rhop*this->patch().magSf());
    }
    else
    {
        return phip/this->patch().magSf();
    }
}

phi为速度通量U()&magSf()，在这里只是程序判定phi是否为密度流量的。如果是，那就要在输出的时候除以面元上的密度值（rhop）。因而这个Un的计算只是把边界面的速度通量除以面矢量，还原成面上的法向速度罢了。

2.2 和

lInf和fieldInf的设置不是必须的。头文件原话：Additionally an optional mechanism to relax the value at the boundary to a specified far-field value is provided which is switched on by specifying the relaxation length-scale \c lInf and the far-field value \c fieldInf. 不设置时，lInf默认负无穷，fieldInf默认为零。lInf是松弛用的，松弛系数定义如下：

所以不设置lInf和fieldInf时，松弛系数默认是零。

2.3 计算过程

要理解边界条件是怎么被应用于边界单元面上值，首先需要知道五个主要的函数：updateCoeffs(), valueInternalCoeffs(), valueBoundaryCoeffs(), gradientInternalCoeffs()和gradientBoundaryCoeffs()。这里先留个坑，下次再详细讲解。在advective这个边界条件中，程序会根据时间离散格式的不同来计算修正的系数（valueFraction），这个valueFraction会流到后面四个函数中作为计算的要素之一。valueFraction的值为：

当valueFraction为1时，对应fixedValue边界条件；为0时对应fixedGradient（不一定zeroGradient）。 α是wave coefficient，计算公式为：

deltaCoeffs()是face-cell distance coeffcient，是和距离有关的系数，与网格和梯度离散格式有关。

时间格式的不同会影响valueFraction的计算方式，体现出方程中时间导数的影响；而deltaCoeffs()则体现了物理量对面元法向矢量的求导。由此完成对边界的修正。

3. waveTransmissive

waveTransmissive继承了advective的类，控制方程为：

此时程序里的advectionSpeed()计算的是(U_n+c)，也就是在advective的基础上，再加上当地声速。声速c的计算公式如下：

所以在使用时，需要给定psi和gamma：

    <patchName>
    {
        type            waveTransmissive;
        phi             phi;
        psi             psi;
        gamma           1.4;
    }

而advective只要type和phi即可。

参考资料：

http://caefn.com/openfoam/bc-advective-wavetransmissive

https://cpp.openfoam.org/v9/waveTransmissiveFvPatchField_8H.html

码农公寓