Point-by-point field – pressure( density ) free energy boundary condition in 3D

[araeli]

Hello everyone,

I have two questions regarding my current project. I’m working on implementing a capillary tube model at the micron scale using a free-energy model, and I would appreciate some clarification on the following points:

1) For the convergence criterion, I need to access the value of phi at each point, so I can compute the difference between vectors at the previous time step. However, I’m not entirely clear on how to access this data to copy it into a double vector in the general case. Any guidance on this would be greatly appreciated.

2) In defining the inlet and outlet, I would like to impose a density condition. I attempted to use coupling setParameter for this purpose, but it hasn’t been effective. Additionally, I’m unsure how to initialize the FreeEnergy InletOutletCoupling funtionalities, as it seems to be based only on velocity. Could anyone provide insights or examples on how to work with this? (The objective is to reproduce eqs 24-25 in https://doi.org/10.1002/2017WR021409)

Thank you in advance for your help!
Alice

[savis]

Hi Alice,

1) I’m not familiar with this, so someone might be able to offer a better way. Phi is accessed by calculating the ‘density’ on the order parameter super lattice (the one using FreeEnergyBGKdynamics). So, for each block lattice on this super lattice you should iterate through the cells and then call computeRho(), appending each value to the vector.

2a) Just to check, when you call setFreeEnergy[Inlet/Outlet]Boundary have you set the type argument to “density”, rather than “velocity”? You can then call defineRho to set the desired density on the inlet and outlet nodes. I think that should be sufficient.

2b) As for the inletOutletCoupling, this is needed to pass the velocity from the density lattice to the order parameter lattice(s) for the collision. This is especially necessary with density constraints because the velocity is constantly being updated. You can see how it is used in the microfluidics2d example. One thing to note is that, like ForceCoupling, the density lattice should be labelled B and the order parameter lattice labelled A (and C if using a ternary fluid system).

It has been a while since I have looked at this so there might be some inaccuracies.

Sam

[araeli]

Dear Sam,

1) Thank you for your suggestion regarding the convergence criterion. I will try this approach as I am still determining the best method for my specific physical problem.

2a-b) I set the density (and I believe the label is correct) using the following:

    setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice1,omega,superGeometry,3,"density",1);
    setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice1,omega,superGeometry,4,"density",1);

    setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice2,omega,superGeometry,3,"density",2);
    setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice2,omega,superGeometry,4,"density",2);

However, the boundary values are not being enforced as expected. For example, in a piston-like simulation, I tried setting an inlet density of 10 (non-physical, just for visualization purposes in Paraview), like so:

  AnalyticalConst3D<T,T> rhoinVal( inletDensity );
        AnalyticalConst3D<T, T> rhooutVal(T {1});
        AnalyticalConst3D<T,T> phiminusone( T {-1} );
        AnalyticalConst3D<T, T> phione(T {1});
        sLattice1.defineRho(superGeometry,3,rhoinVal);
        sLattice1.defineRho(superGeometry,3,rhooutVal);

        sLattice2.defineRho(superGeometry,3,phiminusone);//3 inlet material number
        sLattice2.defineRho(superGeometry,4,phione);// 4 outlet material number

        sLattice1.setProcessingContext<Array<momenta::FixedDensity::RHO>>(
                    ProcessingContext::Simulation);
        sLattice2.setProcessingContext<Array<momenta::FixedDensity::RHO>>(
                    ProcessingContext::Simulation);

I expected a rho value of 10 at the output, but it seems the boundary condition is not applied. I tried setting this before the coupling definitions, after them, and in the simulation loop, but without success.

I defined the couplings in the order: outlet, chemicalPotential, force, and inlet, as shown below:

    SuperLatticeCoupling coupling1(
                InletOutletCoupling3D{},
                names::A{}, sLattice2,
                names::B{}, sLattice1
                );

    coupling1.restrictTo(superGeometry.getMaterialIndicator({4}));

    SuperLatticeCoupling coupling2(
                ChemicalPotentialCoupling3D{},
                names::A{}, sLattice1,
                names::B{}, sLattice2);

    coupling2.template setParameter<ChemicalPotentialCoupling3D::ALPHA>(alpha);
    coupling2.template setParameter<ChemicalPotentialCoupling3D::KAPPA1>(kappa1);
    coupling2.template setParameter<ChemicalPotentialCoupling3D::KAPPA2>(kappa2);
    coupling2.restrictTo(superGeometry.getMaterialIndicator({1}));

    SuperLatticeCoupling coupling3(
                ForceCoupling3D{},
                names::A{}, sLattice2,
                names::B{}, sLattice1);

    coupling3.restrictTo(superGeometry.getMaterialIndicator({1}));

    ///inlet coupling
    SuperLatticeCoupling coupling4(
                InletOutletCoupling3D{},
                names::A{}, sLattice2,
                names::B{}, sLattice1
                );
    coupling4.restrictTo(superGeometry.getMaterialIndicator({3}));

Thank you for your suggestions. I hope the attached code provides some helpful context. I am sure that I am missing something fundamental, but I am struggling to identify exactly what it is.

Best regards,
Alice

[savis]

Hi Alice,

It looks like you may have made a mistake when defining rho on the inlet and outlet. You are using a material of ‘3’ for both, so the value that you set for the inlet gets overridden.

Regarding the couplings, you can combine coupling1 into coupling4 by applying it to both materials 3 and 4:
coupling4.restrictTo(superGeometry.getMaterialIndicator({3,4}));
This coupling should come after coupling3 because it needs the updated velocity.

Secondly, you should replace coupling1 with a DensityOutletCoupling. This ensures the density on the outlet remains fixed.

SuperLatticeCoupling coupling1(
    DensityOutletCoupling2D{},
    names::A{}, sLattice1,                                                                                                                                                           
    names::B{}, sLattice2
    );
coupling1.template setParameter<DensityOutletCoupling2D::RHO>(rhoOutVal);
coupling1.restrictTo(superGeometry.getMaterialIndicator({4}));

[Author]

Posts