[mnobrega]

Comments and questions are welcomed on the article Scheufler & Roenby, TwoPhaseFlow: A framework for developing two phase flow solvers in OpenFOAM, November 2023, OpenFOAM Journal, https://doi.org/10.51560/ofj.v3.80

[kerim]

Dear All,
https://doi.org/10.51560/ofj.v3.80
This DOI cannot be found in the DOI System.

[mnobrega]

Dear Kerim,
As mentioned in the initial message this paper is still in production. It will be published in the next few days.

[mnobrega]

Dear Kerim,

The paper was just published in the journal.


Good readings.

[kerim]

Thanks a lot!

[Graciella]

Hello, thanks for this wonderful solver framework.

[roenby]

A minor typo in the paper:

We refer a couple of places in the paper to the scheme pointCellLeastSquares. This should be pointCellsLeastSquares.

Cheers,
Johan

[aprakash]

Greetings dear authors

I am trying to leverage your compressible 2-phase solver for my research. I have a very naive question regarding equation 1. In the derivation of equation 1, it is my understanding that at some point we take \rho_k = \rho_k (p), i.e., express phase density as a function of pressure which in turn introduces the \psi_k -- phase compressibility terms. My question is that later on in equation 3, we also express rho_k that is phase density as a function of p and T. Why is T dependence not taken into account to derive the equation for \alpha i.e. equation 1 ?

* Is it because you assume pressure dominated compressibility or a weak T dependence of rho ?
* Is it because you want to simplify the transport equation for \alpha ?
* Or is it because you wish to disentangle thermal effects and treat them in the energy conservation equation ?

Is it fine to neglect T dependence of \rho in equation 1 for thermally driven phase changes ?

Maybe what is happening is that the \alpha (\partial . U) term in the RHS of equation (1) brings in the thermal influence because this term is also there in the energy conservation eqaution (6). I would deeply appreciate if you could clear my confusion.

[ancolli]

Hello,
Impressive work!
I'm currently encountering some issues when trying to use dynamicRefineFvMesh together with the reconstruction scheme isoSurface. For example, when running the damBreakWithObstacle case, I obtain the following:

Selected 0 split points out of a possible 563.
GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 5.27418e-09, No Iterations 18
GAMG:  Solving for pcorr, Initial residual = 0.00690608, Final residual = 3.53723e-09, No Iterations 12
time step continuity errors : sum local = 1.77161e-11, global = -5.035e-12, cumulative = -1.64867e-07
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  void Foam::volPointInterpolation::interpolateInternalField<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>&) const at ??:?
#4  void Foam::volPointInterpolation::interpolate<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>&) const at ??:?
#5  Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh> > Foam::volPointInterpolation::interpolate<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::word const&, bool) const at ??:?
#6  Foam::interpolationCellPoint<Foam::Vector<double> >::interpolationCellPoint(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#7  Foam::advection::isoAdvection::timeIntegratedFlux() at ??:?
#8  void Foam::advection::isoAdvection::advect<Foam::DimensionedField<double, Foam::volMesh>, Foam::DimensionedField<double, Foam::volMesh> >(Foam::DimensionedField<double, Foam::volMesh> const&, Foam::DimensionedField<double, Foam::volMesh> const&) in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
#9  ? in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
#10  ? in /lib/x86_64-linux-gnu/libc.so.6
#11  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#12  ? in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
Segmentation fault (core dumped)

Am I missing something, or do you know how to resolve this issue?
Thanks a lot!

[ancolli]

Originally Posted by ancolli

Hello,
Impressive work!
I'm currently encountering some issues when trying to use dynamicRefineFvMesh together with the reconstruction scheme isoSurface. For example, when running the damBreakWithObstacle case, I obtain the following:

Selected 0 split points out of a possible 563.
GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 5.27418e-09, No Iterations 18
GAMG:  Solving for pcorr, Initial residual = 0.00690608, Final residual = 3.53723e-09, No Iterations 12
time step continuity errors : sum local = 1.77161e-11, global = -5.035e-12, cumulative = -1.64867e-07
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  void Foam::volPointInterpolation::interpolateInternalField<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>&) const at ??:?
#4  void Foam::volPointInterpolation::interpolate<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>&) const at ??:?
#5  Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh> > Foam::volPointInterpolation::interpolate<Foam::Vector<double> >(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::word const&, bool) const at ??:?
#6  Foam::interpolationCellPoint<Foam::Vector<double> >::interpolationCellPoint(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#7  Foam::advection::isoAdvection::timeIntegratedFlux() at ??:?
#8  void Foam::advection::isoAdvection::advect<Foam::DimensionedField<double, Foam::volMesh>, Foam::DimensionedField<double, Foam::volMesh> >(Foam::DimensionedField<double, Foam::volMesh> const&, Foam::DimensionedField<double, Foam::volMesh> const&) in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
#9  ? in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
#10  ? in /lib/x86_64-linux-gnu/libc.so.6
#11  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#12  ? in ~/OpenFOAM/alejandro-v2206/platforms/linux64GccDPInt32Opt/bin/interFlow
Segmentation fault (core dumped)

Am I missing something, or do you know how to resolve this issue?
Thanks a lot!

I found the problem.
The correct member function in isoSurface.C should be:

// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

void Foam::reconstruction::isoSurface::reconstruct(bool forceUpdate)
{
    const bool uptodate = alreadyReconstructed(forceUpdate);

    if (uptodate && !forceUpdate)
    {
        return;
    }
    
    // Interpolating alpha1 cell centre values to mesh points (vertices)
    if (mesh_.topoChanging())
    {
        // Introduced resizing to cope with changing meshes
        if (ap_.size() != mesh_.nPoints())
        {
            ap_.resize(mesh_.nPoints());

        }
        if (interfaceCell_.size() != mesh_.nCells())
        {
            interfaceCell_.resize(mesh_.nCells());
        }
    }
    
    // Interpolating alpha1 cell centre values to mesh points (vertices)
    ap_ = volPointInterpolation::New(mesh_).interpolate(alpha1_);
    
    DynamicList<List<point>> facePts;

    interfaceLabels_.clear();

    forAll(alpha1_,cellI)
    {
            cutCell_.calcSubCell(cellI, iso_);

            if(cutCell_.cellStatus() == 0)
            {
                normal_[cellI] = cutCell_.faceArea();
                centre_[cellI] = cutCell_.faceCentre();
                if(mag(normal_[cellI]) != 0)
                {
                    interfaceCell_[cellI]=true;
                }
                else
                {
                    interfaceCell_[cellI]=false;
                    normal_[cellI] = vector::zero;
                    centre_[cellI] = vector::zero;
                }
            }
            else
            {
                normal_[cellI] = vector::zero;
                centre_[cellI] = vector::zero;
                interfaceCell_[cellI]=false;
            }

    }

}

I followed the definition in isoAlpha.C, which was working correctly.

[roenby]

Hi Alejandro


Can you provide a little more context? Was it a wrong setting in your case setup or did you find and fix a bug in the code? In the latter case, would you mind opening an issue here https://github.com/DLR-RY/TwoPhaseFlow and report the bug and fix so others can benefit from?


Cheers,
Johan

[ancolli]

Originally Posted by roenby

Hi Alejandro


Can you provide a little more context? Was it a wrong setting in your case setup or did you find and fix a bug in the code? In the latter case, would you mind opening an issue here https://github.com/DLR-RY/TwoPhaseFlow and report the bug and fix so others can benefit from?


Cheers,
Johan

Thanks for following up! It was indeed a bug (a small missing piece of code). I’ve already submitted a pull request with the minor modification needed. Thanks again!

[roenby]

Great, thanks! What is the pull request number?

[ancolli]

Sorry, I'm not sure what I did, but it's now #61.
I hope it's okay now!
Regards, and sorry again!

[aprakash]

In the TEqn.H for compressibleInterFlow, we have the Kinetic energy term

   
 fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
      - fvm::laplacian(turbulence.alphaEff(), T)
      + (
            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
           + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErr*K // KE term 
        )
       *(
           alpha1()/mixture.thermo1().Cv()()
         + alpha2()/mixture.thermo2().Cv()()
        )
     ==
        fvOptions(rho, T)
    );

I understand this and I can derive this as well. However I don't think a kinetic energy term has been included for the multiRegionPhaseChangeFlow solver's TEqn.H e,g, for the vapor phase (and the Liquid phase):

    fvScalarMatrix TVEqn
    (
        fvm::ddt(rhoCpV,TV)
        + fvm::div(rhoCpVPhi,TV)
	    - fvm::Sp(fvc::ddt(rhoCpV) + fvc::div(rhoCpVPhi), TV)
        - fvm::laplacian(kappaVEff,TV)
        - dpdt*alpha2
        + phaseChange.TSource2()
            ==
        fvOptions(rhoCpV, TV)
    );

Is there a reason why this was omitted ? I too can derive this expression but I have Kinetic energy terms also.
Thanks

[higorr]

I have recently started working with OpenFOAM and the multiRegionPhaseChangeFlow solver to simulate droplet evaporation. In my current setup, the droplet evaporation case is modeled as a 2D case with the frontAndBack patches set to empty, and the lateral boundaries defined as outlet patches with a fixed temperature boundary condition (fixedValue) of 500 K applied. The droplet did evaporate, but in a highly non-uniform way, forming strange and unexpected structures. I may be misunderstanding some fundamental aspect of heat transfer or solver behavior, and I would appreciate any guidance on this.

Additionally, I modified the stefanProblem testsuite by setting the air temperature to 500 K and applying a zeroGradient boundary condition on the initially heated wall. However, I observed that the heat flux flowed from the colder water region toward the hotter air region, which seems counterintuitive. I may be misunderstanding some fundamental aspect of heat transfer or solver behavior, and I would appreciate any clarification.

Finally, the solver uses three different temperature files in the setup, and I’m a bit confused about the purpose of each.

Thank you!