CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Solve poisson equation just add a source term (https://www.cfd-online.com/Forums/openfoam-solving/58351-solve-poisson-equation-just-add-source-term.html)

nandiganavishal November 3, 2008 18:32

Hi I am new to OpenFoam, I w
 
Hi
I am new to OpenFoam, I want to solve a poisson equation of the form laplacian (T) = exp(x)+sin(y) (some source term).I am using the laplacianFoam for this, but I am not able to incorporate the source term in the equation. I have read the programmers guide and it stated that i need to add a volVectorField to carry out this operation. I am not able to comprehend this and could you tell me where exactly I have to add this vectorfield and how to go about coding for the same.

Thanks

Vishal

smehdi609 November 6, 2008 14:19

Hi, I do these stuff in this
 
Hi,
I do these stuff in this way; first get reference to your x and y position and define S as the source term:
volScalarField yPos = mesh.C().component(vector::Y);
volScalarField xPos = mesh.C().component(vector::x);
volScalarField S
(
IOobject
(
"S",
runTime.timeName(),
mesh,
IOobject:: NO_READ,
IOobject:: NO_WRITE,
),
mesh
);

Then do whatever you want with your source term for example:

forAll(S, counter)
{
S[counter] = exp(xPos[counter]) + sin(yPos[counter]);
}
Now solve it:
solve(fvm::laplacian(T) == S);
This is a little bit messy, someone may do it in a better way.
Best,
Mahdi

smehdi609 November 7, 2008 16:57

Hi, do the following changes
 
Hi,
do the following changes:

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

forAll(S,counter)
{
S[counter] = Foam::sin(M_PI*xPos[counter]) + Foam::sin(M_PI*yPos[counter]));
}
The first change is to be able to read and write your source term and see if you evaluate it right or wrong. Also to be able to set the dimension. You may also be able to do it in a better way. The second change was to have a source term which is
compatible with my boundary conditions. I use zero boudary condition in a 1 by 1 square. Also, I added the "Foam::" because otherwise I cannot compile it. Maybe that's not the case in your version. I use 1.4.1. Also M_PI is pi number.
have fun,
Mahdi

nandiganavishal November 10, 2008 17:06

I am also attaching the entire
 
I am also attaching the entire code here. Kindly let me know where I am making an error.

Regards

Vishal

nandiganavishal November 10, 2008 17:09


 


nandiganavishal November 10, 2008 17:13

http://www.cfd-online.com/Ope
 
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif poissonFoam.tar.gz

nandiganavishal November 10, 2008 18:06

A small change my solver looks
 
A small change my solver looks like this.. in the attachment it s wrong,,


for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
solve (fvm::poissonian(DT, T)== S);
}

nandiganavishal November 10, 2008 19:19

Hi Mahdi, Thanks a lot for
 
Hi Mahdi,

Thanks a lot for the statements.

1. You were right I didn't compile them, and hence was getting the same results.

2. Regarding the second comment, shouldn't the object be S instead of T

FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object T;(Shouldn't it be S here)
}

3. Regarding the source term factor, then should i have to multiply 1000 to S while solving the equation, that is

solve
(
fvm::laplacian(DT, T)
==
1000 *S (Instead of S)

);

Thanks Again for all the help.
Regards

Vishal

smehdi609 November 10, 2008 19:47

Hi, 2. Yeah you are right it
 
Hi,
2. Yeah you are right it should be S, but it does not make any difference.
3. You can do it here or in the creatFile.H
Enjoy
Mahdi

nandiganavishal November 10, 2008 19:59

Hi Mahdi, Thanks a ton for
 
Hi Mahdi,

Thanks a ton for all the suggestions. The program now is working perfectly fine.

Thanks a lot.

Now, I would like to incorporate boundary conditions that are not fixed. Rather they are also functions of x and y.

Like my boundary conditions are

inlet (AT x =0) - Dirichlet B.C's = -y^3 + sin(y) (something of this type).
Similarly (at y=0) (lowerwall) = -x^3 + sin (x)..

Neumann B.C's at x =1 (outlet) = -3-2*sin(x)(some arbitary function)
Neumann B.C's at y =1 (upperwall) = 4+sin(y)


Now, How exactly should I call these in the temperature file (in the 0 folder).

I searched the forum (couldn't comprehend what exactly they were saying.) . Can you give your suggestions for these.

Regards

Vishl

smehdi609 November 10, 2008 20:56

Hi, For example you can do t
 
Hi,
For example you can do this:

label inletPatchID = mesh.boundaryMesh().findPatchID("inlet");
fvPatchScalarField& inletT = T.boundaryField()[inletPatchID];
const vectorField& faceCentresInlet = mesh.Cf().boundaryField()[inletPatchID];
forAll(inletT, pointI)
{
inletT[pointI] = faceCentresInlet[pointI].y()*(1.0 - faceCentresInlet[pointI].y());
}
best,
Mahdi

nandiganavishal November 10, 2008 21:27

Hi Mahdi, So Inlet function
 
Hi Mahdi,

So Inlet function now is Dirichlet B.C given as (y*(1-y)). This is what you have done right.

For Gradient (Neumann) B.C how exactly should we go about it, as we have to consider the flux form right for this case. Could you elaborate a little bit regard to this.

Regards

Vishal

smehdi609 November 11, 2008 17:06

Hi, You should add them to y
 
Hi,
You should add them to your creatField.H file. It does not matter what you set in 0 directory, just set a fixed value of zero for example. Because, then you will modify it by those lines in creatField.H file. About Neumann bounday condition, I never used it. I searched the forum, there is one called wallHeatTransfer, but I have no idea how to set it.
Cheers,
Mahdi

nandiganavishal November 11, 2008 17:30

Hi Mahdi, I incorporated th
 
Hi Mahdi,

I incorporated them and tried it again. But this time I got the non uniform values for inlet alone and repeated the same procedure for outlet, lowerwall and upperwall, but didn't see any change. I know I would have made a simple error. Could you point out where exactly I am doing a mistake. I however incorporated the values as a list and tried it out and it is working fine. (I know this is a crude way of doing it).

Here is my poissonfoam.c file.

#include "fvCFD.H"


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

int main(int argc, char *argv[])
{

# include "setRootCase.H"

# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"

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

label inletPatchID = mesh.boundaryMesh().findPatchID("inlet");
fvPatchScalarField& inletT = T.boundaryField()[inletPatchID];
const vectorField& faceCentresInlet = mesh.Cf().boundaryField()[inletPatchID];
forAll(inletT, pointI)
{
inletT[pointI] = faceCentresInlet[pointI].y()*(1.0 - faceCentresInlet[pointI].y());
}

label outletPatchID = mesh.boundaryMesh().findPatchID("outlet");
fvPatchScalarField& outletT = T.boundaryField()[outletPatchID];
const vectorField& faceCentresOutlet = mesh.Cf().boundaryField()[outletPatchID];
forAll(outletT, pointI)
{
outletT[pointI] = 5.0 - faceCentresOutlet[pointI].y()*(2.0 - faceCentresOutlet[pointI].y());
}

label upperwallPatchID = mesh.boundaryMesh().findPatchID("upperwall");
fvPatchScalarField& upperwallT = T.boundaryField()[upperwallPatchID];
const vectorField& faceCentresUpperwall = mesh.Cf().boundaryField()[upperwallPatchID];
forAll(upperwallT, pointIII)
{
upperwallT[pointIII] = faceCentresUpperwall[pointIII].y()*(4.0 - faceCentresUpperwall[pointIII].y());
}

label lowerwallPatchID = mesh.boundaryMesh().findPatchID("lowerwall");
fvPatchScalarField& lowerwallT = T.boundaryField()[lowerwallPatchID];
const vectorField& faceCentresLowerwall = mesh.Cf().boundaryField()[lowerwallPatchID];
forAll(lowerwallT, pointIV)
{
lowerwallT[pointIV] = (2.0 -faceCentresLowerwall[pointIV].y())*(1.0 - faceCentresLowerwall[pointIV].y());
}

T.write();



Info<< "\nCalculating temperature distribution\n" << endl;

for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readSIMPLEControls.H"

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
solve (fvm::laplacian(DT, T)== 1000 * S);
}

# include "write.H"

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

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

return(0);
}


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

My T file is giving values for inlet condition but not for the outlet, lowerwall, etc... Its just taking it as zero gradient.

Thanks

Vishal

smehdi609 November 11, 2008 18:04

I cannot find the error
 
I cannot find the error

maalan July 14, 2013 14:21

Hi there!!

I am trying to solve the laplace's equation for a scalar phi over an arbitrary domain!! do you know how I should modify the laplacianFoam application?? I have noticed in the .C file the laplacian application has 2 inputs and I just need one of them: laplacian(phi). By the moment, the steps I have followed have been:

1) create my phi volScalarField in the same manner as T in the original laplacianFoam, and comment the rest

2) leave the original libraries in the top of the laplacianFoam.C file and modify the diffusion equation by

solve
(
fvm::laplacian(phi);
);

It should be very easy, but I don't see the mistake!
Thanks a lot!

Shakil Masum August 28, 2014 13:02

nandiganavishal and Mahdi
 
Hi nandiganavishal,

what did you do to the initial condition for the source "S"? did you add another file in "0" folder named "S"? If yes, is it possible to get that file? Because, i have a constant Source. the problem is I have to create a file in 0 folder named S with internal field and boundary field. since it is a constant value, i do not need to have any information for internal and boundary field but I am not sure about it. Any suggestion would be appreciated.

thanks
Shakil

alimea January 21, 2017 04:24

Quote:

Originally Posted by smehdi609 (Post 185086)
Hi,
do the following changes:

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

forAll(S,counter)
{
S[counter] = Foam::sin(M_PI*xPos[counter]) + Foam::sin(M_PI*yPos[counter]));
}
The first change is to be able to read and write your source term and see if you evaluate it right or wrong. Also to be able to set the dimension. You may also be able to do it in a better way. The second change was to have a source term which is
compatible with my boundary conditions. I use zero boudary condition in a 1 by 1 square. Also, I added the "Foam::" because otherwise I cannot compile it. Maybe that's not the case in your version. I use 1.4.1. Also M_PI is pi number.
have fun,
Mahdi

Hi
could you plz give me procedure to solve this equation:
laplacian(y)=y0
in a rectangular domain that y0 is a velocity distribution { u(x,y) } that is obtained by solving the other problem in this domain and I should import it to this solution.
I saw the probelem of

laplacian(y)=0

in tutorial but I have to change it to laplacian(y)=y0.
How can I do that?
How should I import a u(x,y) as a y0?

vinguva November 14, 2022 09:12

Seems like openfoam can't solve a general Poission Equation with the fvm::laplacian operator. Things are better with the fvc operator, but it is explicit.



fvm::laplacian(field) == S gives total garbage.


All times are GMT -4. The time now is 07:25.