Using an iterative method with openfoam solver!

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

 March 12, 2014, 02:46 Using an iterative method with openfoam solver! #1 Member   Join Date: Apr 2013 Posts: 32 Rep Power: 5 Hi every one! Can we define a volume field and then use it just for an iterative method on the wall cells? Please point a way to do this! I tried defining the volume field first. But when i am updating that using nsite[iCell]=....for only wall cells it is giving some error!

March 12, 2014, 05:12
#2
Senior Member

Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,423
Rep Power: 25
Quote:
 Originally Posted by rajcfd Hi every one! Can we define a volume field and then use it just for an iterative method on the wall cells? Please point a way to do this! I tried defining the volume field first. But when i am updating that using nsite[iCell]=....for only wall cells it is giving some error!
Can you be more specific? How do you iterate over near-wall cells, what error do you get?

March 12, 2014, 07:09
#3
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by alexeym Can you be more specific? How do you iterate over near-wall cells, what error do you get?

Hi alexeym...

In the code my aim is to find the wall temperature by iterative method when wall flux is given..Here is the code...

Info<<"Calculating the wall heat temperature"<<endl;

scalar count = 0.0;

dimensionedScalar Tnew1

(
"Tnew1",
dimensionSet (0,0,0,1,0,0,0),
scalar (380.0)
);

volScalarField Tnewn = (Tb/Tb)*Tnew1;

dimensionedScalar Tnew2

(
"Tnew2",
dimensionSet (0,0,0,1,0,0,0),
scalar (300.0)
);

volScalarField Tnewp = (Tb/Tb)*Tnew2;

dimensionedScalar Twnew1

(
"Twnew1",
dimensionSet (0,0,0,1,0,0,0),
scalar (350.0)
);

volScalarField Twnew = (Tb/Tb)*Twnew1;

volScalarField Tw =(Tb/Tb)*Twall;

label1: ///loop

count =0;

volScalarField nsite = nref*pow(((Tw-Tsat)/delTref),1.805); //site density correlation

volScalarField A2F = min((3.141*pow(Dref,2.0)*nsite), scalar(1.0)); //quenching heat flux area

volScalarField A1F = max((scalar(1.0)-A2F), scalar(0.0001)); //convective heat flux area

volScalarField mev = ((3.14/6.0)*rhoa*pow(Ddep,3.0)*fdep*nsite);

/*......Quenching, convective, evaporative heat flux.....*/

volScalarField qe = mev*latHeat; //evaporative heat flux

volScalarField qc = A1F*hc*(Tw-Tb); // convective heat flux

volScalarField qq = A2F*hq*(Tw-Tb); //quenching heat flux

volScalarField qtotal = qe+qc+qq; //total heat flux

Info<<"\nCalculated the individual contributions"<<endl;

volScalarField error = mag(qwall-qtotal);

Info<<"\nEntering the loop"<<endl;

forAll(patches, patchi)
{

const fvPatch& thePatchItselfWall = patches[patchi];

if (isA<wallFvPatch>(thePatchItselfWall))

{
forAll (thePatchItselfWall,iFace)

{
label iCell = thePatchItselfWall.faceCells()[iFace];

vector iCellCentre = mesh.cellCentres()[iCell];

//pointed located at the centre line but at same axial position

vector bulkPoint(0,iCellCentre.y(),iCellCentre.z());

label bulkCell = mesh.findCell(bulkPoint);

if(error[iCell]>scalar(1.0))

{ count = count+1.0;

if(qtotal[iCell]>qwall.value())

{
Tnewn[iCell] = Tw[iCell];
}

else
{
Tnewp[iCell] = Tw[iCell];
}

Twnew[iCell] = ((Tnewn[iCell]+Tnewp[iCell])/2.0);

Tw[iCell] = Twnew[iCell];
}

else
{
Info<<"\nUpdating values. Wait";
}

}

}

}

if(count!=0)

{

goto label1;

}

else
{
Info<<"\nIterations Finished..Exiting the loop\n";

goto label2;
}

label2:

Info<<"Exiting the loop finally"<<endl;
...

While compiling it is not showing any error but when I am running a case...It is entering the loop once or twice and then showing floating point error!!!!!!!!
Can you point out the error?...If you want I can send you the part of code also...

Thanks a lot for helping me out!

regards,
Raj

 March 12, 2014, 09:04 #4 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,423 Rep Power: 25 Well, it'll surely be better if you just attach your code to a message. Right now (though the code is a mess, and inserted in the body of the message without CODE tag) I have the following questions: 1. What is Tb? 2. Why use Tb/Tb*Tnew instead of 0.0*Tb + Tnew, in this case you exclude the error in case if Tb == 0 at some point. 3. count should not be scalar as it is int. 4. Can you provide output with the floating point error you get?

March 12, 2014, 12:46
#5
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by alexeym Well, it'll surely be better if you just attach your code to a message. Right now (though the code is a mess, and inserted in the body of the message without CODE tag) I have the following questions: 1. What is Tb? 2. Why use Tb/Tb*Tnew instead of 0.0*Tb + Tnew, in this case you exclude the error in case if Tb == 0 at some point. 3. count should not be scalar as it is int. 4. Can you provide output with the floating point error you get?

Hi alexeym...
Sorry for the messy part. Tb is the liquid temperature and I have used it in 2 to make it Volume field as Tnew is a dimensionedScalar. Your point is absolutely right but nowhere temperature is zero, so I missed that part. Also, the output error is as follows..

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 Uninterpreted:
#3 __pow_finite in "/lib/i386-linux-gnu/libm.so.6"
#4 pow in "/lib/i386-linux-gnu/libm.so.6"
#5 Foam:ow(Foam::Field<double>&, Foam::UList<double> const&, double const&) at ??:?
#6 void Foam:ow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) at ??:?
#7 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:ow<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::dimensioned<double> const&) at ??:?
#8 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:ow<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, double const&) at ??:?
#9
at ??:?
#10 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#11
at ??:?
Floating point exception (core dumped)
Attached Files
 wallHeat.H (4.3 KB, 11 views)

 March 12, 2014, 12:51 #6 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,423 Rep Power: 25 In the output you've provided it is mentioned that the problem in in pow call. AFAIK in C standard fractional power of negative number throws FPE. So, my question is: can Tw be lower than Tsat?

March 12, 2014, 12:54
#7
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by alexeym In the output you've provided it is mentioned that the problem in in pow call. AFAIK in C standard fractional power of negative number throws FPE. So, my question is: can Tw be lower than Tsat?

Hi...
No Tw cannot be lower than Tsat.....

March 12, 2014, 13:01
#8
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by rajcfd Hi... No Tw cannot be lower than Tsat.....

Hi....I have made a C++ program for the same...I am attaching that also...This program is what I am trying to implement in my solver code
Attached Files
 wallHeat2.doc (5.9 KB, 9 views)

March 12, 2014, 17:28
#9
Senior Member

Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,423
Rep Power: 25
Well, you've got plenty of pow calls and this is a quote from c++98 standard (http://www.cplusplus.com/reference/cmath/pow/):

Quote:
 If x is finite negative and y is finite but not an integer value, it causes a domain error.
I.e. if you use pow(x, 2.0) and x is negative, you'll get the domain error. I don't know what is Dref, Ddep. Also pow(x, 2) is rather slow compared to x*x.

Can you change pow(x, 2) and pow(x, 3) to x*x and x*x*x, change pow(x, 1.8...) to something else and check if you get the same error?

March 13, 2014, 00:37
#10
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by alexeym Well, you've got plenty of pow calls and this is a quote from c++98 standard (http://www.cplusplus.com/reference/cmath/pow/): I.e. if you use pow(x, 2.0) and x is negative, you'll get the domain error. I don't know what is Dref, Ddep. Also pow(x, 2) is rather slow compared to x*x. Can you change pow(x, 2) and pow(x, 3) to x*x and x*x*x, change pow(x, 1.8...) to something else and check if you get the same error?

Hi alexeym...
You are absolutely right. Changing the pow(x,1.8..) to x*x is working. The error it was showing earlier is not coming but it is taking long time cacculating the Tw (wall temperature). So, I will look at the the algo again......

Thanks a lot for your help...Will get in touch with you again...

March 13, 2014, 03:36
#11
Senior Member

Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,423
Rep Power: 25
Quote:
 Originally Posted by rajcfd Changing the pow(x,1.8..) to x*x is working.
That means Tw - Tsat is negative somewhere (so Tw CAN be lower than Tsat ).

March 13, 2014, 04:44
#12
Member

Join Date: Apr 2013
Posts: 32
Rep Power: 5
Quote:
 Originally Posted by alexeym That means Tw - Tsat is negative somewhere (so Tw CAN be lower than Tsat ).

Hi....
Yeah Tw-Tsat is negative somewhere. But finally after the iterations it should come more than Tsat. Actually I am using binary method, so it should iterate between positive and negative values till it reaches the solution!....

I will have a look at the code and will be in touch with you...Thanks a lot for your incredible help.....:-)

regards,
Raj

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Amit Mathur OpenFOAM 16 October 6, 2013 11:09 anon_a OpenFOAM Programming & Development 6 May 24, 2012 04:10 tupe OpenFOAM 0 October 3, 2011 05:55 wwwjuventuscom OpenFOAM 2 January 4, 2011 20:01 alvin11 OpenFOAM 0 October 28, 2010 03:37

All times are GMT -4. The time now is 00:56.