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

help to initialize alpha1

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

Reply
 
LinkBack Thread Tools Display Modes
Old   July 25, 2011, 14:18
Default help to initialize alpha1
  #1
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
Hello, i wrote a function who initializa alpha1, based in sample from funkysetFields, but it is giving me floating point error when access to the center of the face.

Any idea??

//// inicializamos alpha1 antes de comenzar los cálculos

const pointField& pp = mesh.points();

forAll(mesh.cells(),cellI) {

const cell& cFaces = mesh.cells()[cellI];
scalar temp1=0;
scalar temp2=0;
scalar temp3=0;
scalar count=0;
forAll(cFaces, i)
{
const labelList& f = mesh.faces()[cFaces[i]];
label nPoints = f.size();
scalar maxz=-1000;
scalar minz=1000;
Info<< "cara leida " << endl;
for (label id = 0; id < nPoints; id++)
{
if (pp[f[id]].component(2) >= maxz) maxz = pp[f[id]].component(2);
if (pp[f[id]].component(2) <= minz) minz = pp[f[id]].component(2);
}
scalar projecZ = maxz-minz;
Info<< "projeccion leida " << endl;
vector nn = mesh.Cf()[mesh.cells()[cellI][cFaces[i]]];
Info<< "vector centros " << nn << endl;
Info<< "calculo proj1 " << projecZ << " z face " << nn.component(2) << endl;
temp1 += (nn.component(2) + 0.5 * projecZ); Info<< "temp1 " << temp1 << endl;
temp2 += (nn.component(2) - 0.5 * projecZ); Info<< "temp2 " << temp2 << endl;
temp3 += (0.5 - nn.component(2) / (projecZ+0.00000001)); Info<< "temp3 " << temp3 << endl;
count++;
}
Info<< "leidas caras " << endl;
temp1 /= count; temp2 /= count; temp3 /= count;
if (temp1 < 0 ) alpha1[cellI] = 1;
else if (temp2 > 0) alpha1[cellI] = 0;
else alpha1[cellI] = temp3;
Info<< "calculo alpha1 " << endl;
}

Info<< "escribo alpha1 " << endl;

alpha1.write();
pablodecastillo is offline   Reply With Quote

Old   July 26, 2011, 05:24
Default
  #2
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
Now it is working perfect, the solution was;

//// initia alpha1 before loop time

const pointField& pp = mesh.points();

forAll(mesh.cells(),cellI) {

const cell& cFaces = mesh.cells()[cellI];
scalar temp1=0;
scalar temp2=0;
scalar temp3=0;
scalar count=0;
scalar alphatemp=0;
forAll(cFaces, i)
{
const labelList& f = mesh.faces()[cFaces[i]];
label nPoints = f.size();
scalar maxz=-1000;
scalar minz=1000;

for (label id = 0; id < nPoints; id++)
{
if (pp[f[id]].component(2) >= maxz) maxz = pp[f[id]].component(2);
if (pp[f[id]].component(2) <= minz) minz = pp[f[id]].component(2);
}
scalar projecZ = maxz-minz;

vector nn = mesh.Cf()[cFaces[i]];

temp1 = (nn.component(2) + 0.5 * projecZ); //Info<< "temp1 " << temp1 << endl;
temp2 = (nn.component(2) - 0.5 * projecZ); //Info<< "temp2 " << temp2 << endl;
temp3 = (0.5 - nn.component(2) / (projecZ+0.00000001)); //Info<< "temp3 " << temp3 << endl;
if (temp1 < 0 ) alphatemp += 1;
else if (temp2 > 0) alphatemp += 0;
else alphatemp += temp3;
count++;
}

alpha1[cellI] = alphatemp/count;
}

Info<< "writing alpha1 " << endl;

alpha1.write();


So if you include this file before the loop you will avoid to use funkysetfields.

Pablo
pablodecastillo 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
Formulation in compressibleInterFoam scttmllr OpenFOAM Running, Solving & CFD 52 May 2, 2015 11:42
Problem with interFoam; Wave/wiggle alpha1 behavior JonW OpenFOAM 3 February 23, 2013 21:41
Implicit solver for gamma volume fraction equation gives alpha1 > 1. vitor.geraldes@ist.utl.pt OpenFOAM 0 April 19, 2010 08:29
Fixing the alpha1 min/max under/overshoots in interFoam.C solver. idrama OpenFOAM 1 January 27, 2010 20:34
alpha or alpha1 in interFoam & interDyMFoam wavytracy OpenFOAM Bugs 3 September 10, 2009 02:51


All times are GMT -4. The time now is 09:23.