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

how to properly add a forcing term to laplaces equation?

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

Reply
 
LinkBack Thread Tools Display Modes
Old   May 18, 2018, 05:47
Default how to properly add a forcing term to laplaces equation?
  #1
KTG
New Member
 
Abe
Join Date: May 2016
Posts: 29
Rep Power: 4
KTG is on a distinguished road
Hi all,


I am trying to solve poisson's equation by adding a forcing term to laplacianFoam and removing the time derivative. My code runs, but the solution around the forcing is way off - I get large numbers that depend on the size of the mesh, but the solutions at the Dirichlet boundaries make sense. This lead me to believe that when I write:


Code:
           

 fvScalarMatrix psixEqn
            (
                fvm::laplacian(DT, psix) + f1
            );
openfoam does not multiply by the cell volume but rather inserts the forcing term f1 directly into the matrix problem? After a lot of internet searching, I modified the forcing this way in my createFields.H:




Code:
Info<< "Creating forcing term F1"<<endl;
volScalarField F1 (
        IOobject (
            "F1",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ
       ),
       mesh,
       dimensionedScalar("F1",dimensionSet(0,-3,0,0,0,0,0), 0.0),
       calculatedFvPatchScalarField::typeName
);

F1.ref() = f1/mesh.V();
which gave results closer to what I expected, though it is still off depending on the mesh size. I also feel like I should be multiplying by the cell volume rather than dividing, but that gave really wild results. Clearly there is something I am not getting here - there is probably an easier way to do this that will make it work properly...help.



Thanks
KTG is offline   Reply With Quote

Old   May 21, 2018, 00:52
Default
  #2
KTG
New Member
 
Abe
Join Date: May 2016
Posts: 29
Rep Power: 4
KTG is on a distinguished road
This is such a basic thing, its very strange that I can't find a casefile anywhere that shows how to do it! Its seems like I should be using fvc::sp(f1) or fvc::su(f1), but when I use these, I get an error:


Code:
error: no matching function for call to ‘Sp(Foam::volScalarField&)’
                 fvm::laplacian(DT, psix) == fvc::Sp(f1)

where the ^ points to the last ) after f1


I have defined f1 this way:


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

Maybe volScalarField is the problem? I don't know what else to put there, and am not finding any clues in fvcSup.H (if that is even the right file).
KTG is offline   Reply With Quote

Old   May 21, 2018, 11:38
Default
  #3
Member
 
Join Date: Aug 2015
Posts: 74
Rep Power: 4
clapointe is on a distinguished road
What version of OpenFOAM are you using? I quick search of the code (version 4) returns many uses of fvm::Sp, but none of fvc::Sp, which I think is the reason you are getting the error.

Caelan
clapointe is offline   Reply With Quote

Old   May 22, 2018, 02:54
Default
  #4
KTG
New Member
 
Abe
Join Date: May 2016
Posts: 29
Rep Power: 4
KTG is on a distinguished road
Yes, I think that you are correct for the reason for the error - maybe it is outdated, but I read that fvc:: should explicitly add the source term to the right hand side. My original question stands though. I used this code:


Code:
            fvScalarMatrix psixEqn
            (
                fvm::laplacian(DT, psix) == fvm::Su(f1,psix)
            );

and am getting the same results as before - the forcing terms seems to vary with the size of the mesh, and is negative when I put it on the right hand side. When I set f1=1, it comes out to about 1800 in the solution for psix. Do you know if openfoam is just setting the source term equal to f1, or actually multiplying it by the cell volume?
KTG is offline   Reply With Quote

Old   May 22, 2018, 09:33
Default
  #5
Member
 
Join Date: Aug 2015
Posts: 74
Rep Power: 4
clapointe is on a distinguished road
It looks like the strength is scaled by cell volume -- take a look at (for Openfoam 4x): https://github.com/OpenFOAM/OpenFOAM...e/fvm/fvmSup.C.

Caelan
clapointe is offline   Reply With Quote

Old   May 22, 2018, 10:29
Default
  #6
KTG
New Member
 
Abe
Join Date: May 2016
Posts: 29
Rep Power: 4
KTG is on a distinguished road
Yes, I think you are right, which makes sense in terms of integrating over the cell in the finite volume method. What confuses me is that the result is the same when I just insert the forcing term directly into the equation without fvm::Sp(f1):


Code:
            fvScalarMatrix psixEqn
            (
                fvm::laplacian(DT, psix) == f1
            );


Anyway it seems like using fvOptions is the preferred way to do things now, so maybe I should do that instead. I think some of my initial confusion was due to poor mesh quality anyway - thanks for the help!
KTG is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
How to add implicit terms to the momentum equation? MaryBau OpenFOAM Programming & Development 11 December 26, 2014 22:46
Calculation of the Governing Equations Mihail CFX 7 September 7, 2014 06:27
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 15:45
What is the Boussinesq Term in Momentum Equation CFDtoy Main CFD Forum 0 August 11, 2008 09:56
bouyancy term in epsilon equation Michael Main CFD Forum 1 June 25, 1999 10:20


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