[ophr]

Hello people.


I am trying to implement a convection boundary condition by copying and modifying the mixedFvPatchField class. I chose this class due to it having almost all foundations necessary for it, as described by Alex Roqué in his foaming time wordpress website, in the downloads section.


All that remains to do is to calculate a value fraction and replace the original valueFraction_ private data member with it. This is due to the fact that, in the original mixedFvPatchField class, the value fraction is prescribed by the user. However, in the natural convection, it is a function of the convection coefficient in the patch, the heat conductivity coefficient in the cell and the face-to-center distance of each patch cell. One could calculate this fraction by hand in a simple mesh, regular mesh; but for an unstructured mesh, for instance, that would be impossible.


In order to calculate the new value fraction, I've removed the valueFraction_ variable from the class. In place, I've added two scalarField private data members kappa_ and h_, which are defined by the user in the boundary description. They refer to, respectively, the heat conductivity and convection coefficients. Then, they are used to calculate the value of the value fraction f, which replaces the valueFraction_ from the old formulation. Excerpts of code from the convectionFvPatchField.C, showing the dictionary reading constructor and the new evaluate function, are as follows:

template<class Type>
Foam::convectionFvPatchField<Type>::convectionFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const dictionary& dict
)
:
    fvPatchField<Type>(p, iF, dict, false),
    refValue_("refValue", dict, p.size()),
    refGrad_("refGradient", dict, p.size()),
    kappa_("kappa", dict, p.size()),
    h_("h", dict, p.size())
{
    evaluate();
}

template<class Type>
void Foam::convectionFvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
    if (!this->updated())
    {
        this->updateCoeffs();
    }

    scalarField beta = kappa_/h_ * this->patch().deltaCoeffs();
    scalarField f = 1.0/(1.0 + beta);

    WarningInFunction << "\nValue fraction: " << f
                      << "\nPatch-to-center distance: " << 1/this->patch().deltaCoeffs()
                      << "\nkappa: " << kappa_
                      << "\nh: " << h_
                      << "\nkappa_function: "<< this->kappa()
                      << "\nh_function: "<< this->h() << "\n" << endl;

    Field<Type>::operator=
    (

        f*refValue_
      +
        (1.0 - f)*
        (
            this->patchInternalField()
          + refGrad_/this->patch().deltaCoeffs()
        )
    );

    fvPatchField<Type>::evaluate();
}

The WarningInFunction line is just there for debugging purposes. When reading the log from my test case, it shows the following information for each patch:

--> FOAM Warning :
    From function void Foam::convectionFvPatchField<Type>::evaluate(Foam::UPstream::commsTypes) [with Type = double]
    in file convection/convectionFvPatchField.C at line 165

Value fraction: 1(0.0004997449)
Patch-to-center distance: 1(0.0005)
kappa: 1(9.0677014e-317)
h: 1(9.0676066e-317)

As you can see, both the kappa and h values aren't read and stored correctly by the class. They should read, respectively, 41.61 and 5. I have no idea why they are stored in that way.


As a test case, I am solving the plane wall with convection, as described in chapter 5.5 of Heat and Mass Transfer by Theodore Bergman, Frank Incropera et al. I am using the buoyantPimpleFoam solver with a fluid in which heat is transferred mostly by conduction (Prandtl = 1.7e-2). The initial condition at the domain is T = 280.15K and with a fluid at refValue = 301.15K.


Testing it with the convectionFvPatchField class I am writing, it shows a parabolic distribution of temperatures. That is correct. However, the concavity of the parabola should be facing upwards, and not downwards as in the test. A downwards facing concavity has the maximum temperature at the center of the domain, and not at its edge. Since the wall is being heated by the convection, the profile should be the exact opposite. This also suggests there is another unknown numerical error abound.


Can someone point me some directions? How can I have my new values be read correctly? Am I doing something else wrong apart from that?


Thanks for the patience,


Ophir