CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   Cryogenics, compressible multiphase and heat transfer solver (

mturcios777 February 17, 2012 15:19

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:


    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)
      - 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)
      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
      + alpha2*Dp2Dt/Cp2/rho2



    // 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?

kwardle March 16, 2012 16:17

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.

mturcios777 March 16, 2012 16:31

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!

kwardle March 22, 2012 11:11

Well here is the TEqn.H that I used.


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::div(phi, T)
 - fvm::laplacian(kByCpf, T)





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.

All times are GMT -4. The time now is 07:02.