CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   creating a new field by volume averaging of an existing field and reuse it in Solver (https://www.cfd-online.com/Forums/openfoam-programming-development/210796-creating-new-field-volume-averaging-existing-field-reuse-solver.html)

Gbaz November 4, 2018 15:13

creating a new field by volume averaging of an existing field and reuse it in Solver
 
Hello Foamers,

Using VolfieldValue, I created a new field based on volume averaging of an existing field. Now, I don't know how to ask Solver to call and read the new field and use it for its next time step.

Your help is appreciated.

Thanks in advance!

clapointe November 4, 2018 16:24

We're going to need some more information to help you. What solver/application is this? What code do you already have? How do you want to use the new field that you've created?

Caelan

Gbaz November 5, 2018 00:58

Hi Caelan,

Thanks for your help.

I have a code based on PimpleFoam, modified for solidification. The soldification temperature (Tsol) is a function of total liquid fraction (alphAvg). Therefore, I defined a new field (using volFieldValue in controlDict) that calculates the volume average of alpha. The output file can be seen properly in my postprocessing.

Now, my problem is how to read (call) this output FROM Solver (heattransfer/Pimplefoam) to be used for the next time step.
As an option, I just simply created Tsol as another volScalarfield in CreateFields.H. Then tried to relate Tsol to alphAvg there. but it seems Tsol is NOT UPDATED by alphAvg at each time step and only uses the initial alphAvg. somehow I think Solver can't read alphAvg in postprocessing folder.

Thanks for your time.

clapointe November 5, 2018 12:26

Ok, I think I'm getting it. My first thought is this : creating a field or doing field manipulation in createFields will only happen once because createFields is only called once at the beginning to initialize fields. If you want the field to be updated, you'll need to put code inside of the runtime loop.

Caelan

Gbaz November 5, 2018 13:36

I think you are absolutely right.


what I did in controlDict for defining the new field alfAvg is:

functions
alfAvg
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
enabled true;
writeFields no;
valueOutput true;
regionType all;
writeControl timeStep;
writeInterval 1;
operation volAverage;
fields
(
alpha
);
}


and in createField.H, I tried to call this new field by:

Info<< "Reading field alfAvg\n" <<endl;
volScalarField alfAvg
(
IOobject
(
"alfAvg",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

then I tried to use alfAvg in energy equation (TEqn) for next steps of the solution...and as you mentioned its called only once at the beginning.

How can I define something similar in the run loop? basically where is run loop? Energy equation solver (TEqn)?

any change needed in pimpleFoam.c or elsewhere?

your reply is highly appreciated.

cheers

clapointe November 5, 2018 13:46

I see. So using the function object here can be thought of doing something "on top" of the base solver -- it doesn't know to communicate with the rest of the solver (ie both ways -- it is a one way communication). So what we can do is create the field like you have in createFields and then compute an average inside the loop (eg each timestep).

The run loop is inside pimpleFoam.C. You'll see a while statement. Then, after the time integration we can put something like :

Code:

scalar alfAvg = gSum(alpha*mesh.V())/mesh.V().size().
This will create a scalar alfAvg that is the average alpha in the domain (scaled by cel volume). Then, you can use alfAvg in Teqn.H. Beware that the specifics of the scalar (ie scalar vs dimensionedScalar, etc) and getting it to work with the rest of the solver can be an issue, but it's a solvable one.

Caelan

Gbaz November 6, 2018 13:14

Dear Caelan,

Thanks to your post, I could modified the code. However, still I think the alfAvg is not engaged with other fields.

What I did is this:

for volume fraction, alfAvg, I removed its functionObject of volFieldValue in controlDict and also from createField.H.

instead, as you suggested, I added this line in pimpleFoam.C:

scalar alfAvg = gSum(alpha*mesh.cells().size())/(mesh.V().size())

from the outputs, I can see that the alfAvg is calculated correctly. However, I defined a new scalar Tsol (solidification temperature) based on alfAvg as:

dimensionedScalar Tsol = 0.5*alfAvg;

later on, in Teqn file, Tsol will be used in energy equation solving as I wish. but from outputs, I can see that Tsol is not updated based on alfAvg since it remains equal to its initial value!!

as an option, I tried to define Tsol=0.5*alfAvg inside Teqn (and not in pimpleFoam), but it didn't change.

meanwhile after compiling the solver, I get this warning... unused variable'alfAvg'....

ps. I added the line of scalar alfAvg = gSum..... to runTime loop at different locations of pimpleFoam (inner loops, outer loop) but I didn't help.

I don't know how to call and use alfAvg inside Solver.
:confused:

by the way, in postProcessing, I have an output file for alfAvg. Is it proper to ask Solver to read alfAvg from their for Tsol?

Regards,
Gbaz

clapointe November 6, 2018 13:19

It is incredibly hard to follow what you've done without the code tags/ providing more code (eg the Teqn). Also providing output/ error from a log file helps. Additionally, info statements are helpful in general for debugging.

Caelan

Gbaz November 6, 2018 13:49

Thanks for your advise. I rewrite it as this. "Tsol" is in TEqn.

You are right. I forgot !

This is runTime loop:

Code:

while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "TEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
     
        }


        runTime.write();

        scalar alfAvg = gSum(alpha*mesh.cells().size())/(mesh.V().size())
                dimensionedScalar Tsol = 0.5*(Tl+Ts)*pow(alfAvg,10);



        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }


    Info<< "End\n" << endl;

    return 0;
}


and this is part of TEqn, where I need Tsol:

Code:

fvScalarMatrix TEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*fvm::ddt(T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
    );

    TEqn.relax();
    TEqn.solve();


I hope this helps otherwise please let me know what should I add.

clapointe November 6, 2018 13:59

For future reference, "code tags" are used to isolate chunks of code within a comment,

Code:

like this
The # button will insert them as you create a comment, and you place your code in between them.

As for your problem, I don't even see Tmelt (I'm assuming this is your quantity of interest) is used in Teqn.

Caelan

Gbaz November 6, 2018 14:27

with your guidance, I edited my previous thread. Hope you can see it in a better way now. Appreciate your comment.

anon_q November 6, 2018 14:31

Code:

runTime.write(); ///<<<<<< This should be after Tsol

scalar alfAvg = gSum(alpha*mesh.cells().size())/(mesh.V().size())
dimensionedScalar Tsol = 0.5*(Tl+Ts)*pow(alfAvg,10);

I mean:

Code:



scalar alfAvg = gSum(alpha*mesh.cells().size())/(mesh.V().size())
dimensionedScalar Tsol = 0.5*(Tl+Ts)*pow(alfAvg,10);

runTime.write(); ///<<<<<< This should be after Tsol


clapointe November 6, 2018 14:32

Much clearer, thanks. I see now where Tsol is going. I'd recommend putting the code to calculate Tsol inside Teqn.H. Then creating a volScalarField (or uniformDimensionedField) of Tsol (name it something else). Then you ensure computation of alfAvg/Tsol right before you need it. You might need to create a volScalarField or uniformDimensionedField of Tsol to subtract it from T.

Caelan

Gbaz November 6, 2018 14:50

Thanks Linda. I did so and now applying Caelan's comment.

Gbaz November 6, 2018 15:05

Caelan:I'd recommend putting the code to calculate Tsol inside Teqn.H.
Gbaz: I removed it from pimplefoam.c and added it to TEqn as this:


Code:


       
Code:

       
dimensionedScalar Tsol=0.5*(Tl+Ts)*pow(alfAvg,10);

fvScalarMatrix TEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*fvm::ddt(T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
    );

    TEqn.relax();
    TEqn.solve();



Caelan:Then creating a volScalarField (or uniformDimensionedField) of Tsol (name it something else).

Gbaz: In CreateFields.H, I added this:

Code:

    volScalarField Tsol
  (
    IOobject
    (
      "Tsol",
      runTime.timeName(),
      mesh,
      IOobject::MUST_READ,
      IOobject::AUTO_WRITE
    ),
       
    mesh
    );

Name it as Tmelt is correct?:confused:

Caelan:You might need to create a volScalarField or uniformDimensionedField of Tsol to subtract it from T.

Gbaz: Is this something more than what I did in createfields.H ?

Gbaz November 6, 2018 15:33

Applying above comments, outputs are still showing that Tsol is not updated. This means alfAvg is not also updated and remains equal to initial value of alfAvg.

by the way, I tried to remove initial value file of alfAvg in 0 direcroty. but running failed with this message:

FOAM FATAL ERROR:
cannot find file "C:/PROGRA~1/BLUECF~1/ofuser-of5/run/SolidiFoam/0/alfAvg"

:confused:

clapointe November 6, 2018 16:03

Sorry for the confusion -- I didn't realize you'd created a volScalarField named Tsol. So, you don't need "dimensionedScalar" before Tsol in Teqn.H. ...unless Ts and Tl volScalarFields, which it looks like they could be.

If they are volScalarFields, then your Teqn.H looks like it'll be :

Code:

Tsol=0.5*(Tl+Ts)*pow(alfAvg,10);

fvScalarMatrix TEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*fvm::ddt(T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
    );

    TEqn.relax();
    TEqn.solve();

If not, a messy/rudimentary solution would be :

Code:


forAll(Tsoll, celli)
{
    Tsol[celli]=0.5*(Tl+Ts)*pow(alfAvg,10);
}

fvScalarMatrix TEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*fvm::ddt(T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
    );

    TEqn.relax();
    TEqn.solve();

Or a cleaner version, using Tsol as a dimensionedScalar (and removing its declaration as a volScalarField from createFields.H)

Code:


dimensionedScalar Tsol("Tsol", dimTemperature, 0.5*(Ts+Tl)*pow(alfAvg,10));

fvScalarMatrix TEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*fvm::ddt(T)
      + hs*4.0*exp(-pow(4.0*(T-Tsol)/(Tl-Ts),2))/Foam::sqrt(pi)/(Tl-Ts)*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
    );

    TEqn.relax();
    TEqn.solve();

Note that I quickly tested all options and they compiled using a the dev branch. I've not tried running anything, though.

Also note that for your calculation of alfAvg you dropped the mesh.V() volume normalization part.

Replying to your most recent post, this is because createFields pulls intial conditions from your 0 directory to create the fields needed for the solver to run. It does not write files. Additionally, when we create the scalar alfAvg, we are just creating it for use within the solver -- it is a singular value that is not written. You'll likely want to add it to the Teqn.H file as well to make sure it's updated before Tsol. Otherwise you'll need to initialize it to some value (eg in craeteFields).

Caelan

Gbaz November 7, 2018 10:37

Thanks Caelan for your help.
:)

Applying all your comment, now I want to trace Tsol values over the the domain. basically it should be changing as alfAvg is changed as

Code:

Tsol=0.5*(Tl+Ts)*alfAvg;

the output file of postprocessing/alfAvg/0/volFieldValue.dat shows that alfAvg is updating properly (from initial value of 1 to final value of 0.84) but Tsol is always equal to 0.5*(Tl+Ts) means alfAvg have been taken as 1 which is the initial value! does it mean Tsol is not updated!

by the way, what's the easiest way to monitor alfAvg and Tsol on screen (or log file) as they fly?

cheers,
Gbaz

clapointe November 7, 2018 10:40

Like you said, this means that alfAvg is not being updated. Where did you place the alfAvg code? And do you declare it in createFields?

An easy way to see values is as follows :

Code:


//see alfAvg

Info << "\nalfAvg : " << alfAvg << nl << endl;

//see Tsol

Info << "\nTsol : " << Tsol << nl << endl;

Place those lines after you compute alfAvg/Tsol (eg in Teqn.H, if you put the calculations there), and the information will be printed to the terminal or the log file.

Caelan

Gbaz November 7, 2018 13:53

Caelan: Where did you place the alfAvg code? And do you declare it in createFields?

Gbaz: alfAvg is in pimpleFoam, at the end of runTime loop as:

Code:

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "TEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
     
        }



        scalar alfAvg = gSum(alpha*mesh.V())/(mesh.V().size())
                Info << "\nalfAvg:"<<alfAvg<<nl<<endl
               

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

and declared in createFields.H as:

Code:

Info<< "Reading field alfAvgM\n" <<endl;
    volScalarField alfAvgM
    (
      IOobject
      (
      "alfAvgM",
      runTime.timeName(),
      mesh,
      IOobject::MUST_READ,
      IOobject::AUTO_WRITE
      ),
      mesh
    );





Caelan: Place those lines after you compute alfAvg/Tsol (eg in Teqn.H, if you put the calculations there), and the information will be printed to the terminal or the log file.


Gbaz: I did so as you see ablve, but still can't see them on screen. This is my screen:

Code:

Time = 0.548

Courant Number mean: 0.000396649 max: 0.00334194
PIMPLE: iteration 1
DILUPBiCG:  Solving for T, Initial residual = 0.000295952, Final residual = 1.74438e-009, No Iterations 1
DICPCG:  Solving for p_rgh, Initial residual = 0.0325929, Final residual = 0.000182117, No Iterations 73
time step continuity errors : sum local = 4.28928e-006, global = 4.26282e-006, cumulative = 0.0199114
DICPCG:  Solving for p_rgh, Initial residual = 0.0296605, Final residual = 7.80347e-009, No Iterations 91
time step continuity errors : sum local = 4.26286e-006, global = 4.26282e-006, cumulative = 0.0199156
PIMPLE: iteration 2
DILUPBiCG:  Solving for T, Initial residual = 7.07245e-006, Final residual = 1.58461e-011, No Iterations 1
DICPCG:  Solving for p_rgh, Initial residual = 0.0298471, Final residual = 0.000182628, No Iterations 73
time step continuity errors : sum local = 4.28937e-006, global = 4.26295e-006, cumulative = 0.0199199
DICPCG:  Solving for p_rgh, Initial residual = 0.0296423, Final residual = 7.23257e-009, No Iterations 91
time step continuity errors : sum local = 4.26298e-006, global = 4.26295e-006, cumulative = 0.0199242
PIMPLE: iteration 3
DILUPBiCG:  Solving for T, Initial residual = 2.07291e-006, Final residual = 1.33683e-011, No Iterations 1
DICPCG:  Solving for p_rgh, Initial residual = 0.029479, Final residual = 0.000153555, No Iterations 73
time step continuity errors : sum local = 4.28505e-006, global = 4.26284e-006, cumulative = 0.0199284
DICPCG:  Solving for p_rgh, Initial residual = 0.0296139, Final residual = 8.03922e-009, No Iterations 91
time step continuity errors : sum local = 4.26287e-006, global = 4.26284e-006, cumulative = 0.0199327
ExecutionTime = 5.803 s  ClockTime = 6 s



All times are GMT -4. The time now is 16:27.