OpenFOAM无反射边界条件源码学习_Kino Chan的博客

1. 无反射边界条件

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

2 advective

2.1 控制方程解析

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

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

时，为fixedValue。而当：

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

$U_n$ 的计算由以下代码进行：

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());
        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即可。 
参考资料： 
Non-Reflecting Boundary Conditions in OpenFOAM | CFD WITH A MISSION 
OpenFOAM: src/finiteVolume/fields/fvPatchFields/derived/waveTransmissive/waveTransmissiveFvPatchField.H File Reference
                    求解流声分解法的shen方程时，需要对变量的边界作无反射处理。OpenFOAM提供的无反射边界条件有advective和waveTransmissive这两种，但这两种似乎都不能满足笔者的需要，可能要进行某些修改。在这之前，先来学习一下这两种无反射边界条件的代码以及它俩的区别，这对自己的运用起到帮助。这里会参考日本博主的博文。他用的是英文，而且只起了个开头，后面重要的没写（等他好久了都没补上）。我写点中国人看得懂的，顺便把后面的内容也补上。无反射边界条件基于有限体积法，OpenFOAM的边界条件就是
要求：具有C++11支持和OpenFOAM v6的g++ 。 其他编译器未经测试，但可以正常工作。
 第1步，按照官方安装OpenFOAM v6。 之后进行测试以确保安装正常（例如，通过运行一些串行和并行教程案例）。
 第2步，无论是在便携式计算机还是HPC群集上，只需键入三个命令。 
git clone https://github.com/ChenguangZhang/sdfibm.git
cd sdfibm/src
 求解器二进制文件是./src/sdfibm 。 建议您将其软链接到系统范围的路径（例如sudo ln -s ~/sdfibm/src/sdfibm /usr/local/bin/sdfibm ），然后只需键入sdfibm即可运行仿真而无需指定整个二进制路径。
 坐标系的约定是，我们
				因研究需要，特写一篇非单一变量（ρU, rhoU）边界条件的实现过程。在解可压NS方程时，rhoCentralFoam(解析)对动量方程的ρU进行直接的插值求解。换句话说，以ρU作为一个守恒变量，在方程中先进行求解，再分别更新ρ和U来。
我们知道， 有限体积法的边界条件就是计算边界面上的通量。在OpenFOAM中单一变量ρ和U的边界条件都必须要分别设置好，这意味着它们边界面上的通量如何计算将由我们来进行决定，其中最常用的有fixedValue（第一类边界条件）和zeroGradient（第二类边界条件）。*
				转载链接：http://xiaopingqiu.github.io/2016/04/02/Boundary-conditions-in-OpenFOAM2/
本篇在上一篇的基础上来解读 OpenFOAM 中的基础边界条件。基础边界条件一般包括三类，一是Dirichlet 边界，二是 Neumann 边界，三是混合 Dirichlet 和 Neumann 的边界。 
fi...
zeroGradient：边界物理量梯度为0。（Neumann边界）
fixedValue：边界物理量为定值。（Direchlet边界）
empty：在利用2D模型模拟3D情况时，侧边同城设定为empty边界，在计算过程中，此边界并不参与求解。
注：当fixedvalue的值与internalField的值相同时，fix...
				注：如有翻译不妥，还请见谅 
翻译自：http://openfoamwiki.net/index.php/HowTo_Adding_a_new_boundary_condition
首先请看：http://openfoamwiki.net/index.php/Contrib_groovyBC
如果你没有合适的边界条件，请进行下面的步骤： 
选择一个和你需要的边界条件相近的已有边...
				之前参考了两端固定的有反射边界的代码：
http://wenku.baidu.com/link?url=awvj-UI_m5QZHiN_iCrWZfPLnSBneXFq8SMP52PGUDwNfcN7Q_GfSnk0GFXGE0QCb0uY17W-4LrK96jgyuWGfXW2l6Qpsugx4rLpqSnQFcO
function reflection_1dfdtd
clear
				在 OpenFOAM 中，可以使用 `wallFlux` 函数来计算某个边界周围的通量。该函数需要在求解器中引入 `wallFunctions` 库，并在求解器的代码中调用。下面是一个示例代码片段：
#include "wallFunctions/wallFunctions.H"
// ...
fvScalarMatrix UEqn
    fvm::div(phi, U)
  + turbulence->divDevReff(U)
  + wallFlux(U)
在上面的代码中，`wallFlux(U)` 表示计算边界 `U` 附近的通量。这个函数会自动识别边界条件并计算相应的通量，无需手动指定边界条件。计算结果将返回一个 `fvMatrix` 类型的对象，可以在求解器中使用。
需要注意的是，`wallFlux` 函数只能在求解器中使用，无法在其他地方调用。如果需要在其他程序中计算通量，可以参考 `wallFlux` 函数的实现方式，在代码中手动计算通量。