非稳态计算向压力方程中添加源项

Tong

今天发现这个问题折腾了小半天也没解决，恳请大家帮忙分析一下：

非稳态计算过程中需要连续的向模型中注入流体，由于尺度的问题无法使用“显式”的入口，而是使用一个“隐式”的方法：添加一个体积源项到压力方程中实现注入；

体积源项定义：

volScalarField wellSource(
            IOobject("wellSource",
                     runTime.timeName(),
                     mesh,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE),
            mesh,
            dimensionedScalar("Scalar", dimVol / dimTime, Zero)
    );

需要添加体积源项的网格对应的wellSource有数值，否则为零，注入量是一个固定不变的值。

整体求解流程类似PISOFoam 前后代码就没有列出来，压力方程源项添加如下：

fvScalarMatrix pEqn((fvm::laplacian(rTU, p) == fvc::div(phiHbyA))-(wellSource/mesh.V()));

但是在计算过程出现问题是：每一个时间步求解一次的压力方程后压力都会增加，对应的库朗数也一直增加（如下）：

Courant Number mean: 0.000925584 max: 0.19292
Courant Number mean: 0.00103634 max: 0.331209
Courant Number mean: 0.00132796 max: 3.74959
Courant Number mean: 0.00177495 max: 4.40821
Courant Number mean: 0.00422436 max: 24.4105
Courant Number mean: 0.0173964 max: 47.9946
Courant Number mean: 0.0974709 max: 287.174
Courant Number mean: 0.679524 max: 1792.34
Courant Number mean: 39.8774 max: 147650
...........
Courant Number mean: 8.97372e+11 max: 6.71492e+15

看上去是每一个时间步的源项和上一个时间步的结果都在累加，从压力方程上每个时间步计算的过程来看的确也是这样的。

从物理过程的角度来说，连续的向模型中注入流体，每个时间步的压力方程都应该有源项存在的，如果源项大小不变，压力应该是稳态的，

不知道要如何避免这种压力不断增加的结果~,希望大家提提意见。

Tong

莫非只需要再第一个时间步添加源项，以后每个时间步只需要添加源项的改变量？

Foamer24

wellSource是常数吗？那mesh.V是什么参数？

东岳

@Tong 在非稳态计算向压力方程中添加源项中说：

非稳态计算过程中需要连续的向模型中注入流体，

具体是什么过程？

由于尺度的问题无法使用“显式”的入口，而是使用一个“隐式”的方法：添加一个体积源项到压力方程中实现注入；

虽然看的不是很懂但是感觉很厉害

但是，你的代码看起来并不“老练”，你压力方程的源项，出现在动量方程中，倒是可以理解。但为何不直接加在动量方程中？即使加在压力方程中，直接一减fvc::div(phiHbyA))-(wellSource/mesh.V())看起来很粗暴的样子，最好换成通量

话说回来，你这个应该是方程植入有问题。并且，你把源项放在压力方程中得到phi之后，速度需要进行减去或者加上

不是大问题，好好检查下你的方程和代码的一致性

yfclark

不知道你的源项具体是怎样的，我在压力方程中植入过源项，基本计算都不稳定，我的源项和相体积分数，压力温度相关，耦合十分严重，后来也是重新推导方程强行引入隐式的压力项求解才稳定的。另外reactingTwoPhaseEulerFoam的算例下有injection算例，通过fvOption植入的质量入注源项，你可以看一下

Tong

@Foamer24
wellSource是常数，mesh.V() 是网格体积，这样主要是为了和pEqn中的量纲保持一致性。

Tong

@东岳
可能是我表述的不太好，其实就是想在某些网格的连续方程上添加一个源项，来实现流体注入:

这个源项添加到压力方程上来实现不知道有没有问题(当然这样看上去比较粗暴，而且效率上不太好)，但是理解上比较直接，当然可以直接：

pEqn.source()+=wellSource;

或者使用fvOption ，不过感觉实际效果应该是一样的

Tong

@yfclark
多谢我去研究研究

Tong

回复一下大家，问题实锤了。源项添加应该是没有问题的。
问题出在：

fvm::laplacian(rTU, p)

中的rTU，因为在动量方程部分区域中添加了各向异性的达西阻力项，这一项是个对称张量场，如果换用各项同性的阻力场整个求解器就没有问题。应该计算达西阻力的时候出了点问题。

CFD小童

@Tong ，我现在在做凝固过程两相区内渗透率的各向异性，加载到动量中也就是以达西阻力项的方式添加，能请教下您关于添加各向异性达西阻力项的思路吗？

Tong

@CFD小童

嗯 OpenFoam 添加达西阻力源项有隐式和显式的两种，详细的可以参考
src\finiteVolume\cfdTools\general\porosityModel 目录下的程序

CFD小童

@Tong 好的，谢谢啊！