From 384bed27f78dd59d3118432d17c21ed5746dec97 Mon Sep 17 00:00:00 2001 From: skf0558 Date: Sat, 15 Aug 2020 12:30:13 +0800 Subject: [PATCH 1/3] =?UTF-8?q?=E8=AE=A4=E9=A2=8638=20The=20MRF=20method.m?= =?UTF-8?q?d.=20##=2038=20=E5=A4=9A=E5=8F=82=E8=80=83=E7=B3=BB=E6=A8=A1?= =?UTF-8?q?=E5=9E=8B=20(MRF)=E6=96=B9=E6=B3=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ​ **MRF**方法允许在没有实际旋转网格的情况下模拟旋转机械。通过使用多参考坐标系，可以用静态网格来模拟问题。与动网格模拟相比，尽管这种简化会带来一定的建模误差，但它降低了复杂性和计算更快，同时在全尺度上保证足够的精度，在合适的条件下该方法是合理的。MRF方法的理论背景参见第65节。 ### 38.1用法 ​ **MRF**方法应用于单元区。这在**MRFProperties**的字典中通过**cellZone**关键字清楚地表明了这一点，参见列表 243。除了在哪使用**MRF**方法，我们还可以使用**MRF**的***active***关键字启用/禁用它。 ​ 另一个重要的输入是非旋转面的列表(恰当的名字是**nonRotatingPatches**)，当固体的旋转时， 根据参照系的旋转应用到所有面patches。然而，可能有这种算例，patches在MRF单元区内，实际上它们是静止的，因此可以选择从**MRF**方法的使用中排出单独的patches。 ​ 最后，在列表243中，指定了旋转参考系性质。这是通过提供旋转轴(**axis**)和空间中的一个点(**origin**)。origin点必须位于旋转轴上。关键字**omega**用于指定旋转速度。 ​ 图97显示了一个搅拌槽网格的单元区域。**MRF**方法所在的单元格区域显示为白色线框。请注意，圆柱形的这个区域与旋转轴对齐。 ​ 如果我们将圆柱体从最底部延伸到最顶部，那么，定子底部到顶部patches需要输入到非旋转patches列表中。 --- ...\257\221\347\253\240\350\212\202\350\256\244\351\242\206.md" | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git "a/\347\277\273\350\257\221\347\253\240\350\212\202\350\256\244\351\242\206.md" "b/\347\277\273\350\257\221\347\253\240\350\212\202\350\256\244\351\242\206.md" index b48643d..564ae82 100644 --- "a/\347\277\273\350\257\221\347\253\240\350\212\202\350\256\244\351\242\206.md" +++ "b/\347\277\273\350\257\221\347\253\240\350\212\202\350\256\244\351\242\206.md" @@ -173,7 +173,7 @@ - 37.3 Arbitrary Coupled Mesh Interface - ACMI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 - 37.4 Avoiding errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 ### 38 The MRF method 192 -- 38.1 Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 +- 38.1 Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 skf0558 - 38.2 Avoiding errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 ### 39 The fvOption framework 194 - 39.1 Controlling space & time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 -- Gitee From e2bfa171defd575a7bbfbee57290ade68e1dad3e Mon Sep 17 00:00:00 2001 From: skf0558 Date: Sun, 16 Aug 2020 10:53:36 +0800 Subject: [PATCH 2/3] 38.MRF --- 38.MRF.md | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 38.MRF.md diff --git a/38.MRF.md b/38.MRF.md new file mode 100644 index 0000000..18da767 --- /dev/null +++ b/38.MRF.md @@ -0,0 +1,61 @@ +# 38 MRF多参考系模型方法 + +​ **MRF**方法允许在没有实际旋转网格的情况下模拟旋转机械。通过使用多参考坐标系，可以用静态网格来模拟问题。与动网格模拟相比，尽管这种简化会带来一定的建模误差，但它降低了复杂性和计算更快，同时在全尺度上保证足够的精度，在合适的条件下该方法是合理的。MRF方法的理论背景参见第65节。 + +## 38.1用法 + +​ **MRF**方法应用于单元区。这在**MRFProperties**的字典中通过**cellZone**关键字清楚地表明了这一点，参见列表 243。除了在哪使用**MRF**方法，我们还可以使用**MRF**的***active***关键字启用/禁用它。 + +​ 另一个重要的输入是非旋转面的列表(恰当的名字是**nonRotatingPatches**)，当固体的旋转时， 根据参照系的旋转应用到所有面patches。然而，可能有这种算例，patches在MRF单元区内，实际上它们是静止的，因此可以选择从**MRF**方法的使用中排出单独的patches。 + +​ 最后，在列表243中，指定了旋转参考系性质。这是通过提供旋转轴(**axis**)和空间中的一个点(**origin**)。origin点必须位于旋转轴上。关键字**omega**用于指定旋转速度。 + +``` +zone1 +{ + cellZone rotor; + active yes; + // Fixed patches (by default they 'move' with the MRF zone) + nonRotatingPatches (); + origin (0 0 0); + axis (0 0 1); + omega table + 2( + (0 0.0) + (0.75 20.0) + ); +} +``` + +​ 表243：在**MRFProperties**字典中为MRF方法指定必要的输入 + +​ 图97显示了一个搅拌槽网格的单元区域。**MRF**方法所在的单元格区域显示为白色线框。请注意，圆柱形的这个区域与旋转轴对齐。 + +​ 如果我们将圆柱体从最底部延伸到最顶部，那么，定子底部到顶部patches需要输入到非旋转patches列表中。 + +![images](89.PNG) + +​ 图97：带有Rushton叶轮的导流式搅拌槽。定子patches显示为灰色，转子显示为红色。白色线框指转子区域的边界。对于转子区的所有单元采用**MRF**方法 + + + +## 38.2 避免错误 + +### 38.2.1 非旋转patches + +**周期性任意匹配网格界面(AMI) 和 MRF** + +可以发现，AMI型的 patches要添加到非旋转patches列表中,参见37.4.1 + +**周期性(Cyclic) patches 和MRF** + +可以观察到，当使用cyclic型patches时，例如当模拟半个搅拌槽时，如图98，需要将cyclic型patches添加到**MRFProperties**文件中的非旋转patches列表中。 + +![images](90.PNG) + +图98：半个Rushton叶轮的导流式搅拌槽，cyclic型patches显示为彩色线框，所有其他的patches显示为彩色面。 + + + + + -- Gitee From 753d3f51c44de34463dcf098c735d03f374f2cd5 Mon Sep 17 00:00:00 2001 From: skf0558 Date: Tue, 18 Aug 2020 08:18:24 +0800 Subject: [PATCH 3/3] 39.fvOption framework --- 39. fvOption framwork .md | 242 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) create mode 100644 39. fvOption framwork .md diff --git a/39. fvOption framwork .md b/39. fvOption framwork .md new file mode 100644 index 0000000..9af3877 --- /dev/null +++ b/39. fvOption framwork .md @@ -0,0 +1,242 @@ +# 39 fvOption架构 + +​ **fvOption**架构处理数值流动模型的源项和约束。**fvOption**框架可以将各种约束或源项插入到现有求解器中，而无需对求解器进行任何修改。此外一般的源项和约束，有许多特殊的源项代表一个特定的物理模型，例如， 多孔区对动量方程的影响。 + +​ **fvOptions**框架是在2013年随着**OpenFOAM-2.2.0**[108]的发布而引入的，并沿用至今。 + +### 动机 + +​ **fvOption**架构的使用在一个例子中得到了最好的解释。方程(57)表示一个通用标量C的对流输运方程，例如，被动示踪剂浓度。在RHS(公式右侧)，我们看到通常的线性化源项。在例子中，如果示踪剂通过入口进入模拟区，那么所有，那么一切都好，并且对示踪剂浓度C来说一个入口边界就足够了。然而，如果示踪剂是通过探针引入的，我们不想用网格[109]来解决这个问题，那么我们需要一个机制来引入示踪剂浓度C到模拟区。这就是**fvOption**显身手的地方，它使我们，在探针尖端的位置，可以通过指定场C的注入速率或一个固定值来解决问题。 + +$\frac{\partial C}{\partial t} + \nabla \cdot(C u) = S_u +S_p C \tag{57}$ + +​ **fvOption**架构提供了一个固定值约束和一个半隐式源项，可以用它来建模示踪剂生成探针。 + +### 应用 + +​ 在下面的列表 244中，在**OpenFOAM**中的一个可压缩流动求解器的能量方程。在这个列表中， + +对**fvOptions**架构的所有调用都用红色、绿色和蓝色标记。可以看到**fvOptions**架构既可以用于某个场量本身的控制方程，也可以用于控制方程求解后，来计算场量 。 + + + + { + volScalarField & he = thermo .he (); + fvScalarMatrix EEqn + ( + fvm::ddt(rho , he) + fvm::div(phi , he) + + fvc::ddt(rho , K) + fvc::div(phi , K) + + ( + he.name () == "e" + ? fvc::div + ( + fvc::absolute(phi / fvc::interpolate( rho), U), + p, + " div(phiv ,p)" + ) + : -dpdt + ) + - fvm :: laplacian(turbulence->alphaEff() , he) + == + fvOptions(rho, he) + ); + EEqn.relax(); + fvOptions.constrain(EEqn); + EEqn.solve(); + fvOptions.correct(he); + thermo.correct(); + } +列表244：可压缩单相求解器*rhoPimpleFoam* 能量方程的**EEqn.H**文件 + +​ 一个使用**fvOptions**架构的求解器，创建一个**fvOptionList**类型的对象，它是一个**fvOptions**列表。这也是为什么对**fvOptions**架构的所有调用都必须传递参数的原因。否则，架构将无法将选项指定给各自的方程。在一个单相模拟多孔区域的算例中，我们可以使用**fvOptions**架构来考虑动量方程中的附加流动阻力，或者对湍流模型方程施加某类源项。因此，必须遍历选项列表以找到并应用适当的**fvOption**到相应的模型方程。 + +​ 在下面的列表245中，用**fvOptionList::constrain()**方法作为**OpenFOAM**使用**fvOptions**的例子。对每个选项遍历所有选项的列表，检查每个选项是否适用所提供的模型方程 。如果适用，决定该选项是否启用，即它是否处于活动状态，如果是，则通过调用**fvOption**类中的**constrain()**方法，使用该选项，即调用源项本身的**constraint()**方法。 + +``` +templatevoid Foam::fv::optionList::constrain(fvMatrix& eqn) +{ + checkApplied(); + forAll (*this, i) + { + option& source = this->operator[](i); + label fieldi = source.applyToField(eqn.psi().name()); + if (fieldi != -1) + { + source.setApplied(fieldi ); + if (source.isActive ()) + { + if (debug) + { + Info<<" Applying constraint" <divDevReff(U) + == + fvOptions(U) +​ Listing 246: pimpleFoam 求解器UEqn.H文件中的动量方程 + +#### 定系数 + +​ **fixedCoeff**模型计算与速度和平方速度。该模型具有两个模型常数: $\alpha$ 和 $\beta$ + +$S = - \rho_{Ref} (\alpha + \beta |u|)u \tag{58} $ + +或者重排一下： + +​ $S = - \rho_{Ref} (\alpha |u| + \beta |u|^2)u $ + +​ 在**OpenFOAM**中α和β是矢量量。另外还有一个参考坐标系统为多孔区，可考虑各向异性孔隙度。各向同性孔隙度是一种特殊情况，通过选择α和β的各分量相等来表示。 + +#### 幂律 + +​ **powerLaw**模型计算动量贡献S，它与模型常数$C_0$和速度的$C_1$次方成正比。**powerLaw**模型不支持各向异性。 + +#### DarcyForchheimer + +​ **DarcyForchheimer**模型与**fixedCoeff**模型非常相似。它也有两个与线性和平方速度成比的项。这个模型是达西模型$U = -\frac{\kappa}{\mu}\frac{dp}{dx}$的组合，对层流是有效，并且可以扩展更高的雷诺数，被称为**Forchheimer**模型$-\frac{dp}{dx} = \frac{\kappa}{\mu}U + \frac{\rho}{\kappa_2} U^2$。达西模型与κ成正比，即多孔介质的渗透率，[κ]=m2。**Forchheimer**模型引入了$k_2$，惯性渗透率，[$k_2$]=m。 + +​ **OpenFOAM**实现**DarcyForchheimer**模型有两个模型常数。这个达西系数d是渗透率d倒数$d=1/κ$，而**Forchheimer**系数f是惯性渗透率的倒数$f=1/k_2$ + +$S = -(\mu d + \frac{1}{2}\rho |u|f)U \tag{60}$ + +或者重排为： + +​ $S = -\rho (\nu d |u| + \frac{1}{2}|u|^2 f)\rm e_u$ + +​ 在OpenFOAM中，d和f是向量量，也是**fixedCoefff**模型的系数。除了多孔区域的参考坐标系外，还可以考虑各向异性孔隙度。各向同性孔隙度通过选定d和f的各分量相等，作为各向异性的一个特例。 + +### 39.3.2相位稳定 + +​ 在**OpenFOAM-6**[110]中引入了一个**fvOption**源项，当相分数低于某个阈值时，可以使输运方程数值计算稳定。当相分数低于某一阈值，它将隐式源项添加到输运方程所有单元。 + +​ 用相分数场表示的输运方程的求解在数值上变得更加复杂，当某些单元的相分数场趋于零时，就更难。添加隐式源项在这些单元中，可以确保离散化的方程组不会变得不适定。当一个输运方程由于相位分数趋于零而变得不适定，线性方程组求解就会因浮点错误失败。 + +​ **列表247**显示了这个源项类相关方法的代码。max()语句确保该源项仅英用于相分数**alpha**小于阈值 **residualAlpha**的单元。 + +``` +template void Foam::fv::PhaseLimitStabilization::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi +) +{ + const GeometricField& psi = eqn.psi(); + uniformDimensionedScalarField & rate = + mesh_.lookupObjectRef(rateName_); + +eqn -= fvm::Sp(max( residualAlpha_ - alpha, scalar(0))*rho*rate,psi); + +} +``` + +​ 列表247：类**PhaseLimitStabilization**的**addSup()**方法 + +​ 从**列表247**中的代码可以看出，用**uniformDimensionedScalarField rate**来计算 + +隐式源项。**uniformDimensionedScalarField rate**可以由**codedFunctionObject**创建然后在网格上注册。仅当**UniformMDimensionedScalarfield rate **在网格注册时，**PhaseLimitStabilization**源项能通过网格的**lookup()**方法访问它。**列表248**展示了如何做到这一点的例子。去注册计算速率场的相关行代码红色高亮显示。如果省略这一行，那么**fvOption**源项将因尝试查找不存在的场而中止。有关对象注册表的详细信息，请参见第**57.7**节 + +``` +functions +{ + rate + { + libs (" libutilityFunctionObjects.so"); + type coded ; + name rate ; + codeRead + #{ + static autoPtr < uniformDimensionedScalarField > pField ; + if (! pField.valid()) + { + pField.set + ( + new uniformDimensionedScalarField + ( + IOobject + ( + "rate.water", + mesh().time().timeName() , + mesh() , + IOobject::NO_READ , + IOobject::NO_WRITE + ), + dimensionedScalar("rateInit", dimensionSet(0,0,-1,0,0,0,0),1.0) + ) + ); + } + uniformDimensionedScalarField& rate = pField(); + rate.checkIn(); + #}; + } +} +``` + +​ 列表248：使用文件**controlDict**中的**codedFunctionObject**计算速率场。 + + + -- Gitee