reactingTwoPhaseEulerFoam计算面上的力

王长宏

各位老师好，我现在在做船舶气泡减阻的数值模拟，在使用reactingTwoPhaseEulerFoam求解器计算船体阻力的时候出现问题.
2206版本的force.C里的部分代码如下：

Foam::tmp<Foam::volSymmTensorField>
Foam::functionObjects::forces::devRhoReff() const
{
    typedef compressible::turbulenceModel cmpTurbModel;
    typedef incompressible::turbulenceModel icoTurbModel;

    if (foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
    {
        const auto& turb =
            lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);

        return turb.devRhoReff();
    }
    else if (foundObject<icoTurbModel>(icoTurbModel::propertiesName))
    {
        const auto& turb =
            lookupObject<icoTurbModel>(icoTurbModel::propertiesName);

        return rho()*turb.devReff();
    }
    else if (foundObject<fluidThermo>(fluidThermo::dictName))
    {
        const auto& thermo = lookupObject<fluidThermo>(fluidThermo::dictName);

        const auto& U = lookupObject<volVectorField>(UName_);

        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
    }

在欧拉欧拉方法的两相流求解器的求解过程中会出现U.water,U.air等场，但在force.C里计算牛顿应力张量的部分只识别U,导致无法计算粘性力。

Force = fp + fv;
Fp = Sfb * p. boundaryField();
Fv = Sfb & devRhoReffb;
devRhoReffb = devRhoReff().boundaryField();
devRhoReff = rho * devReff;
devReff = nuEff * grad(U);
nuEff = nut +nu;
nu = mu / rho;
mu = alpha1 * rho1 * nu1 + alpha2 * rho2 * nu2;
rho = alpha1 * rho1 + alpha2 * rho2;
rho1 = 1000 kg/m^3;
rho2 = 1 kg/m^3;
nu1 = 1 * 10^(-6) m^2/s;
nu2 = 1.48 * 10^(-5) m^2/s;
nut = Cu * k^2 / epsilon, Cu = 0.09；

后续仿照interFoam求解粘性力的部分 rho()*turb.devReff()把force.C里的替换掉了，但忽略了nut.air等量，目前的计算结果显示reactingTwoPhaseEulerFoam的受力比interFoam略小，想问下各位老师按interFoam计算有效粘度和应力张量是否是合理的？十分感谢！

李东岳

TFM这面的受力你得计算方法是正确的么。这部分据我所知就是因为两相流的受力计算方法口径不统一，因此默认没有带力的计算方法。

王长宏

@李东岳谢谢李老师，interFoam里粘性应力计算我是按上面那几行公式算的，所有的场都输出来验证过，最后得到的devRhoReff和rho()*turb.devReff是一样的；受力的计算方法不统一的话我试试别的公式，找个差别小的就用了。

CFD中文网

reactingTwoPhaseEulerFoam计算面上的力