CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Solve poisson equation just add a source term

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

Like Tree2Likes
  • 1 Post By nandiganavishal
  • 1 Post By nandiganavishal

Reply
 
LinkBack Thread Tools Display Modes
Old   November 3, 2008, 19:32
Default Hi I am new to OpenFoam, I w
  #1
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   November 6, 2008, 15:19
Default Hi, I do these stuff in this
  #2
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
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 is offline   Reply With Quote

Old   November 7, 2008, 17:57
Default Hi, do the following changes
  #3
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
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
smehdi609 is offline   Reply With Quote

Old   November 10, 2008, 18:06
Default I am also attaching the entire
  #4
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
I am also attaching the entire code here. Kindly let me know where I am making an error.

Regards

Vishal
nandiganavishal is offline   Reply With Quote

Old   November 10, 2008, 18:09
Default
  #5
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road

nandiganavishal is offline   Reply With Quote

Old   November 10, 2008, 18:13
Default http://www.cfd-online.com/Ope
  #6
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
poissonFoam.tar.gz
mm.abdollahzadeh likes this.
nandiganavishal is offline   Reply With Quote

Old   November 10, 2008, 19:06
Default A small change my solver looks
  #7
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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 is offline   Reply With Quote

Old   November 10, 2008, 20:19
Default Hi Mahdi, Thanks a lot for
  #8
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   November 10, 2008, 20:47
Default Hi, 2. Yeah you are right it
  #9
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
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
smehdi609 is offline   Reply With Quote

Old   November 10, 2008, 20:59
Default Hi Mahdi, Thanks a ton for
  #10
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   November 10, 2008, 21:56
Default Hi, For example you can do t
  #11
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
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
smehdi609 is offline   Reply With Quote

Old   November 10, 2008, 22:27
Default Hi Mahdi, So Inlet function
  #12
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   November 11, 2008, 18:06
Default Hi, You should add them to y
  #13
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
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
smehdi609 is offline   Reply With Quote

Old   November 11, 2008, 18:30
Default Hi Mahdi, I incorporated th
  #14
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9
nandiganavishal is on a distinguished road
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
mm.abdollahzadeh likes this.
nandiganavishal is offline   Reply With Quote

Old   November 11, 2008, 19:04
Default I cannot find the error
  #15
Member
 
M. Mahdi Salehi
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 50
Rep Power: 8
smehdi609 is on a distinguished road
I cannot find the error
smehdi609 is offline   Reply With Quote

Old   July 14, 2013, 14:21
Default
  #16
Member
 
Join Date: Jun 2011
Posts: 65
Rep Power: 6
maalan is on a distinguished road
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!
maalan is offline   Reply With Quote

Old   August 28, 2014, 13:02
Default nandiganavishal and Mahdi
  #17
New Member
 
Shakil Masum
Join Date: Aug 2014
Posts: 10
Rep Power: 3
Shakil Masum is on a distinguished road
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
Shakil Masum 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 solve for transport with large source term yong Main CFD Forum 2 February 2, 2007 10:31
How to add a source term in u or v equation? luckyluke Phoenics 3 December 13, 2004 03:52
Can we solve only convective and source term Mehdi Main CFD Forum 1 February 18, 2004 00:40
How to solve Poisson equation in CFX 4? windhair CFX 3 January 27, 2004 22:09
Solve Poisson equation Linfeng BI Main CFD Forum 4 November 21, 2002 00:46


All times are GMT -4. The time now is 01:17.