CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Implicit source term

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree10Likes
  • 1 Post By Lieven
  • 8 Post By henrik
  • 1 Post By mcdonalds

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 29, 2013, 15:16
Default Implicit source term
  #1
Member
 
Ronald McDonald
Join Date: Jul 2012
Posts: 38
Rep Power: 13
mcdonalds is on a distinguished road
Dear openFoamers,


I'd like to add an implicit source term to my pde. It would look like this:

HTML Code:
        solve
            (
                Foam::fvm::laplacian(sigmaeff, phi) == -iagCell*((Foam::exp(alphaaa*F*(phiao-phi-EaEqgCell)/(Rg*T)))-(Foam::exp(-alphaca*F*(phiao-phi-EaEqgCell)/(Rg*T)))) 
            );
Where phi is my volscalarfield (Notice phi is on both sides of the equation). Other variables are independent or constant.

I looked at OpenFoam documentation, specifically from this site:
http://www.foamcfd.org/Nabla/guides/...sGuidese9.html

And it talks about using Su, Sp, and SuSp. I would like to use those things but I'm not sure how to. In the documentation it talks about a rho and phi, where rho can either be a dimensioned scalar or a volscalarField, and phi is the volscalarfield.

Now in my pde I do not know what rho would be.

Also, I found another link talking about this:
http://openfoamwiki.net/index.php/Ho...sport_equation

Any help you, the community, could provide would be great!

Sincerely,

Benjamin
mcdonalds is offline   Reply With Quote

Old   April 29, 2013, 16:52
Default
  #2
Senior Member
 
Lieven
Join Date: Dec 2011
Location: Leuven, Belgium
Posts: 299
Rep Power: 22
Lieven will become famous soon enough
Hi Benjamin,

If you do a incompressible simulation, rho is simply 1.0.
Now, regarding your source term. The fact that it is inside the exponential function makes it impossible to included implicitly. This is only possible for linear systems. The closest you will probably get is something like
Code:
fvm::Sp(exp(...)/phi,phi)
but I don't know how much effect this will have... It's certainly not fully explicit.

Cheers,

L
granzer likes this.
Lieven is offline   Reply With Quote

Old   April 30, 2013, 14:01
Default
  #3
Member
 
Ronald McDonald
Join Date: Jul 2012
Posts: 38
Rep Power: 13
mcdonalds is on a distinguished road
Dear openFoamers,

Thank you for your response, Lieven. I suppose I cannot implement an implicit source term in that manner, specifically using the Foam::fvm::Sp function of openFoam.

Then I am stumped on how to implement the pde in openFoam. My equation with more clarity is attached to this file using Latex.




where phi_e is the volscalarField and all others are either constant or independent of phi.

When I put this equation into openFoam directly I get some funky results, incorrect as compared to similar simulations. Thus, I think it has to do with the phi on the right side of the equation and the equation as an implicit pde.

Any help would be really great as I feel that this is my last and final hurdle to get my simulation running smoothly.

Sincerely,

Benjamin
Attached Images
File Type: jpg Selection_043.jpg (7.6 KB, 204 views)
mcdonalds is offline   Reply With Quote

Old   May 1, 2013, 07:23
Default
  #4
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18
henrik is on a distinguished road
Dear Benjamin,

Yes, I am not surprised that you are having trouble with that source ... Did you try to underrelax?

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == -iagCell*((Foam::exp(alphaaa*F*(phiao-phi-EaEqgCell)/(Rg*T)))-(Foam::exp(-alphaca*F*(phiao-phi-EaEqgCell)/(Rg*T))))
);
phiEqn.relax(0.5);
phiEqn.solve();

Make sure to b build/solve multiple times!

I order to make it (more) implicit, you need to linearise your source term. There are many ways, but I suggest to use Picard's method.

http://www.cfd-online.com/Wiki/Sourc..._linearization

Once you get through the maths, the code will look


volScalarField Sp = // derived equations
volScalarField Sc = // derived equations

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == fvm::Sp(Sp, phi) + Sc
);
phiEqn.relax(0.5);
phiEqn.solve();

Best Regards,

Henrik

Last edited by henrik; May 1, 2013 at 11:14.
henrik is offline   Reply With Quote

Old   May 8, 2013, 13:16
Default
  #5
Member
 
Ronald McDonald
Join Date: Jul 2012
Posts: 38
Rep Power: 13
mcdonalds is on a distinguished road
Dear openFoamers,

Thank you, Henrik, for your help. I am still trying to linearize but I have another idea.

I would like to do the following:

Etta = phiao - phi - EaEqgCell;

If Etta > .5 then Etta = .5;

Else Etta = phiao - phi -EaEqgCell;

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == -iagCell*((Foam::exp(alphaaa*F*(Etta)/(Rg*T)))-(Foam::exp(-alphaca*F*(Etta)/(Rg*T))))
);
phiEqn.relax(0.5);
phiEqn.solve();

It's the if statement I am having trouble with. How would I do that in my solver? How to create an if statement?

Sincerely,

Benjamin
mcdonalds is offline   Reply With Quote

Old   May 8, 2013, 17:20
Default
  #6
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18
henrik is on a distinguished road
Benjamin,

Did you try to underrelax? Did it make a difference?

Having a switch is nasty - especially if you try to linearise. Anyway, the way to do it is to use the pos and neg functions.

Alternatively, you can loop through the field cell by cell and use raw C++.

Henrik
henrik is offline   Reply With Quote

Old   May 10, 2013, 15:28
Default Pleasantly perplexed
  #7
Member
 
Ronald McDonald
Join Date: Jul 2012
Posts: 38
Rep Power: 13
mcdonalds is on a distinguished road
Dear OpenFoamers,

Hi Henrik. I did try to underrelax and it didn't make a difference, I got an error right away. I did put an if-statement in and it ran without an error. However, it leads to another problem I have.

Currently, my solver looks like this:

HTML Code:
    Info<< "\nCalculating concentration distribution\n" << endl;

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



//**Add laplacians to all four species: O2, N2, H2O, H2
        while (simple.correctNonOrthogonal())
        {
       #include "diffusion.H"
    #include "partialpressure.H"
    #include "equilibriumpotential.H"
    #include "currentdensity.H"
    #include "mapToCell.H"

    #include "Etta.H"
    #include "chargedensity.H"

    #include "electricpotential.H" 

        }




    //#include "mapFromCell.H"



    
        #include "write.H"

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

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

    return 0;
The electricpotential.H file has that crazy laplaican(sigmaeff,phi) = source pde. Now the really really interesting thing is that when I look at phi directly in paraview or even after that equation, let's say between electricpotential.H and the curly bracket }, I get an everchanging phi and it seems correct. However, when I look at phi BEFORE the electricpotential.H equation, I get a phi that changes after one time step to a linear curve between its two boundary conditions and that's it. It remains steady thereafter.

Now I think there is something going on with the "solve" command that allows openFoam to be able to get the old value and put it back into the equation. It seems like outside and BEFORE the solve equation it is not reading the old values back. But AFTER the solve equation it seems to work well.

Now to my question: how do I read back the output of the solve equation, place it before my electricpotential.H file (so that I can edit it) and then place it back into my pde? In essence, I would like to see phi changing when I check on it BEFORE my solve equation of the pde. What do you, the community of openFoamers, think about this obfuscating situation?

Sincerely,

Benjamin


Quote:
Originally Posted by henrik View Post
Benjamin,

Did you try to underrelax? Did it make a difference?

Having a switch is nasty - especially if you try to linearise. Anyway, the way to do it is to use the pos and neg functions.

Alternatively, you can loop through the field cell by cell and use raw C++.

Henrik
mcdonalds is offline   Reply With Quote

Old   May 11, 2013, 13:36
Default
  #8
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18
henrik is on a distinguished road
Benjamin,

underrelaxation MUST change the residuals - If they do not change, you are not under-relaxing!

I am entirely sure that I understand what you are trying to do, but I will provide enough rope to you in order to hang yourself

The solve statement writes the result into field it is operating on. If you put loops into your code, the explicit terms see the new (time) values. The idea here is to iterate until the (non-linear) system of equations converges.

Of course, ddt() needs the old time values to function correctly and that's why they are stored automagically. You can retrieve them by calling .oldTime() on the field.

Henrik
henrik is offline   Reply With Quote

Old   May 29, 2013, 17:47
Default
  #9
Member
 
Ronald McDonald
Join Date: Jul 2012
Posts: 38
Rep Power: 13
mcdonalds is on a distinguished road
Thanks, Henrik, for your help. I changed the fvsolutions folder to:

phi
{
solver smoothSolver;
smoother GaussSeidel;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
}
And it runs without error. I'm not sure why. And it also takes about 6x longer to run, which sucks.

I am onto a new problem, also with source terms. Currently, I have as the source terms:

HTML Code:
              sah2ogCell = jbvagCell/(F*4)*H2OMM;

              sah2gCell = jbvagCell*H2MM/(4*F);
Where sah2ogCell and sah2gCell both are volscalarfields. I then map those to my specified regions with:

HTML Code:
    {
       forAll (sah2o, cellI)
        {
            sah2o[cellI] = sah2ogCell[anodeCellMap[cellI]];
        }

    }

    {
       forAll (sah2, cellI)
        {
            sah2[cellI] = sah2gCell[anodeCellMap[cellI]];
        }

    }
Then I place them in my diffusion equations like so:

HTML Code:
        fvScalarMatrix H2OEqn
         
            (
                fvm::ddt(H2O) == fvm::laplacian(DH2Oeff, H2O) + sah2o
            );


       H2OEqn.solve();
        fvScalarMatrix H2Eqn
         
            (
                fvm::ddt(H2) == fvm::laplacian(DH2eff, H2) - sah2
            );


       H2Eqn.solve();
It compiles fine but I get a floating point exception error after the first iteration. After this if I comment out H2 and make sah2o have an internal field of 1e-14 (previously zero) then that runs without error, but the same does not hold for sah2. Any ideas why I'm getting an error when i try to run it? Is there another why to add these source terms?

Sincerely,

Benjamin


Quote:
Originally Posted by henrik View Post
Benjamin,

underrelaxation MUST change the residuals - If they do not change, you are not under-relaxing!

I am entirely sure that I understand what you are trying to do, but I will provide enough rope to you in order to hang yourself

The solve statement writes the result into field it is operating on. If you put loops into your code, the explicit terms see the new (time) values. The idea here is to iterate until the (non-linear) system of equations converges.

Of course, ddt() needs the old time values to function correctly and that's why they are stored automagically. You can retrieve them by calling .oldTime() on the field.

Henrik
babakflame likes this.
mcdonalds is offline   Reply With Quote

Old   May 30, 2013, 05:54
Default
  #10
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18
henrik is on a distinguished road
Dear Benjamin,

How many iterations does the linear solver take?

Probably a lot. If so, switch to GAMG - that should be much more efficient as laplacians are a prime target for these solvers.

Do you employ outer iterations, ie. do you solve phi multiple times per time step?

If so, Picard's method should speed up convergence.

It's hard to tell what is going wrong in your new problem. First thing to do is to write out the source fields and check that they are correct.

A likely problem is that of boundedness, ie. negative source terms will drive your species negative. I your current implementation there is nothing to prevent this. You should make the source implicit fvm::Sp and possibly clip after the solution. Have a look into reactingFoam how this can be done.

Henrik
henrik is offline   Reply With Quote

Old   September 8, 2013, 07:47
Default Implicit source term
  #11
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AUS
Posts: 137
Rep Power: 14
ahmmedshakil is on a distinguished road
Hi Henrik,
I have a stupid question to ask. What do you mean by "build/solve it multiple times!!"? Is it be done in the solve() function considering the tolerance in the fvSolution automatically?? or I have to add another loop inside the time loop like:
while(runTime.loop())
{
do
(

volScalarField Sp = // derived equations
volScalarField Sc = // derived equations

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == fvm::Sp(Sp, phi) + Sc
);
phiEqn.relax(0.5);
phiEqn.solve();
)while(...condition)
}



Quote:
Originally Posted by henrik View Post
Dear Benjamin,

Yes, I am not surprised that you are having trouble with that source ... Did you try to underrelax?

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == -iagCell*((Foam::exp(alphaaa*F*(phiao-phi-EaEqgCell)/(Rg*T)))-(Foam::exp(-alphaca*F*(phiao-phi-EaEqgCell)/(Rg*T))))
);
phiEqn.relax(0.5);
phiEqn.solve();

Make sure to b build/solve multiple times!

I order to make it (more) implicit, you need to linearise your source term. There are many ways, but I suggest to use Picard's method.

http://www.cfd-online.com/Wiki/Sourc..._linearization

Once you get through the maths, the code will look


volScalarField Sp = // derived equations
volScalarField Sc = // derived equations

fvMatrix phiEqn
(
fvm::laplacian(sigmaeff, phi) == fvm::Sp(Sp, phi) + Sc
);
phiEqn.relax(0.5);
phiEqn.solve();

Best Regards,

Henrik
ahmmedshakil is offline   Reply With Quote

Old   November 14, 2019, 09:33
Default
  #12
New Member
 
Xun Lan
Join Date: Oct 2019
Posts: 21
Rep Power: 6
Lann is on a distinguished road
Hi Benjamin,

Did you succeed to solve your equation with the implicit exponential source term? I met the same equation and stumped here.

Thanks
Lann is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem of SOURCE term gradient in UDS wind Fluent UDF and Scheme Programming 6 December 1, 2022 15:21
transient source term strakakl OpenFOAM 38 November 19, 2013 02:18
[swak4Foam] funkySetFields compilation error tayo OpenFOAM Community Contributions 39 December 3, 2012 06:18
Adding implicit source term to momentum equation fs82 OpenFOAM 6 September 23, 2009 04:29
[Gmsh] Compiling gmshFoam with OpenFOAM-1.5 BlGene OpenFOAM Meshing & Mesh Conversion 10 August 6, 2009 05:26


All times are GMT -4. The time now is 12:29.