openfoam中cyclic周期性边界的问题

余正东

在使用interfoam或multiphaseinterfoam求解器中，计算油滴气泡上升的问题，上下左右前后均使用周期性边界，根据观察结果cyclic函数类型（type），能在边界上传递相分数项、压强U和P，但不能继续向空间中以上区域传播，而且一运算到!边界处时间步长就会变得特别小（可能也是发散），请问这种情况怎么解决啊~~!
0_1510990756944_1a74d12a-08cd-4dff-94b2-3cddbb4d8bb1-图片.png
0_1510990197002_4997256f-4599-490c-9876-d469eb0a7590-图片.png
0_1510990874481_f1bc7158-1b28-4e37-ab6c-955f0ff693fc-图片.png
0_1510990518747_5168e382-3fb0-4655-a1c3-9611656d2fc6-图片.png
0_1510990792142_c024618e-6fe0-424b-bc10-897b6c22bd23-图片.png

xiaofenger

我也曾遇到非常类似的问题，我就有一个液滴在一段小管子（两边cyclic）中移动，速度不够大的时候，就发现液滴移动到一头出不去，需要等很久才能出去，明显是有问题；然而速度达到一定程度，看起来就能自由出去然后从cyclic的另一端进来，不过我也不知道为什么。。。

random_ran

不知道把计算域扩大效果会怎么样？

是要观察油滴气泡的上升速度和其它参量的关系？一定要用cyclic 边界吗？或者说为什么要用周期边界条件？

余正东

@xiaofenger 这个已经困惑我很长时间了，可是我的油滴还有气泡都是自然上升启动的，也就是说到达边界时速度很小，如果可以的话，可以私下交流一下吗，谢谢！

余正东

@random_ran 是的，重在研究相关关系，而且我是油滴气泡群，所以得用周期性边界！你说的扩大计算域，是为了让气泡到达边界时速度更大吗

random_ran

@余正东我不是多相流的专家，恐怕很难直接帮你解决问题。

不过我有用过周期边界条件。计算域如果太小，会产生一些无法解释的现象。而且这些现象和真实世界有很大差异。通常来说计算域要大到一定程度之后，才能消除这种数学概念上的“周期边界条件”所带来的人为干扰。

余正东

@random_ran 不管怎么样，都要谢谢你的解答！我先试一试，等有结果了再与你交流~

东岳

网格多少万？3D？

余正东

@李东岳李老师你好，我是80* 160* 80=102.4万个网格，（3D格式）

东岳

网格太多了我这跑不动，如果可以减少到1万个网格我可以跑跑看看。

余正东

@李东岳好的，感谢李老师的回复，可以简化成二维问题试一试~

xiaofenger

@余正东交流当然可以啊，不过关于cyclic边界条件我也不是特别懂。。。

余正东

@xiaofenger 你好，我也是这方面的新手，我就是想问问你的算例最后怎么跳出时间步长特别小的情况，我今天又重新跑了算例，发现还是有问题！qq：768620698，谢谢！

余正东

@xiaofenger 你好，我还想问一下请问一下你用的是cyclic还是cyclicAMI边界，另外你用的是VOF模型吗？谢谢~

余正东

@YHDTHU 请问一下这个可能是由于MULES算法的原因吗，导致在周期性边界处，不能收敛

yhdthu

@余正东不知道你在边界处有没有加密网格？

余正东

@yhdthu 在边界处局部加密吗

yhdthu

@余正东可以把网格介绍一下

余正东

@yhdthu 我的算例很简单

blocks
(
    hex (0 1 2 3 4 5 6 7) (160 320 160) simpleGrading (1 1 1)
);

程迪

cyclic BC是个大坑，因为OF的理论和编程经常不自洽。感觉网上从来没说清楚过。刚才我才搞明白。

首先，face addressing的owner, neighbour和ldu addressing的lower, upper并不天然等价。face addressing描述的的网格的几何拓扑关系。而ldu addressing需要描述的的是数值拓扑关系，如果没有耦合边界，二者是等价的，但是如果有耦合边界，就有了face addressing描述不了的off-diagonal term。这就有两种处理方式：
在OF+1706中，有两种lduAddressing，包括fvMeshLduAddressing和fvMeshPrimitiveLduAddressing。
- 实际上fvMeshLduAddressing是将face addressing和ldu addressing等价了，内部面的owner就是lower，neighbour就是upper。这也是最最常用的类型，可以称之为普通ldu。
- 而fvMeshPrimitiveLduAddressing（这里的Primitive的意思是not developed or derived from anything else ）中的lower可以添加owner之外的相互作用关系。但是这个类位于overset目录下，可见它主要用于重叠网格的情况，可以称之为overset ldu。
而coupled BC是实现在普通ldu中。依靠两个类：coupledFvPatch和coupledFvPatchField，其中关键是他们各自继承了lduInterface和lduInterfaceField，并实现了相关函数。所以其实interface和coupled BC是几乎同样的概念。
从概念上讲，对于普通LDU，在有interface的情况下虽然在lduMatrix中存了upper,diag,lower矩阵系数，但是还有额外耦合项，也就是说$A=L+D+U+C\ne L+D+U$，$C$就是额外耦合项，由于数据结构设计上的限制，不能存到普通ldu的矩阵中，需要单独处理。
处理方式也比较简单，$C$因为也是代表cell-cell相互作用，所以主要是非对角的，因为对角的部分完全可以包含在$D$中，而且源项部分也可以容纳在原来的源项里，可以分裂为上下两部分，$C=C_L+C_U$。在lduMatrix.Amul()的操作中，计算$A\cdot x$的操作被分为$D\cdot x+(L+U)\cdot x +C_U\cdot x_n $：
- initMatrixInterfaces: $C_L\cdot x_o$，其中$x_o$是$x$中属于耦合边界的内侧单元，计算结果要送到另一侧去；
- all cells: ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];: $D\cdot x$
- all faces: ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];: $L\cdot x$
- all faces: ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];: $U\cdot x$
- updateMatrixInterfaces: $C_U\cdot x_n$，其中$x_n$只包含$x$中属于耦合边界的外侧单元；
- 这种$C\cdot x$叫做patch Contribution。 $C_L, C_U$叫interfaceUpper, interfaceLower。
问：为什么$C_L,C_U$要分开处理？答：因为在存在processor边界时，$x_n$并不在本地，要通过通信获取结果，或者获取系数。

问：$C_L$和$C_U$存在哪儿？答：

class fvMatrix
:
    public refCount,
    public lduMatrix
{
//...
        //- Boundary scalar field containing pseudo-matrix coeffs
        //  for internal cells
        FieldField<Field, Type> internalCoeffs_; //注意这里不是引用，是真货

        //- Boundary scalar field containing pseudo-matrix coeffs
        //  for boundary cells
        FieldField<Field, Type> boundaryCoeffs_; //都初始化为0；
//...
}

coupled BC, interface, constraint BC等的关系：
interface就是coupled BC，constraint BC是OF文档中的分类，其实包含了empty, wedge,symmetry等几何约束类的BC和coupled BC。coupled BC中比较重要的两类是cyclic和processor。

总的来看cyclic边界条件本身的实现应该没有问题，还挺巧妙的。

回到这个问题，我觉得可能是CFL数计算的部分可能在cyclic BC处有bug。
Courant数计算方式来看，普通Courant数和其他一样，额外多了一个Interface Courant数。里面多了个这玩意儿：

    scalarField sumPhi
    (
        mixture.nearInterface()().primitiveField()
       *fvc::surfaceSum(mag(phi))().primitiveField()
    );

    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

    meanAlphaCoNum =
        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();

// mixture.nearInterface()()，定义是：

 Foam::tmp<Foam::volScalarField>
Foam::interfaceProperties::nearInterface() const
{
     return pos(alpha1_ - 0.01)*pos(0.99 - alpha1_);
}

可能是相边界速度太大了吧。估算一下有上万了，你的平均Courant这么小，最大Courant这么大。。。你又是这么均匀的网格。。

dzw05

@程迪大佬可以详细的解释一下cyclic的计算过程中吗？比如这个CL和CU到底对应哪一项，initMatrixInterfaces怎么计算了Cl*x0.