[aliyah.]

Hello everyone,

I am currently trying to implement the following custom boundary condition in OpenFOAM:






The boundary condition relates the gradient of a scalar field (phi) at a certain boundary (y=0) with the gradient of temperature. I have decided to use codedMixed for this purpose in the boundary field of phi.

Below is the code I've written, which seems to work up to a point:

fluid_to_solid
{
    type            codedMixed;
    value           uniform 5e-2;
    name            phitFlux;
    refValue        uniform 5e-2;       
    refGradient     uniform 0;      
    valueFraction   uniform 0;      //1 - Dirictlect ;  0 - Neumann

    codeInclude
#{
	#include "fvCFD.H"
#};

codeOptions
#{
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
  	-I$(LIB_SRC)/finiteVolume/lnInclude \
		-I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude 
#};

codeLibs
#{
	-lmeshTools \
	-lfiniteVolume
#};

code
#{
    // Get current patchi
    const label patchi = patch().index();

    // Accessing the temperature field of the current patch
    const volScalarField& TS=db().lookupObject<volScalarField>("T");

    // I am facing issues while trying to compute snGrad for the temperature field
    // When I use snGrad(), I encounter an error. However, removing .snGrad() (like shown below) lets it run without any problem.
    const scalarField& TG = TS.boundaryField()[patchi];  // How can I add .snGrad() here?

    // This is the phi of the current patch, which is not utilized in the code right now
    const scalarField& phip(*this);

    // Constants of the equation
    const scalar C1 = 1;
    const scalar C2 = 2;

    // Here, this->refGrad() represents the left side of the equation concerning phi.
    this->refGrad() = C1*C2*TG;  
}
#};

}

My specific issue arises when trying to calculate the snGrad of the temperature field. When I add .snGrad() to the line const scalarField& TG = TS.boundaryField()[patchi]; (I mean when I write const scalarField& TG = TS.boundaryField()[patchi].snGrad() I receive an error. Without .snGrad(), the code runs without any issues.

Could someone kindly provide insights or corrections on how I can successfully calculate the snGrad for the temperature field in the above context?

Thank you in advance!

[aliyah.]

I managed to add snGrad using the following code

// Access the entire temperature field
const volScalarField& T = db().lookupObject<volScalarField>("T");

// Calculate the surface-normal gradient for the entire temperature field
const surfaceScalarField gradT = fvc::snGrad(T);

// Extract the values at the boundary of interest
const scalarField& gradT_boundary = gradT.boundaryField()[patchi];

// Here, this->refGrad() represents the left side of the equation concerning phi.
this->refGrad() = C1*C2*gradT_boundary;

But there is a problem still when i am running the code, after the first iteration, it gives me core dumped error

#0  Foam::error::printStack(Foam::Ostream&) at ??:?                                                                                                                                                                                        #1  Foam::sigFpe::sigHandler(int) at ??:?                                                                                                                                                                                                  #2  ? in /lib/x86_64-linux-gnu/libpthread.so.0                                                                                                                                                                                             #3  ? in /lib/x86_64-linux-gnu/libm.so.6                                                                                                                                                                                                   #4  Foam::log(Foam::Field<double>&, Foam::UList<double> const&) at ??:?                                                                                                                                                                    #5  void Foam::log<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?                                    #6  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::log<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?                   #7  Foam::surfaceTensionModels::temperatureDependent::sigma() const at ??:?                                                                                                                                                                #8  Foam::phaseSystem::surfaceTensionCoeff(Foam::phasePairKey const&) const at ??:?                                                                                                                                                        #9  Foam::phaseSystem::surfaceTensionForce() const at ??:?                                                                                                                                                                                 #10  ? in ~/OpenFOAM/ubuntu-v2006/platforms/linux64GccDPInt32Opt/bin/icoReactingMultiphaseInterFoamPhi2BrounNew2                                                                                                                           #11  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6                                                                                                                                                                                  #12  ? in ~/OpenFOAM/ubuntu-v2006/platforms/linux64GccDPInt32Opt/bin/icoReactingMultiphaseInterFoamPhi2BrounNew2                                                                                                                           Floating point exception (core dumped)

[Lorenzo210]

Hi aliyah,


It seems that something could be wrong with the surface tension rather than the BC.
Could you please provide the phaseProperties dict of your test case?


I guessed the error is due the surface tension force because of the following lines:

#7 Foam::surfaceTensionModels::temperatureDependent:: sigma() const at ??:?

#8 Foam:haseSystem::surfaceTensionCoeff(Foam:hase PairKey const&) const at ??:?

#9 Foam:haseSystem::surfaceTensionForce() const at ??:?

Regards,
Lorenzo

[aliyah.]

Originally Posted by Lorenzo210

Hi aliyah,


It seems that something could be wrong with the surface tension rather than the BC.
Could you please provide the phaseProperties dict of your test case?


I guessed the error is due the surface tension force because of the following lines:





Regards,
Lorenzo

Hello Lorenzo,

Thanks for your tip. You were right. I was using a custom surface tension model and this custom surface tension model had logarithmic terms causing issues. When I used the abovementioned boundary condition, it resulted in negative values for the logarithm. I've now addressed this and things are working better.

Appreciate your help.

Best,