[user_of_cfx]

Hi,

I am trying to set up a simulation which can model a fluid which is:

- non-Newtonian liquid (specifically pseudoplastic fluid)
- compressible (I am interested in the pressure waves traveling up and down the domain)

From what I can see, there is no compressible non Newtonian single phase model within OF. The closest I can find is sonicLiquidFoam. The file structure looks pretty straightforward to me. Essentially it looks like the only thing I need to change is to make mu from a scalar to a volScalarField calcMu() -calculated by my viscosity model of choice- and maybe also modify the mu term of the UEqn?

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
);


As I am a complete novice in OF and C++, any advice, links, resources, tips, tricks and information of any kind will be greatly appreciated.

Many thanks,

Christa

[kebsiali]

Hello

I think the choice of the solver to use depends on the case you try to reproduce
a list of compressible Solvers available within openfoam is provided in this link
http://www.openfoam.org/features/standard-solvers.php

there are many viscosity models and these can be chosen from within the case folder, available in the following link

http://www.openfoam.org/features/transport.php

in case you want to create your own viscosity model or rheology law this has nothing to do with the solver itself
you need to create a new viscosity class in the libraries
this can be cleared out following this link
http://www.tfd.chalmers.se/~hani/kur...nFoam%20v2.pdf

[user_of_cfx]

Hi Ali,

Thank you very much for replying to me. I have done quite a bit of research and so far I can run a compressible Newtonian fluid case using rhoPimpleFoam and sonicLiquidFoam, and I have successfully implemented my own transport model (Bingham plastic) in nonNewtonianIcoFoam.

My problem is that all transport models are available for incompressible flows, but what I want to model is very clearly a compressible flow, and assuming incompressibility will just not work.

I believe sonicLiquidFoam is exactly the base solver I need for what I am trying to do, but I need some advice on how to implement the transport models, which are so far only available for incompressible solvers, into the compressible sonicLiquidFoam model.

To be a bit more precise, here is a list of the challenges I am trying to figure out:

1. How can I change mu from a scalar to a volScalarField for use with sonicLiquidFoam?
2. Are there any modifications I need to make to the transport model H and C files so it will run with compressible solvers?
3. How to link the transport model to the compressible solver. i.e. is it enough to include the right header files, or are there more modifications required?

I am currently reading through this document:
http://foam.sourceforge.net/docs/Gui...mmersGuide.pdf

but my limited knowledge of C++ makes this a bit hard for me.

Thanks,

Christa

Edit: I dug around and found this base class for compressible transport models:
https://github.com/OpenFOAM/OpenFOAM...s/compressible
Looks like a good place to start. Still, any help will be appreciated. Cheers, C

[user_of_cfx]

In case anyone wondered, it was quite easy to implement after all.

I modified sonicLiquidFoam.C to call a function fluid.mu() instead of a single scalar mu, with fluid being defined within createFields.H. The dynamic viscosity is calculated in a function calcMu(). I used the SinglePhaseTransportModel.H as guide to build my compessible version (which I imaginatively called compressibleSinglePhaseTransportModel.H) as well as the necessary transportModel and viscosityModel classes and individual vicosity models (I created a Bingham plastic for a start).

So long as one follows the right folder structure for the classes and the Make files and options, which can be found here, it is just a matter of making the relevant changes (such as using dynamic instead of kinematic viscosity) and modifying the units accordingly.

This is what the new solver code looks like:

#include "fvCFD.H"
#include "pimpleControl.H"
#include "compressibleSinglePhaseTransportModel.H" 

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    pimpleControl pimple(mesh);

    #include "readThermodynamicProperties.H"
  //"readTransportProperties.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readTimeControls.H"
        #include "compressibleCourantNo.H"

        fluid.correct();

        solve(fvm::ddt(rho) + fvc::div(phi));

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            fvVectorMatrix UEqn
            (
                fvm::ddt(rho, U)
              + fvm::div(phi, U)
              - fvm::laplacian(fluid.mu(), U)
            );

            solve(UEqn == -fvc::grad(p));

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                volScalarField rAU("rAU", 1.0/UEqn.A());
                surfaceScalarField rhorAUf
                (
                    "rhorAUf",
                    fvc::interpolate(rho*rAU)
                );

                U = rAU*UEqn.H();

                surfaceScalarField phid
                (
                    "phid",
                    psi
                   *(
                       (fvc::interpolate(U) & mesh.Sf())
                     + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                    )
                );

                phi = (rhoO/psi)*phid;

                fvScalarMatrix pEqn
                (
                    fvm::ddt(psi, p)
                  + fvc::div(phi)
                  + fvm::div(phid, p)
                  - fvm::laplacian(rhorAUf, p)
                );

                pEqn.solve();

                phi += pEqn.flux();

                solve(fvm::ddt(rho) + fvc::div(phi));
                #include "compressibleContinuityErrs.H"

                U -= rAU*fvc::grad(p);
                U.correctBoundaryConditions();
            }
        }

        rho = rhoO + psi*p;

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

After doing this I thought it is incredibly easy, but that's not how I felt when I started. So don't be intimidated, there is enough information in the standard OF solvers to use as basis for most modification you may want to make.

Cheers,

Christa

[StefanW]

Hi Christa,

good work you did here!

I am facing now a similar problem.
I want to implement nonNewtonian Fluids in a compressible multiphase solver (compressibleInterFoam) and so far my knowledge of C++ is still upgradable.

Therefore I want to ask if you would mind to publish the whole solver and openFoam parts you had to change.
So I could compare 1 by 1 to the original one you used and may understand what I have to change in the multiphase case.

Did you got any convergence probs with your solver?

It would be great if you could help!

Best Regards

Stefan