CFD Online Discussion Forums - Cryogenics, compressible multiphase and heat transfer solver

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

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Cryogenics, compressible multiphase and heat transfer solver (https://www.cfd-online.com/Forums/openfoam-programming-development/97472-cryogenics-compressible-multiphase-heat-transfer-solver.html)

Cryogenics, compressible multiphase and heat transfer solver

Hello All,

I'm looking at solvers that can possibly do compressible cryogenic multiphase flow, and I want to get my head around multiphase flow first. The solver that looks like it might be the best starting point to modify if multiphaseEulerFoam. The description in the source code states that the solver takes into account compressibility and heat transfer, but a quick reading shows that there is no temperature/enthalpy solving and the models appear to be all incompressible. There is a TEqns.H commented out in the main solver file. The TEqn file itself looks like this:

Code:

{

    volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt));

    volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt));



    fvScalarMatrix T1Eqn

    (

        fvm::ddt(alpha1, T1)

      + fvm::div(alphaPhi1, T1)

      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1)

      - fvm::laplacian(kByCp1, T1)

     ==

        heatTransferCoeff*T2/Cp1/rho1

      - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)

      + alpha1*Dp1Dt/Cp1/rho1

    );



    fvScalarMatrix T2Eqn

    (

        fvm::ddt(alpha2, T2)

      + fvm::div(alphaPhi2, T2)

      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2)

      - fvm::laplacian(kByCp2, T2)

     ==

        heatTransferCoeff*T1/Cp2/rho2

      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)

      + alpha2*Dp2Dt/Cp2/rho2

    );



    T1Eqn.relax();

    T1Eqn.solve();



    T2Eqn.relax();

    T2Eqn.solve();



    // Update compressibilities

    psi1 = 1.0/(R1*T1);

    psi2 = 1.0/(R2*T2);

}

This looks like it was meant for twoPhaseEulerFoam, but never got implemented. With that introduction, I have the following questions/topics for discussion:

1) What experience/recommendations does anyone have as to the best solver to use as a starting point for this type of simulation?

2) From the TEqns.H above, it looks like there were some entities created for thermal conductivity and density, but the files that code was located in are no longer present in the current version of OF (2.1.x for me). Is there any way to see (or does anyone have) an old version of the code beside downloading old source packs and untaring them?

Marco,
You are correct. The heat transfer stuff in multiphaseEulerFoam are just holdovers from the previous twoPhase version and not fully implemented here. Not sure how many phases you will need--multiphaseEulerFoam is of course based on incompressible flow only, but is set up for an arbitrary number of phases. Presumably you want more than two else you would probably be looking to start from compressibleTwoPhaseEulerFoam instead...

Anyway, if you take a look at transportProperties in one of the example cases and phaseModel/phaseModel.H you will see that some of the tools for HX are already there--that is Cp and kappa are read in already.

I recently implemented temperature (on the mixture level) and buoyancy of the Buossinesq variety into this solver. This seemed to work pretty well for what I needed but it may not be what you need. Let me know if you have specific questions about how to do it.

Two phase would probably be enough to start, as I'm looking to model LNG tanks and pumps. Temperature would be useful, I would greatly appreciate having a look at your implementation as a jumping off point.

I will likely need to include a different thermophysical model (real gas behavior with phase change), but I'll leave that after I get some experience/insight into how these systems work.

Many thanks!

Well here is the TEqn.H that I used.

Code:

{

volScalarField rhoCp = fluid.rho()*fluid.Cp();



surfaceScalarField rhoPhi("rhoPhi", phi*1.0/fvc::interpolate(1.0/rho));



surfaceScalarField kByCpf = fvc::interpolate(fluid.kappa()/rhoCp + sgsModel->nut()/fluid.Prt());



 fvScalarMatrix TEqn 

 ( 

   fvm::ddt(T) 

 + fvm::div(phi, T)

 - fvm::laplacian(kByCpf, T) 

 ==

   sources(T)

 ); 



TEqn.relax();



sources.constrain(TEqn);



TEqn.solve();



}

You have to modify multiphaseSystem(.C,.H) to create extra functions which look exactly as the one for rho which spit out the mixture Cp and kappa. I have also included a turbulent Prandtl number but am not certain I have done this correctly. That needs to be added to phaseModel to be read in with the other transportProperties for each phase. Recall this is for multiphaseEulerFoam; it is also set up to use field sources as described here. As you said, it may make more sense for you to start from twoPhaseEulerFoam in which case the concept is the same.
Hope this is useful.