[mkraposhin]

slicedSurfaceScalarField works as pointer to data of array (field) that was specified during construction.

Originally Posted by strakey

Matvey - Thanks for the reply. I think I understand this better now. I downloaded the OF6 version and compiled it for OF7 with some minor changes. I also downloaded the Yeqn.H and heqn.H from the 2112 branch and this also compiled fine under OF7. The MULES modifications you made here did have a big impact on diffusion of temperature. I also have it outputting the hLambdaCoeff as a volScalarField with an fvc::average() procedure and the coefficients are near 1 now as expected. I still, however see too much diffusion in the species, much more than energy. From the Yeqn.H code, it looks like all species equations are being modified by the same set of lambdaCoeffs, which are 0 regardless of what I set my inertSpecies to. I'm not sure how a slicedSurfaceScalarField works though. I'm going to insert some more output fields to see what exactly is being computed for lambdaCoeffs for each of the species (4 in my case) but it seems as though lambdaCoeffs should be different for each species field?


I also tried backward, vanLeer and other things to minimize diffusion but these all have a small effect.


Lastly, I don't know if you have seen this work (https://www.appliedccm.com/wp-conten...M-2019-DWS.pdf), but this seems similar to what you are doing here with your hybrid solver.


Thanks
Pete

[mkraposhin]

Originally Posted by strakey

Matvey - Thanks for the reply. I think I understand this better now. I downloaded the OF6 version and compiled it for OF7 with some minor changes. I also downloaded the Yeqn.H and heqn.H from the 2112 branch and this also compiled fine under OF7. The MULES modifications you made here did have a big impact on diffusion of temperature. I also have it outputting the hLambdaCoeff as a volScalarField with an fvc::average() procedure and the coefficients are near 1 now as expected. I still, however see too much diffusion in the species, much more than energy. From the Yeqn.H code, it looks like all species equations are being modified by the same set of lambdaCoeffs, which are 0 regardless of what I set my inertSpecies to. I'm not sure how a slicedSurfaceScalarField works though. I'm going to insert some more output fields to see what exactly is being computed for lambdaCoeffs for each of the species (4 in my case) but it seems as though lambdaCoeffs should be different for each species field?



Update: It looks like the lambdaCoeffs (allFacesLambda) being used to modify the transport equations in Yeqn.H is that for the last species in the loop calling the MULES limiter,

forAll(Y, iCmpt)
    {
        volScalarField& Yi = Y[iCmpt];
        if ( Yi.name() != inertSpecie )
        {
            Info << "Limiting for: " << Yi.name() << endl;
            surfaceScalarField& rhoPhiYCorr = phiY[iCmpt];
        
            mulesWithDiffusionImplicitLimiter
            (
                rho,
                Yi,
                phi_own,
                phi_nei,
                allFacesLambda,
                rhoPhiYCorr,
                diffusiveFlux[iCmpt],
                mDCf[iCmpt],
                SuSp[iCmpt]
            );
        
        if (Yi.name() == "CH4")  //pete 2-28-2022
        {
        lambdavolH2O=fvc::average(lambdaCoeffs);
        }
        }
    }

it just keeps getting overwritten.



For my case, the lambdaCoeff for the last species is 0 everywhere, so 0 is being used for all of the species fields everywhere. I think hLambdaCoeff is just left at its initial value (1.0).



I also tried backward, vanLeer and other things to minimize diffusion but these all have a small effect.


Lastly, I don't know if you have seen this work (https://www.appliedccm.com/wp-conten...M-2019-DWS.pdf), but this seems similar to what you are doing here with your hybrid solver.


Thanks
Pete

Hi, Pete!


This is how MULES algorithm for multicomponent flow was designed. To make approximation of transport equations consistent between different phases, it takes minimum value of lambdaCoeff (the same approach is used for V-version of limited schemes in OpenFOAM, such as vanLeerV, MinmodV, etc):

\phi_f = \phi_f^{Upwind} + \lambda \times (\phi_f^{TVD} - \phi_f^{Upwind}) .

However, if \phi = 0, then MULES sets lambda to 0 and all other fields are transported with upwind. This is a very diffusive, but at the same time robust approach.





I think, this is rather technical issue, but it's solution requires some efforts.




Maybe you will find time to solve it

[strakey]

Matvey - I'm working on this but haven't found a good way to fix this problem without going into the MULESTemplates.C code and I was hoping to not have to mess with that.


Also, I think I found a bug in your Y.eqn (dev-2112) code on line 187:

owner = mesh.owner()[iFace];
neighbour = mesh.owner()[iFace];

shouldn't the second line be:

neighbour = mesh.neighbour()[iFace];

I made that change but then the code threw an "index out of range error" on line 232. Not sure what is happening there.


Thanks
Pete

[mkraposhin]

Yes, you are right. This is a bug, thank you.

To circumvent indexing problem, wrap next code into if (mesh.isInternalFace(iFace)):

    if (mesh.isInternalFace(iFace))
    {
            owner = mesh.owner()[iFace];
            neighbour = mesh.neighbour()[iFace];
            deltaY = mag(Y[iCmpt][owner] - Y[iCmpt][neighbour]);
            if (deltaY > maxDeltaY.primitiveField()[iFace])
            {
                 maxDeltaY.primitiveFieldRef()[iFace] = deltaY;
            }
    }

Originally Posted by strakey

Matvey - I'm working on this but haven't found a good way to fix this problem without going into the MULESTemplates.C code and I was hoping to not have to mess with that.


Also, I think I found a bug in your Y.eqn (dev-2112) code on line 187:

owner = mesh.owner()[iFace];
neighbour = mesh.owner()[iFace];

shouldn't the second line be:

neighbour = mesh.neighbour()[iFace];

I made that change but then the code threw an "index out of range error" on line 232. Not sure what is happening there.


Thanks
Pete

[strakey]

Thanks, that did the trick with the indexing error. Now I am seeing realistic values of hLambdaCoeffs and they are much better behaved than the lambdaCoeffs (only less than 1 in regions of steep phi gradients).


What I am testing now is moving the section of code: "Solve for Components" below the calculation of hLambdaCoeffs and using hLabmdaCoeffs on the Y equations (instead of lamdaCoeffs) as well as h. That way everything is using the same limiter and it seems to get around the phi=0 issue with MULES. I'll keep you posted.
Pete

[mkraposhin]

Please write if you get any progress

[Leckie]

Originally Posted by mkraposhin

C = sqrt(Cp/Cv/(drho/dp)) becomes sqrt(Cp/Cv/psi) = sqrt(Cp/Cv*p/rho) in case of perfect gas EOS. Where psi = (drho / dp) = 1/(R*T) and rho = psi*p
This formulation is true for all compressible implicit OF solvers. If you will change EOS, for example - rho = A*(p/p0)^k, then you have to change equation for pressure. for example - see sonicLiquidFoam

Hello, I know it's been some years since the last post on this thread. I'm still learning OpenFoam and I'm actually doing some work at Mach 6-7 based on the rhoCentralFoam/forwardStep tutorial and I'm having issues when applying different initial conditions aside of P=1Pa and T=1K... What changes do I need to make in order to simulate flow at Mach 6-7 with P=200Pa & T=220K?