CFD Online Discussion Forums - [swak4Foam] viscous force calculation with MRF

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Community Contributions (https://www.cfd-online.com/Forums/openfoam-community-contributions/)

- - [swak4Foam] viscous force calculation with MRF (https://www.cfd-online.com/Forums/openfoam-community-contributions/149388-viscous-force-calculation-mrf.html)

viscous force calculation with MRF

Hi, I'm having trouble calculating viscous force over a patch with swak4Foam (actually two patches). For reference I'm using the built-in libforces calculation. This is the relevant swak4Foam function:

Code:

viscousForces

{

    type patchExpression;

    variables (

        "rho=998.3;"

        "viscous_force=-rho*nu*snGrad(U)*area();"

        "impeller1{patch'impeller1}=sum(-998.3*nu*snGrad(U)*area());"

    );

    accumulations (

        average

    );

    patches (

        impeller0

    );

    expression "sum(viscous_force)+impeller1";

    verbose true;

    outputControlMode   timeStep;

    outputInterval      1;

}

The implementation is copied from here https://openfoamwiki.net/index.php/C...ectionalForces. Now I suspect it might have something to do with the calculation being in MRF. How could I take that into account?

Here is some sample output from the functions:

Code:

Expression pressureForces on impeller0:  average=(4758.3887 -5556.7538 -51217.2)

Expression viscousForces on impeller0:  average=(-0.42171331 0.93947787 2.6970982)

forces rotorForces output:

    sum of forces:

        pressure : (4758.3887 -5556.7538 -51217.2)

        viscous  : (-0.62912771 1.5910869 5.3923825)

        porous   : (0 0 0)

    sum of moments:

        pressure : (-503.3495 -384.21045 -1364.6771)

        viscous  : (0.12301857 0.063222675 12.807648)

        porous   : (0 0 0)

The pressure forces seem to match nicely, but viscous forces do not. The reason I want to implement this using swak4Foam is that I want to automatically calculate some turbomachinery data, and moment calculations are needed in order to calculate power and efficiency.

Hi,

Are you running a turbulent case? If so, the standard forces functionObject takes the turbulent viscosity into account when calculating the viscous forces.

Regards,
Tom

Yes it's a turbulent case and I already tried using nu+nut in the calculations, still different results. Should it be just nut?

Output with "viscous_force=-rho*(nu+nut)*snGrad(U)*area();":

Code:

Expression pressureForces on impeller0:  average=(-2624.4874 1612.6975 33782.874)

Expression viscousForces on impeller0:  average=(-0.59524669 1.5039751 8.371256)

forces rotorForces output:

    sum of forces:

        pressure : (-2624.4874 1612.6975 33782.874)

        viscous  : (-0.64199103 1.4823397 8.9218996)

        porous   : (0 0 0)

    sum of moments:

        pressure : (143.58526 249.0284 559.52474)

        viscous  : (0.13726107 0.070810696 12.710413)

        porous   : (0 0 0)

The forces functionObject seems to use devRhoReff, which is defined in the turbulence model. Most turbulence models will use nuEff for this, which is indeed nut+nu.

The way the devRhoReff is calculated seems to be different to the way you have defined in the variables section, I think you should investigate on that to find the mismatch.

If you type this in your terminal:

Code:

cd $FOAM_SRC

grep -r devRhoReff .

You will get a list of files that use this function. Than it is up to you to find the relevant equation for the viscous forces calculation with your turbulence model.

Good luck,

Regards,
Tom

So it seems that for kOmegaSST the calculation is done by this function in kOmegaSST.C (incompressible):

Code:

tmp<volSymmTensorField> kOmegaSST::devReff() const

{

    return tmp<volSymmTensorField>

    (

        new volSymmTensorField

        (

            IOobject

            (

                "devRhoReff",

                runTime_.timeName(),

                mesh_,

                IOobject::NO_READ,

                IOobject::NO_WRITE

            ),

           -nuEff()*dev(twoSymm(fvc::grad(U_)))

        )

    );

}

I have no idea how to implement this with swak4Foam. Gradient operator grad(U) will not work because it's a patchexpression (?), snGrad(U) works, but it doesn't produce a tensor so twoSymm or dev won't work.