[sturgeon]

Hi all. I have done a lot of research on this but can't find anything to solve my problem, but I am pretty new to OpenFOAM so hopefully this isn't a simple google away... anyway

I am trying to model a problem to compare to an analytic solution and this requires a very specific heat source term:

V = Q*exp(-z/h)

z is the vertical coordinate, Q and h are constants, so essentially a heat source that decays with height, but only in the middle of the domain. Using heat flux at the wall won't be the same, and I've looked into using fvOptions but all I've successfully done is create a constant heat source throughout the region. This was using a semiImplicitSource - I feel like using a codedSource might be my best option for a position dependent source, but I'm having no luck in getting it to work (probably due to lack of C++ skills). Is there an easier way to achieve what I'm doing? I found a few examples of fvOptions using time dependent sources and I tried to adapt them but none of them worked.

Cheers
sturgeon

[Tobi]

Hi and welcome to the forum,

you are on the right way. Personally I would implement it via the fvOptions and using a codedFixed function for that. Right now, I have no access to my server but I was thinking that I had published a fvOptions + coded tutorial. Just check them out later (don't know what my server is doing right now).

Another option would be to implement it direct to the solver and manipulate the matrix by using the setValues() function and a cellZone. However, based on the fact that you said you have lack c++ knowledge, the fvOptions is the easiest and fastest way. You just have to keep in mind that you should divide everything by the volume (or multiply) ... ähm I guess divide.

[sturgeon]

Originally Posted by Tobi

Hi and welcome to the forum,

you are on the right way. Personally I would implement it via the fvOptions and using a codedFixed function for that. Right now, I have no access to my server but I was thinking that I had published a fvOptions + coded tutorial. Just check them out later (don't know what my server is doing right now).

Another option would be to implement it direct to the solver and manipulate the matrix by using the setValues() function and a cellZone. However, based on the fact that you said you have lack c++ knowledge, the fvOptions is the easiest and fastest way. You just have to keep in mind that you should divide everything by the volume (or multiply) ... ähm I guess divide.

Thank you for the reply, glad to know I am on the right track. I don't mind learning c++ and am reasonably comfortable with programming in general so if that's a better option I could pursue that. Your tutorial sounds useful, could you give me a link?

Cheers
sturgeon

[Tobi]

Right now I did not find a tutorial with fvOption published on my website. But in the header file there is an example that you can use. http://www.holzmann-cfd.de/index.php/en/tutorials-en

[sturgeon]

Originally Posted by Tobi

Right now I did not find a tutorial with fvOption published on my website. But in the header file there is an example that you can use. http://www.holzmann-cfd.de/index.php/en/tutorials-en

I do apologise but I have been looking over your site for a while and I can't find the example, can you be more specific? Thank you

EDIT:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



harmonic
{

    type            scalarCodedSource;

    active          true;
    selectionMode   all;

    scalarCodedSourceCoeffs
    {
        fields      (T);
        redirectType    sourceCoeffs;
		selectionMode   all;

        codeInclude
        #{

        #};

        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};

        codeAddSup
        #{
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
	    const volVectorField& U = mesh().lookupObject<volVectorField>("U");
	    const volScalarField z = U.mesh().C() & vector(0,1,0);

            scalarField& hSource = eqn.source();
            forAll(V, i)
            {
            	hSource[i] -= 0.00098*exp(-z[i]/1000)*V[i];
            }
            Pout << "***codeAddSup***" << endl;
        #};

        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};

        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }

    sourceTimeCoeffs
    {
			selectionMode   all;
        // Dummy entry
    }
}


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

This is what I have created so far, adapted from another fvOptions I found. It works in the sense that it executes without error, although honestly I am not sure what all the bits of the code do. The code I copied was also for a heat source, so I am not sure why vector fields are involved, but when I tried trimming them out and adapting the code I kept meeting errors. The temperatures produced seem realistic although my velocities start approaching relativistic levels fairly quickly, so I imagine I have issues elsewhere.

Aside from assigning the source to the entire domain rather than a selected portion of it, is anyone able to tell me if this is in line with what I am trying to achieve? Thanks

[Tobi]

As I said, I could not find the tutorial. Maybe I did not publish it because I do not publish all the tutorials. However, my second hint was the source code header:

   energySource
    {
        type            scalarCodedSource;

        active          yes;

        name    sourceTime;

        scalarCodedSourceCoeffs
        {
            selectionMode   all;

            fields          (h);

            codeInclude
            #{

            #};

            codeCorrect
            #{
                Pout<< "**codeCorrect**" << endl;
            #};

            codeAddSup
            #{
                const Time& time = mesh().time();
                const scalarField& V = mesh_.V();
                scalarField& heSource = eqn.source();
                heSource -= 0.1*sqr(time.value())*V;
            #};

            codeSetValue
            #{
                Pout<< "**codeSetValue**" << endl;
            #};

            // Dummy entry. Make dependent on above to trigger recompilation
            code
            #{
                $codeInclude
                $codeCorrect
                $codeAddSup
                $codeSetValue
            #};
        }

        sourceTimeCoeffs
        {
            $scalarCodedSourceCoeffs;
        }
    }

[sturgeon]

Originally Posted by Tobi

As I said, I could not find the tutorial. Maybe I did not publish it because I do not publish all the tutorials. However, my second hint was the source code header:

Thanks - this is one of the files I found online and was trying to adapt.

[Tobi]

Then you should specify your question more. Right now there is nothing one can get out of your questions. So what is the problem?

[sturgeon]

Originally Posted by Tobi

Then you should specify your question more. Right now there is nothing one can get out of your questions. So what is the problem?

When I originally posted I couldn't get any fvOptions to work, so I suppose there is none now. Although perhaps you could help me understand something:

In an effort to understand exactly what is going on, I am trying to figure out how to see the internal variables during execution, ie, in:

        codeAddSup
        #{
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
	    const volVectorField& U = mesh().lookupObject<volVectorField>("U");
	    const volScalarField z = U.mesh().C() & vector(0,1,0);

            scalarField& hSource = eqn.source();
            forAll(V, i)
            {
            	hSource[i] -= 0.00098*exp(-z[i]/1000)*V[i];
            }
            Pout << "***codeAddSup***" << endl;
        #};

I would love to be able to get readouts or V, C, U, z and hSource so I understand exactly what is happening. I have looked up how to print to console but the code never seems to be executed. I suppose the subtleties of C++ are escaping me. I come from MATLAB where debugging is a simple case of clicking to set a breakpoint and mousing over the variables.

Thanks for your help so far Tobi.

[Tobi]

Are you aware of doxygen? If not, use it. The class fvOptions holds a reference to the object of the fvMesh class named mesh_ which you can access:

const volVectorField& U = mesh().lookupObject<volVectorField>("U"); 
const volScalarField z = U.mesh().C() & vector(0,1,0);

I would do it as follows:

forAll(C, cellI)
{
    heSource[cellI] -= 0.00098*exp(-C[i].z()/1000)*V[i];
}

If the outputting works (don 't know - never tested it):

forAll(C, cellI)
{
    Pout<< "For cellID " << cellI 
        << " we add a source of " << 0.00098*exp(-C[cellI].z()/1000)*V[cellI] 
        << " while the z coordinate is " << C[cellI].z() << nl;

    heSource[cellI] -= 0.00098*exp(-C[cellI].z()/1000)*V[cellI];
}

[Tibo99]

Hi!

I got similar issue.

I applied only one source term on epsilon to the k-e turbulence model for 1D (x=0.001m, y=0.15m, z=0.001m) ABL simulation. Ny=150, so each cells center are at y1=0.0005m, y2=0.0015m, etc.

Y+ is around 50. The y component I used in the formula is the center of the cells so, y-coordinate of internal cells (distance from the wall to the center of each cells,see text file, 3rd column).

From the results, the calculation by the equation "epsilonSource[i]" are good since I printed out the results in a text file and verified the validity but the source term doesn't seem to be applied to the domain(see picture). The boundary wall seem to work (because in the epsilon patch I choose the BC "epsilonWallFunction") and then after, the values are null. The simulation stop after 11 step of iteration...

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1906                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

epsilonSource
{
    type            		scalarCodedSource;
    selectionMode   		all;
    active			true;

    fields          		(epsilon);
    name			epsilonSource;
	
    
    codeInclude
    	#{
		#include "OFstream.H"
	#};
    
    codeAddSup
    	#{
	const Time& time = mesh_.time();
				
        const scalar us_ = 0.76;
        const scalar D1k_ = 1.44;
        const scalar D2k_ = 1.92;
        const scalar kappa_ = 0.41;
        const scalar sigmaEps_ = 1.17;
        const scalar Cmu_ = 0.09;
        const scalar z0_ = 3.3474e-05;
        
        const scalarField& y = mesh_.C().component(1);
        
	scalarField& epsilonSource = eqn.source();
				
        // Apply the source
        forAll(y, i)
        {
          
            {
               epsilonSource[i] = (pow(us_,4)/(pow((y[i]+z0_),2)))
               *(((D2k_-D1k_)*sqrt(Cmu_)/pow(kappa_,2))-(1/sigmaEps_));
            }   
                               
			   			
			std::ofstream file;
			file.open ("Epsilonterm.txt", std::ofstream::out | std::ofstream::app);
			file << time.value()<< "	" << epsilonSource[i] << "	" << y[i] << endl << "\n";
			file.close();
						
        
        };
        
        
    #};
    
    codeCorrect
    #{
    #};
    
    codeConstrain
    #{
    #};
}
// ************************************************************************* //

If you look in the text file, 3rd column, the first 4th y-coordinate seems to be wrong and I don't know why this happen.

Do you guys have any clue?

Thank you very much!

Best Regards,

P.S.: I have tried your method and no luck.

[uosilos]

Originally Posted by Tobi

Are you aware of doxygen? If not, use it. The class fvOptions holds a reference to the object of the fvMesh class named mesh_ which you can access:

const volVectorField& U = mesh().lookupObject<volVectorField>("U"); 
const volScalarField z = U.mesh().C() & vector(0,1,0);

I would do it as follows:

forAll(C, cellI)
{
    heSource[cellI] -= 0.00098*exp(-C[i].z()/1000)*V[i];
}

If the outputting works (don 't know - never tested it):

forAll(C, cellI)
{
    Pout<< "For cellID " << cellI 
        << " we add a source of " << 0.00098*exp(-C[cellI].z()/1000)*V[cellI] 
        << " while the z coordinate is " << C[cellI].z() << nl;

    heSource[cellI] -= 0.00098*exp(-C[cellI].z()/1000)*V[cellI];
}

Hello Tobi,


how can I access a topoSet defined cellSet inside the loop??


I have tried

forAll(mesh.cellZones()[ignitionCellsCZ],i)

and it says this cellZone has not been declared in this scope


Moreover, what is the difference between mesh and mesh_??


Thanks