#201806
在 interFoam 的 pEqn.H 里面，有这样一段：
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
这一段程序是出现在 pEqn.H 里面求解 p 方程，并进行通量非正交修正后，对 U 变量实行的一种修正方式。
那么为什么要这么做？
相信大家对下面这个方程也不再陌生了
这个方程是利用 U 的预测值和 p 的求解值来获得 U 的真实值，在上一期 code of the week 我们提到了通量的修正。但是，当通量修正后，速度并没有自动修正，而要重新处理。
而速度处理有两种方式，一种是直接利用上面这个式子获得，在 incompressible 的一些求解器里面，使用的是这样一种方式，比如 icoFoam：
U -= rUA*fvc::grad(p);
那开头那一种方式是怎么回事呢？
进行修正了的通量是 phi，未进行修正的通量是 phiU，他们的差值，自然就是
式中， \dfrac{1}{A_p} 就是 rAUf。当前， phi 和 phiU 都已求得，那如何从守恒的通量 \text{phi} 获得守恒的 U 呢？
我们来看看 reconstruct 都干了啥？ 文件是 fvcReconstruct.C，
template<class Type>
tmp
<
GeometricField
<
typename outerProduct<vector,Type>::type, fvPatchField, volMesh
>
>
reconstruct
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
typedef typename outerProduct<vector, Type>::type GradType;
const fvMesh& mesh = ssf.mesh();
tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
(
new GeometricField<GradType, fvPatchField, volMesh>
(
IOobject
(
"volIntegrate(" + ssf.name() + ')',
ssf.instance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
& surfaceSum((mesh.Sf()/mesh.magSf())*ssf),
zeroGradientFvPatchField<GradType>::typeName
)
);
GeometricField<GradType, fvPatchField, volMesh>& reconField =
treconField();
// Note:
// 1) Reconstruction is only available in cell centres: there is no need
// to invert the tensor on the boundary
// 2) For boundaries, the only reconstructed data is the flux times
// normal. Based on this guess, boundary conditions can adjust
// patch values
// HJ, 12/Aug/2011
GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
surfaceSum((mesh.Sf()/mesh.magSf())*ssf);
reconField.internalField() =
(
inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField())
& fluxTimesNormal.internalField()
);
reconField.boundaryField() = fluxTimesNormal.boundaryField();
treconField().correctBoundaryConditions();
return treconField;
}
其中，核心的一段是
inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
& surfaceSum((mesh.Sf()/mesh.magSf())*ssf),
这一段直观地、返回的量是
为什么可以这样？设一对称张量 \mathcal D
有
上式的意思可以表述为：对称张量将速度向量 \boldsymbol U 映射到与 \boldsymbol n_f 平行的向量上，其幅值为 |\boldsymbol n_f| \left( \boldsymbol S_f \cdot \boldsymbol U \right) \approx |\text{phi}| ，即 \boldsymbol n_f \text{phi} 。OpenFOAM 认为式 (7) 是成立的，因此，运用了这一方法来重构速度。
有些人对此两种方法进行了对比，认为差别并不大。Henry Weller 也表示过 fvc::reconstruct() 可能会导致较大的误差。openfoamwiki 认为使用这个的原因是消除棋盘压力[2]。因此，暂时未有定论。对于哪一种方法更好，可能是 depends on cases。