[nandiganavishal]

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]

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]

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]

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

Regards

Vishal

[nandiganavishal]

[nandiganavishal]

poissonFoam.tar.gz

[nandiganavishal]

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]

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]

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]

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]

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]

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]

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]

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]

I cannot find the error

[maalan]

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]

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]

Originally Posted by smehdi609

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]

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.