[LasseJacobsen]

Hi all,

In relation to my Master's Thesis on CHT modelling of fin-and-tube heat exchangers using Porous Media simplifications, I am implementing a custom solver to account for anisotropy in terms of the thermal diffusivity within the porous region.

My solver builds on chtMultiRegionSimpleFoam (OF v. 2206), and the following features are implemented:

• Porosity is defined as a volScalarField and read from a dictionary
• The diffusive term in the energy equation is modified by using a volTensorField as input for the anisotropic thermal diffusivity (in favor of the standard isotropic volScalarField representation)
• Interpolation between parallel and serial heat-transfer modes has been added for the PRINCIPAL directions of the thermal diffusivity tensor

For my own purposes, the thermal anisotropy is perfectly aligned with the coordinate system of my mesh.

However, for generalization purposes I would like for the user to be able to transform the 'principal diffusivity tensor' using an input transformation/rotation tensor.

In general, thermal diffusivity is represented as a 2nd order tensor - let's call it alpha. I want to define the principal directions of this tensor, and subsequently transform it using another 2nd order tranformation tensor - let's call it T. Mathematically, the transformed diffusivity tensor alpha' is computed as:



In OpenFOAM I have tried the following implementations:

• alpha' = transform(T, alpha)
• alpha' = T * alpha * T.T()

Either way I am getting what I think are type errors during compilation. But from what I can gather and should be types volTensorField and tensor, respectively - and hence be valid with the transform() function.

I will admit that I am indeed not an experienced C++ programmer - so any help will be much appreciated. Below I have added what I asses to be the required context - including all lines that have been added or changed in the chtMultiRegionSimpleFoam solver so far. Note that the current implementation is working, but no transforms are made during computation of the diffusivity tensor, as seen in the last code block. Note that the diffusivity tensor is named alpha_app in the code.

The inputs for interpolation between heat-transfer modes and the transformation tensor are defined like so:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      pmThermoProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Dimensioned scalar for solid phase thermal conductivity
kappa_s     kappa_s [1 1 -3 -1 0 0 0] 59;

// Non-dimensional input for interpolation between serial and parallel
// diffusive heat-transfer in the PM region for principal directions
// Input must be in range, 0 <= x <= 1, where:
// 1 = Perfectly parallel heat-transfer between solid-fluid phases
// 0 = Perfectly serial heat-transfer between solid-fluid phases
pmHeatTransferMode (1 1 0);

// Non-dimensional input for rotation tensor of anisotropic heat-transfer
// properties. Input the identity tensor if principal directions of heat-transfer
// modes correspond perfectly with the coordinate system of the mesh
rotationTensor (1 0 0 0 1 0 0 0 1);

// ************************************************************************* //

The ptrLists, since multiple regions must be supported, are initialized in createFluidFields.H:

PtrList<volScalarField> epsilonFluid(fluidRegions.size());              // Porosity (non-dimensional scalar field)
PtrList<IOdictionary> pmThermoProperties(fluidRegions.size());          // PM thermal property dictionary
PtrList<dimensionedScalar> kappaSolid(fluidRegions.size());             // Solid PM phase thermal conductivity (dimensioned scalar)                    
PtrList<vector> pmHeatTransferMode(fluidRegions.size());                // PM diffusive heat-transfer mode in principal directions (non-dimensioned vector, 0<=x<=1)
PtrList<tensor> rotationTensor(fluidRegions.size());                    // Rotation tensor for principal thermal diffusivity tensor
PtrList<volTensorField> alphaApplied(fluidRegions.size());              // Anisotropic thermal diffusivity (dimensional tensor field)
PtrList<volScalarField> alphaSerial(fluidRegions.size());               // Thermal diffusivity for serial heat-transfer (dimensional scalar field)             
PtrList<volScalarField> alphaParallel(fluidRegions.size());             // Thermal diffusivity for parallel heat-transfer (dimensional scalar field)

The various new fields and scalars are read and initialized in the forAll(fluidRegions, i) like so:

    // Read dictionary for porosity term (specified with setCells for cellZones)
    Info<< "    Adding to epsilonFluid\n" << endl;
    epsilonFluid.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "porosity",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            fluidRegions[i]
        )
    );
    
    // Add PM thermal properties dictionary
    Info<< "    Adding to pmThermoProperties\n" << endl;
    pmThermoProperties.set
    (
        i,
        new IOdictionary
        (
             IOobject
             (
                    "pmThermoProperties",
                    runTime.constant(),
                    fluidRegions[i],
                    IOobject::MUST_READ_IF_MODIFIED,
                    IOobject::NO_WRITE
             )
        )
    );

    // Read solid phase thermal conductivity from PM thermo dict
    kappaSolid.set
    (
        i,
        new dimensionedScalar(pmThermoProperties[i].lookup("kappa_s"))
    );
    // Read heat-transfer mode vector
    pmHeatTransferMode.set
    (
        i,
        new vector(pmThermoProperties[i].lookup("pmHeatTransferMode"))
    );
    // Read rotation tensor from PM thermo dict
    rotationTensor.set
    (
        i,
        new tensor(pmThermoProperties[i].lookup("rotationTensor"))
    );


    // Add tensor field for applied thermal diffusivity
    Info<< "    Adding to alphaApplied\n" << endl;
    alphaApplied.set
    (
        i,
        new volTensorField
        (
            IOobject
            (
                "alphaApplied",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::AUTO_WRITE
            ),
            fluidRegions[i]
        )
    );
    alphaSerial.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "alphaSerial",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE
            ),
            fluidRegions[i],
            dimensionedScalar("alphaSerial", dimensionSet(1,-1,-1,0,0,0,0), 0)
        )
    );
    alphaParallel.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "alphaParallel",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE
            ),
            fluidRegions[i],
            dimensionedScalar("alphaParallel", dimensionSet(1,-1,-1,0,0,0,0), 0)
        )
    );

In setRegionFluidFields.H I define:

    // Extract values from added pointer lists
    const volScalarField& epsilon = epsilonFluid[i];                            // Porosity field
    const dimensionedScalar& kappa_s = kappaSolid[i];                           // Solid phase thermal conductivity
    volTensorField& alpha_app = alphaApplied[i];                                // Applied thermal diffusivity
    volScalarField& alpha_serial = alphaSerial[i];                              // Serial thermal diffusivity
    volScalarField& alpha_parallel = alphaParallel[i];                          // Serial thermal diffusivity
    const vector& heatTransferMode = pmHeatTransferMode[i];                     // Serial-parallel heat-transfer mode scaling
    const tensor& rotTensor = rotationTensor[i];                                // Rotation tensor for principal thermal diffusivity tensor

Lastly, the anisotropic thermal diffusivity tensor is computed, and utilized in the energy equation in EEqn.H:

    // Compute thermal diffusivity
    // alpha_eff = rho_f * alpha_f
    // alpha_applied = alpha_i / c_pf
    // kappa_XY = epsilon * (alpha * cpf) + (1 - epsilon) * kappa_s
    // kappa_Z = (alpha * cpf * kappa_s) / (epsilon * kappa_s + (1 - epsilon * alpha * cpf))

    // Compute perfectly serial and parallel heat-transfer mode thermal diffusivities
    alpha_serial = (((turb.alphaEff() * thermo.Cp() * kappa_s) 
                    / (epsilon * kappa_s + (1.0 - epsilon) * turb.alphaEff() * thermo.Cp())) 
                    / thermo.Cp());
    alpha_parallel = ((epsilon * (turb.alphaEff() * thermo.Cp()) 
                    + (1.0 - epsilon) * kappa_s) / thermo.Cp());

    // Interpolate between heat-transfer modes and setup principal thermal diffusivity tensor
    tensor i (1, 0, 0, 0, 0, 0, 0, 0, 0);
    tensor j (0, 0, 0, 0, 1, 0, 0, 0, 0);
    tensor k (0, 0, 0, 0, 0, 0, 0, 0, 1); 
    alpha_app = i * (heatTransferMode.x() * alpha_parallel + (1 - heatTransferMode.x()) * alpha_serial)
              + j * (heatTransferMode.y() * alpha_parallel + (1 - heatTransferMode.y()) * alpha_serial)
              + k * (heatTransferMode.z() * alpha_parallel + (1 - heatTransferMode.z()) * alpha_serial);

    fvScalarMatrix EEqn
    (
      fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(alpha_app, he)
     ==
        rho*(U&g)
      + rad.Sh(thermo, he)
      + fvOptions(rho, he)
    );

[Mike2018]

Hi Lasse,

Did you find a solution to this?

I am modelling something very similar using OpenFOAM 2212.

It is an airfin printed circuit heat exchanger adapting the CHTmultregionsimplefoam, with a porous media approach.

I have regions with bends, and the solid thermal conductivity requires a different coordinate system.

Note that for the friction factor application in fvOptions, you can set different coordinate systems (i.e. cylindrical). However, it is not so trivial for the conductivity (within thermophysical properties.)

Please reach out with progress and we can help each other.

Best,
Mike

[LasseJacobsen]

Hi Mike,

In short I never did manage to solve the initially defined problem, that to my recollection in short revolved around defining the principal directions of porous-media conductivity and subsequently rotating the conductivity tensor by a user-defined input.

For this problem though, I am sure some trivial solution exists that is just beyond my OpenFOAM programming skills at least.

However, the problem you are describing seems to be a bit more complex if I understand it correctly. From what I gather you may actually have multiple referenceframes (RFs) in the same domain for which you'd need to do different transformations or?

E.g., see image below where there is a straight section, requiring some (or maybe no transformation to an input principal conductivity tensor), whereas the curvy section would need another transformation.

If that is the case, likely you'll need to somehow define the conductivity tensor for the porous-media domain of interest as a function of space (which I unfortunately cannot give any great hints on how to succesfully achieve).

ApplicationFrameHost_PPAVaLC2po.png

[Mike2018]

Hi Lasse,

Thank you for your response.

I too believe there may be a solution to this but it is outside my ability as of yet.

After further investigation into my problem, it seems the conductivity anisotropy can be neglected as the heat transfer is dominated by the volume heat source.

You are correct in your assessment of the problem. You can set a global coordinate system for the conductivity (kappa), and you can also set anisotropy with

thermoType
{
...
transport constAnIso;
...
}

mixture
{
...
transport
{
kappa ($KAPX $KAPY $KAPZ);
}
...
}

coordinateSystem
{
x (1 0 0);
y (0 1 0);
#includeEtc "caseDicts/general/coordinateSystem/cartesianXY"
}

(or maybe use cylindrical coordinates)
..........

However, there is no obvious way to set this field in specific cell zones with a given coordinate system. With fvOptions for DarcyForchiemer you can access these features for the momentum source term.

There might be some way to utilise setFields and manually correct the kappa field and make it a function of space as you suggested.

I left these suggestions here for any researchers looking in

Have a nice day!