[TomasDenk]

Hello,

I would like to solve radiative heat transfer in two adjacent regions with chtMultiRegionFoam. I was able to mimic the solution procedure already present for fluid regions and include it in solid regions. However, I can't figure out how to implement boundary condition to couple radiation intensity.

Assumptions:
-vfDOM radiation model (we can also assume the same azimuthal and polar angles)
-fixed solid/fluid or solid/solid interface (constant shape)
-equal refractive index, i.e. no need to redistribute ray intensity that passes through the interface
-grey radiation

In the greyDiffusiveRadiationMixedFvPatchScalarField.C (OpenFOAM 3.0.x) the essential part of code:

forAll(Iw, faceI)
    {
        if ((-n[faceI] & myRayId) > 0.0)
        {
            // direction out of the wall
            refGrad()[faceI] = 0.0;
            valueFraction()[faceI] = 1.0;
            refValue()[faceI] =
                (
                    Ir[faceI]*(scalar(1.0) - temissivity[faceI])
                  + temissivity[faceI]*physicoChemical::sigma.value()
                  * pow4(Tp[faceI])
                )/pi;

            // Emmited heat flux from this ray direction
            Qem[faceI] = refValue()[faceI]*nAve[faceI];
        }
        else
        {
            // direction into the wall
            valueFraction()[faceI] = 0.0;
            refGrad()[faceI] = 0.0;
            refValue()[faceI] = 0.0; //not used

            // Incident heat flux on this ray direction
            Qin[faceI] = Iw[faceI]*nAve[faceI];
        }
    }

How should I treat the direction out of the interface? I think I should read intensities I_?? of particular rays in the other region and set it to refValue()[faceI] instead of reflected intensity and black body radiation. How can I implement it?

[TomasDenk]

Looking at some other coupling BCs, I've managed to get the first approximation of desired boundary condition. Now, my code looks like this:

// Get the coupling information from the mappedPatchBase
    const mappedPatchBase& mpp =
        refCast<const mappedPatchBase>(patch().patch());
    const polyMesh& nbrMesh = mpp.sampleMesh();
    const label samplePatchI = mpp.samplePolyPatch().index();
    const fvPatch& nbrPatch =
        refCast<const fvMesh>(nbrMesh).boundary()[samplePatchI];
    const greyRadiationIntensityCouplingMixedFvPatchScalarField& nbrField =
        refCast <const greyRadiationIntensityCouplingMixedFvPatchScalarField>
        (nbrPatch.lookupPatchField<volScalarField, scalar> (dimensionedInternalField().name()) );
    tmp<scalarField> nbrIw(new scalarField(nbrField.size(), 0.0));
    nbrIw() = nbrField.patchInternalField();
    mpp.distribute(nbrIw());
    
    const radiationModel& radiation =
        db().lookupObject<radiationModel>("radiationProperties");

    const fvDOM& dom(refCast<const fvDOM>(radiation));

    label rayId = -1;
    label lambdaId = -1;
    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);

    const label patchI = patch().index();

    if (dom.nLambda() != 1)
    {
        FatalErrorIn
        (
            "Foam::radiation::"
            "greyRadiationIntensityCouplingMixedFvPatchScalarField::updateCoeffs"
        )   << " a grey boundary condition is used with a non-grey "
            << "absorption model" << nl << exit(FatalError);
    }

    scalarField& Iw = *this;
    const vectorField n(patch().nf());

    radiativeIntensityRay& ray =
        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));

    const scalarField nAve(n & ray.dAve());

    ray.Qr().boundaryField()[patchI] += Iw*nAve;

    scalarField& Qem = ray.Qem().boundaryField()[patchI];
    scalarField& Qin = ray.Qin().boundaryField()[patchI];

    const vector& myRayId = dom.IRay(rayId).d();

    forAll(Iw, faceI)
    {
        if ((-n[faceI] & myRayId) > 0.0)
        {
            // direction out of the wall
            refGrad()[faceI] = 0.0;
            valueFraction()[faceI] = 1.0;
            // assumes the same rays, i.e. phi and theta angles, for both sides
            refValue()[faceI] = nbrIw()[faceI];

            // Emmited heat flux from this ray direction
            Qem[faceI] = refValue()[faceI]*nAve[faceI];
        }
        else
        {
            // direction into the wall
            valueFraction()[faceI] = 0.0;
            refGrad()[faceI] = 0.0;
            refValue()[faceI] = 0.0; //not used

            // Incident heat flux on this ray direction
            Qin[faceI] = Iw[faceI]*nAve[faceI];
        }
    }

It uses intensity from coupled region as refValue for every ray. However, this works only for the same setting of rays in both radiation regions.

Can someone post a code that would retrieve radiation model from neighboring region (on the other side of the interface) to allow me to test for nPhi and nTheta?

Would it make sense to use Qin from the other side and distribute intensities uniformly on my side when nbr_nPhi <> nPhi or nbr_nTheta <> nTheta?