[MMRC]

Good day dear community,



I want to model the effect of a bunch of trees in atmospheric flow. I know there some tools already in openfoam ESI for canopy but I want to follow strictly Dalpe & Masson theory on this.



I have only elaborated a class for the source terms on turbulent variables: k and epsilon & omega. The files are attached below. These variable grows in value in this type of physics, so I assumed by '-=' it is a source term, and not a sink. Either way changing + to - or viceversa does not help when simulating.



The idea is to apply these sourceterms to a cellSet via a fvOption tool.



They compile without issues (wclean .all + wmake libso) in my custom user library for openfoam ESI V2012. The issue happens when running simpleFoam in OpenFOAM ESI 2012 in a singlecore, terminal prompts this error:

mmros@mmros:~/OpenFOAM/mmros-v2306/run/simulations/forestPractices/p1-TurbSourceOnly$ simpleFoam
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : simpleFoam
Date   : Jul 03 2024
Time   : 12:46:03
Host   : mmros
PID    : 230848
I/O    : uncollated
Case   : /home/mmros/OpenFOAM/mmros-v2306/run/simulations/forestPractices/p1-TurbSourceOnly
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: convergence criteria
    field p     tolerance 1e-08
    field U     tolerance 1e-08
    field "(k|epsilon)"     tolerance 1e-12

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
RAS
{
    RASModel        kEpsilon;
    turbulence      on;
    printCoeffs     on;
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              0;
    sigmak          1;
    sigmaEps        1.3;
}

No MRF models present

Creating finite volume options from "constant/fvOptions"

Selecting finite volume options type treeCanopySource
    Source: forestTurbSource
    - selecting cells using cellSet wdCellSet
    - selected 160 cell(s) with volume 31277.4802
    Applying treeCanopySource to: epsilon and k with Turb. model constants of:  Cmu:  0.09C1: 1.43C2: 1.92sigmaE: 1.3sigmaK: 1 The canopy values of Cd: 0.8 and LAD: 2 with betaP: 1

Starting time loop

Time = 1

DILUPBiCGStab:  Solving for Ux, Initial residual = 0.007754258079, Final residual = 1.248271136e-05, No Iterations 3
DILUPBiCGStab:  Solving for Uy, Initial residual = 1, Final residual = 2.908675916e-05, No Iterations 3
DILUPBiCGStab:  Solving for Uz, Initial residual = 1, Final residual = 0.001796598423, No Iterations 3
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.04614133528, No Iterations 30
GAMG:  Solving for p, Initial residual = 0.0002426027944, Final residual = 1.159132676e-05, No Iterations 32
GAMG:  Solving for p, Initial residual = 1.156102622e-05, Final residual = 8.983712988e-07, No Iterations 28
time step continuity errors : sum local = 8.612708198e-08, global = 9.50297687e-09, cumulative = 9.50297687e-09
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::tmp<Foam::fvMatrix<double> > Foam::fv::optionList::source<double>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&, Foam::dimensionSet const&) at ??:?
#4  Foam::tmp<Foam::fvMatrix<double> > Foam::fv::optionList::operator()<double>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&) at ??:?
#5  Foam::RASModels::kEpsilon<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::correct() at ??:?
#6  ? in ~/OpenFOAMv2012/OpenFOAM-v2012/platforms/linux64GccDPInt32Opt/bin/simpleFoam
#7  ? in /lib/x86_64-linux-gnu/libc.so.6
#8  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#9  ? in ~/OpenFOAMv2012/OpenFOAM-v2012/platforms/linux64GccDPInt32Opt/bin/simpleFoam
Segmentation fault (core dumped)

I don't understand the error that is prompting nor what fix does it suggest, could some one take a look on this and explain better, please?

I have attached the files of the class with the make files, + the fvOption file to ease the test.



I bet the issue maybe in my source terms but I'm not quite sure how to make them implicit if it's the case.



Big thanks to anyone who can take a look on this,

[MMRC]

Hello community,



I think my answer is in this tutorial that covers the same topic and clarifies how to define implicit terms: https://www.tfd.chalmers.se/~hani/ku..._Segersson.pdf

[MMRC]

Hello community,


I changed the way of handling the topic, instead of a loop inside a cellSet, which was quite tricky to find a way to apply implicit terms in fvMatrix form into a specific group of cells, I will rely now on on having scalar fields of my Cd and LAD mapped in the mesh, so whether these values are non zero, then a sink/source term will appear in the flow.



The tutorial showed that addSup could be overload so in one class we can code source
term for a vector (velocity) and a scalar (k,epsilon...etc)


I can not compile my new class implementation. New errors are related to:
1) In TKE:

In file included from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.H:329,
                 from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.C:28:
../fvOptions/sources/derived/treeCanopySource/treeCanopySourceTemplates.C:89:67: note:   mismatched types ‘Foam::DimensionedField<Type, GeoMesh>’ and ‘Foam::Field<double>’
   89 |         alpha()*rho()*calcPlantCanopyTerm(U)*betaP_*pow(mag(U),2) - fvm::Sp(betaD_*calcPlantCanopyTerm(U),k)
      |         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedField.C:590,
                 from /home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedField.H:420,
                 from /home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/finiteVolume/lnInclude/fvMesh.H:61,
                 from /home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/fvOptions/lnInclude/cellSetOption.H:104,
                 from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.H:122,
                 from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.C:28:
/home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedFieldFunctions.C:370:1: note: candidate: ‘template<class Type1, class Type2, class GeoMesh> Foam::tmp<Foam::DimensionedField<typename Foam::typeOfSum<arg1, arg2>::type, GeoMesh> > Foam::operator-(const Foam::DimensionedField<TypeR, GeoMesh>&, const Foam::DimensionedField<Type1, GeoMesh>&)’
  370 | operator op                                                                    \
      | ^~~~~~~~
/home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedFieldFunctions.C:643:1: note: in expansion of macro ‘PRODUCT_OPERATOR’
  643 | PRODUCT_OPERATOR(typeOfSum, -, subtract)
      | ^~~~~~~~~~~~~~~~
/home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedFieldFunctions.C:370:1: note:   template argument deduction/substitution failed:
  370 | operator op                                                                    \
      | ^~~~~~~~
/home/mmros/OpenFOAMv2012/OpenFOAM-v2012/src/OpenFOAM/lnInclude/DimensionedFieldFunctions.C:643:1: note: in expansion of macro ‘PRODUCT_OPERATOR’
  643 | PRODUCT_OPERATOR(typeOfSum, -, subtract)

And for epsilon, there is this alert:

In file included from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.H:329,
                 from ../fvOptions/sources/derived/treeCanopySource/treeCanopySource.C:28:
../fvOptions/sources/derived/treeCanopySource/treeCanopySourceTemplates.C: In instantiation of ‘void Foam::fv::treeCanopySource::treeCanopySourceEpsilon(const AlphaFieldType&, const RhoFieldType&, Foam::fvMatrix<double>&, Foam::label) [with AlphaFieldType = Foam::geometricOneField; RhoFieldType = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>; Foam::label = int]’:
../fvOptions/sources/derived/treeCanopySource/treeCanopySource.C:293:32:   required from here
../fvOptions/sources/derived/treeCanopySource/treeCanopySourceTemplates.C:63:12: error: no matching function for call to ‘Sp(Foam::tmp<Foam::Field<double> >, const volScalarField&)’
   63 |     fvm::Sp(rho() * alpha() * calcPlantCanopyTerm(U)* (Ce4 * betaP_ * pow(mag(U), 2) - Ce4 * betaD_ * k), epsilon)

I followed as reference the previous tutorial and the current fvOption atmPlantCanopyTurbSource code.



I attached a new compressed folder where the source files + fvOption dict + boundary dicts are included to ease the test


The idea is to mostly run this sourceterm with simpleFoam, and sometimes maybe with bouyantBoussinesqSimpleFoam


Big tanks to anyone that can take a look in this.