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/)
-   -   Updating Boundary Conditions Each Iteration (https://www.cfd-online.com/Forums/openfoam-solving/88762-updating-boundary-conditions-each-iteration.html)

thomas99 May 25, 2011 08:24

Updating Boundary Conditions Each Iteration
 
Hello All,

I have a situation where I would like to calculate some boundary conditions from the internal field. The problem is that, for each time step, the boundary conditions should be iteratively updated as part of the solution. I have searched the forum and noticed that many people have recommended the groovyBC or swak4Foam add ons for recalculating a boundary condition from the flow field but these tools seem to do this at the end of a time step and not as part of the iterative process. Can I use one of these tools to somehow calculate the BCs iteratively?

Alternatively, I noticed there is a BC built into openFoam called calculated that might work but I can't really find too much information about this in the forum and the tutorials don't really make it clear to me exactly how it works.

Does anyone have any idea how I can iteratively calculate boundary conditions?

Regards,
Tom

marupio May 25, 2011 09:02

I don't think groovyBC does iterative solutions. You may need to write your own boundary condition, and put the iterative solver into it. See

http://openfoamwiki.net/index.php/Ho...dary_condition

for an overview of writing new boundary conditions.

thomas99 May 25, 2011 09:31

Hi Marupio,

Thanks for the reply and the link.

I was hoping that I wouldn't have to do that but I'll give it a try.

Best,
Tom

bigphil May 25, 2011 09:35

Hi thomas99,

As marupio said you could write your own boundary condition,

alternatively, you could just include a header file which updates the boundary contact inside your solver solution loop,
for example I sometimes use the "fixedGradient" boundary condition in my solid stress solver to apply a traction, this means I must update this fixedGradient every solution iteration.
The header file might look something like the following (where U is the volVectorField I solve for):
Code:

label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest_name");
if
(
  U.boundaryField()[patchID].type()
  == fixedGradientFvPatchVectorField::typeName
)
{
          fixedGradientFvPatchVectorField& Upatch =
            refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);
         
          Upatch.gradient() = (.......calculate something here.....);
}


Hope it helps,
Philip

gschaider May 25, 2011 10:03

Quote:

Originally Posted by thomas99 (Post 309182)
Hello All,

I have a situation where I would like to calculate some boundary conditions from the internal field. The problem is that, for each time step, the boundary conditions should be iteratively updated as part of the solution. I have searched the forum and noticed that many people have recommended the groovyBC or swak4Foam add ons for recalculating a boundary condition from the flow field but these tools seem to do this at the end of a time step and not as part of the iterative process. Can I use one of these tools to somehow calculate the BCs iteratively?

Depends on what you mean by iterative. groovyBC is recalculated every time other BCs are recalculated too (sol also in a loop inside the time-loop)

Quote:

Originally Posted by thomas99 (Post 309182)
Alternatively, I noticed there is a BC built into openFoam called calculated that might work but I can't really find too much information about this in the forum and the tutorials don't really make it clear to me exactly how it works.

Forget this one. It means "I was calculated once, I won't change again. And I don't know if I'm Dirichlet or Neuman so I will behave like a Dirichlet"

thomas99 May 26, 2011 04:20

Hi Philip and Bernhard,

Thanks for the information. I think the inclusion of a header file is a good idea and I'll try that first.

Regards,
Tom

mathslw August 29, 2012 12:21

Quote:

Originally Posted by thomas99 (Post 309300)
Hi Philip and Bernhard,

Thanks for the information. I think the inclusion of a header file is a good idea and I'll try that first.

Regards,
Tom

Hi Tom,

Do you succeed in updating the boundary in the loop?
I want to update the boundary condition after each iteration for steady state calculation.
Thanks!

Wei

marupio August 29, 2012 12:26

I'm not infront of my FOAM machine right now, but I seem to recall an issue where the boundary condition would update only once per timestep. I was doing multiple iterations per timestep, and I had to call some function to force it to update. I forget the details... but if you encounter the situation where it only seems to update once, then let me know, I'll dig for the answer.

-Dave

mathslw August 29, 2012 12:30

Quote:

Originally Posted by marupio (Post 379349)
I'm not infront of my FOAM machine right now, but I seem to recall an issue where the boundary condition would update only once per timestep. I was doing multiple iterations per timestep, and I had to call some function to force it to update. I forget the details... but if you encounter the situation where it only seems to update once, then let me know, I'll dig for the answer.

-Dave

Hi Dave,

I am using the steady state solver, so I want to update the valure after each iteration.
Do you know about that?

Thanks!

Wei

mathslw August 29, 2012 13:08

Quote:

Originally Posted by bigphil (Post 309195)
Hi thomas99,

As marupio said you could write your own boundary condition,

alternatively, you could just include a header file which updates the boundary contact inside your solver solution loop,
for example I sometimes use the "fixedGradient" boundary condition in my solid stress solver to apply a traction, this means I must update this fixedGradient every solution iteration.
The header file might look something like the following (where U is the volVectorField I solve for):
Code:

label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest_name");
if
(
  U.boundaryField()[patchID].type()
  == fixedGradientFvPatchVectorField::typeName
)
{
          fixedGradientFvPatchVectorField& Upatch =
            refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);
         
          Upatch.gradient() = (.......calculate something here.....);
}


Hope it helps,
Philip

Hi Philip,

What should I do if I used the fixedValue and want to update it?
Thanks!

Wei

bigphil August 29, 2012 13:29

Quote:

Originally Posted by mathslw (Post 379357)
What should I do if I used the fixedValue and want to update it?
Thanks!

Hi Wei,

For a fixedValue boundary condition, you can do the following (U is the volVectorField which is solved):
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}

// set value on patch to what ever you want
U.boundaryField()[patchID] == vector(1,2,3);

Note the double equal signs, this is how you set the value on a patch. You don't need to cast the patch if it is fixedValue.

Hope it helps,
Philip

mathslw August 29, 2012 14:01

Quote:

Originally Posted by bigphil (Post 379363)
Hi Wei,

For a fixedValue boundary condition, you can do the following (U is the volVectorField which is solved):
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch to what ever you want
U.boundaryField()[patchID] == vector(1,2,3);

Note the double equal signs, this is how you set the value on a patch. You don't need to cast the patch if it is fixedValue.

Hope it helps,
Philip

Hi Philip,

Your suggestion is very helpful!

If I use U.boundaryField()[patchID] == vector(1,2,3) , the boundary would be uniform.
Can I use a loop to set non-uniform boundary conditions?
Besides, to set the new boundary, I have to make use of the velocity on the boundary face, velocity on the boundary cell and velocity gradient on the boundary cell. Do you know how to access these values?

Many thanks!

Wei

bigphil August 30, 2012 08:50

Hi Wei,

You can set non-uniform boundary values based on anything you want:
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
  vector faceOldVel = U.boundaryField()[patchID][facei];
  vector faceCellOldVel =
        U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
  vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
 
  // set whatever you want here
  U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip

mathslw August 30, 2012 12:01

Quote:

Originally Posted by bigphil (Post 379492)
Hi Wei,

You can set non-uniform boundary values based on anything you want:
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
  vector faceOldVel = U.boundaryField()[patchID][facei];
  vector faceCellOldVel =
        U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
  vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
 
  // set whatever you want here
  U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip

Hi Philip,

Thanks for your reply!
I will try it in my program to see how it is works!

Best Regards!

Wei

mathslw August 30, 2012 12:13

Quote:

Originally Posted by bigphil (Post 379492)
Hi Wei,

You can set non-uniform boundary values based on anything you want:
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
  vector faceOldVel = U.boundaryField()[patchID][facei];
  vector faceCellOldVel =
        U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
  vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
 
  // set whatever you want here
  U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip

Hi Philip,

As far as I know, the U.boundaryField()[patchID].snGrad()[facei] returns the patch-normal gradient. How to obtain the velocity gradient on the patch cell, which is a tensor. Can I define a
tensorField gradU = fvc::grad(U),
then use the
tensor faceCellGradientU =
gradU.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
to obtain the velocity gradient on the patch cell?

Many thanks!

Wei

mathslw September 9, 2012 17:19

Quote:

Originally Posted by bigphil (Post 379492)
Hi Wei,

You can set non-uniform boundary values based on anything you want:
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
  vector faceOldVel = U.boundaryField()[patchID][facei];
  vector faceCellOldVel =
        U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
  vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
 
  // set whatever you want here
  U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip

Hi Philip,

Sorry to bother u again.
If I use

U.boundaryField()[patchID] == vector(0.8, 0, 0);

I could update the boundary condition;
But if I use

forAll(U.boundaryField()[patchID],i)
{
U.boundaryField()[patchID][i] == vector(0.8, 0, 0);
}

The boundary would not be updated.

Do you know the reason?

Thanks!

Wei

bigphil September 10, 2012 10:03

Hi Wei,

To answer your first question, to calculate the full gradient (tensor) at the boundary face you use the fvc::grad operator as you have shown:
tensor faceCellGradientU =
gradU.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];

However, be careful because most fvc::grad schemes do not calculate the total gradient fully at boundary faces, they essentially calculate the normal gradient like snGrad and then they set the tangential gradient to be the same as the tangential gradient at the centre of the boundary cell.

As regards the forAll loop not updating the boundary I am not sure but here is a work-around:
Code:

vectorField newPatchValues(U.boundaryField()[patchID].size(), vector::zero);

forAll(U.boundaryField()[patchID],i)
{
    newPatchValues[i] = vector(0.8, 0, 0);
 }

U.boundaryField()[patchID][i] == newPatchValues;

Also, did you try chaining the "==" to "=" inside the forAll loop?

mathslw September 10, 2012 13:07

Quote:

Originally Posted by bigphil (Post 381068)
Hi Wei,

To answer your first question, to calculate the full gradient (tensor) at the boundary face you use the fvc::grad operator as you have shown:
tensor faceCellGradientU =
gradU.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];

However, be careful because most fvc::grad schemes do not calculate the total gradient fully at boundary faces, they essentially calculate the normal gradient like snGrad and then they set the tangential gradient to be the same as the tangential gradient at the centre of the boundary cell.

As regards the forAll loop not updating the boundary I am not sure but here is a work-around:
Code:

vectorField newPatchValues(U.boundaryField()[patchID].size(), vector::zero);
 
forAll(U.boundaryField()[patchID],i)
{
    newPatchValues[i] = vector(0.8, 0, 0);
 }
 
U.boundaryField()[patchID][i] == newPatchValues;

Also, did you try chaining the "==" to "=" inside the forAll loop?

Hi Philip,

You are right, I tried to change the "==" to "=" inside the forAll loop and it works.
In your code, you write "U.boundaryField()[patchID][i] == newPatchValues", is it U.boundaryField()[patchID]== newPatchValues? without the [i]?

Many thanks!

Wei

bigphil September 10, 2012 14:17

Quote:

Originally Posted by mathslw (Post 381094)
In your code, you write "U.boundaryField()[patchID][i] == newPatchValues", is it U.boundaryField()[patchID]== newPatchValues? without the [i]?

Yes, sorry that was a typo.

Best regards,
Philip

ancolli June 26, 2015 17:49

Quote:

Originally Posted by bigphil (Post 379492)
Hi Wei,

You can set non-uniform boundary values based on anything you want:
Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
  vector faceOldVel = U.boundaryField()[patchID][facei];
  vector faceCellOldVel =
        U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
  vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
 
  // set whatever you want here
  U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip


I have the following error when i try to compile the code for a scalar T:

DCP.C: In function ‘int main(int, char**)’:
DCP.C:75:66: error: no match for ‘operator[]’ (operand types are ‘Foam::tmp<Foam::Field<double> >’ and ‘Foam::label {aka int}’)
newPatchValues[i] = -T.boundaryField()[patchID].snGrad()[i];
^
make: *** [Make/linux64GccDPOpt/DCP.o] Error 1


any idea?

bigphil June 26, 2015 18:09

Quote:

Originally Posted by ancolli (Post 552323)
I have the following error when i try to compile the code for a scalar T:

DCP.C: In function ‘int main(int, char**)’:
DCP.C:75:66: error: no match for ‘operator[]’ (operand types are ‘Foam::tmp<Foam::Field<double> >’ and ‘Foam::label {aka int}’)
newPatchValues[i] = -T.boundaryField()[patchID].snGrad()[i];
^
make: *** [Make/linux64GccDPOpt/DCP.o] Error 1


any idea?

fvPatchField::snGrad() return a tmp copy of a scalarField (see OpenFOAM C++ doxygen documentation); so to get your code to compile, change this line:
Code:

      newPatchValues[i] = -T.boundaryField()[patchID].snGrad()[i];
to this:
Code:

      newPatchValues[i] = -T.boundaryField()[patchID].snGrad()()[i];
However, I suggest that you store this snGrad field before you loop over each face, so as the entire patch snGrad would not be re-calculated for every faces, e.g. something like:
Code:

const scalarField patchTSnGrad = T.boundaryField()[patchID].snGrad();
forAll(newPatchValues, i)
{
      newPatchValues[i] = -patchTSnGrad[i];
}

Philip

ancolli June 26, 2015 19:05

Thank you very much. I have some more questions.

1) If I need to perform mathematical operations (like log, exp, etc) to calculate new T from grad and old T. I have errors because of the different kind of variables. How can i solve it?

2) wich solver I have to select in order to verify the convergence of the BC?. Or the only change is modifying while (simple.correct Non Orthogonal())?

bigphil June 26, 2015 19:09

Quote:

Originally Posted by ancolli (Post 552337)
Thank you very much. I have some more questions.

1) If I need to perform mathematical operations (like log, exp, etc) to calculate new T from grad and old T. I have errors because of the different kind of variables. How can i solve it?

If you give the mathematics here that you are trying to implement and also your initial non-compiling implementation, then I can give you pointers on the implementation.

Quote:

Originally Posted by ancolli (Post 552337)
2) wich solver I have to select in order to verify the convergence of the BC?. Or the only change is modifying while (simple.correct Non Orthogonal())?

I'm not sure I understand your question; typically for verification, you find some benchmark (analytical or numerical) and compare your results.

Philip

ancolli June 26, 2015 20:41

Quote:

Originally Posted by bigphil (Post 552339)
If you give the mathematics here that you are trying to implement and also your initial non-compiling implementation, then I can give you pointers on the implementation.

could you please post some tutorials or web page?

Quote:

Originally Posted by bigphil (Post 552339)
I'm not sure I understand your question; typically for verification, you find some benchmark (analytical or numerical) and compare your results.

Philip

I am learning from laplacianFoam in steady state (the easy way). My question it is not referred to verification, instead I want to know an easy way to control convergence when BC change after each iteration. In laplacianFoam, the new nonlinear BC iterate only a number of times given by: simple.correct Non Orthogonal() setting in fvSolution

meth November 11, 2015 23:06

hi Pillip,

Your posts on updating boundary condition on each iteration was very helpful to me. I am implementing a new solver using the icoFoam solver to handle the fluid-solid interaction (Vortex induced vibration). The vibration of the solid is taken into the account by assuming the fluid is in a non-inertial frame vibrate on opposite to the solid motion while the solid stay still. The solid motion due to the fluid forces are solve inside the solver and update the boundary condition accordingly in each iteration.

My question:

Since fluid is in a non-inertial frame the acceleration of the frame needs to add into the navier stokes equations. For that I defined a field variable ddy. This accelaration is a uniform value across the fluid. I can update the boundary condition as you said. Can you please tell me how to update the internal field as well? :( I have defined the internal field initially as a uniform value.

Many thanks,

Methma

suman91 September 23, 2016 02:33

Hi,

This is an old thread but am trying to do the same thing, so posting my doubt here.

I am using OpenFOAM 2.4.0 and using pisoFoam solver. I want to change the fixedGradient B.C at each time loop. So I followed what was suggested here:

Code:

fixedGradientFvPatchVectorField& wallPatch = refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);

wallPatch.gradient() = dU_dy;


dU_dy is some "double" variable for which I calculate something.

But I get this error:
"fixedGradientFvPatchVectorField’ was not declared in this scope"

Please suggest where am I going wrong. Thank you.

gschaider September 23, 2016 06:02

Quote:

Originally Posted by suman91 (Post 619002)
Hi,

This is an old thread but am trying to do the same thing, so posting my doubt here.

I am using OpenFOAM 2.4.0 and using pisoFoam solver. I want to change the fixedGradient B.C at each time loop. So I followed what was suggested here:

Code:

fixedGradientFvPatchVectorField& wallPatch = refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);

wallPatch.gradient() = dU_dy;

dU_dy is some "double" variable for which I calculate something.

But I get this error:
"fixedGradientFvPatchVectorField’ was not declared in this scope"

Please suggest where am I going wrong. Thank you.

That's the technical wonder of Polymorphism at work: your source doesn't know the class but the executable does. Try including fixedGradientFvPatchFields.H

suman91 September 23, 2016 06:47

Hi Bernhard,

I have in included in myfile:

#include "fixedGradientFvPatchField.H"


Still am gettinhe the error. Just for your reference an pasting the start part of my file here. insertModification.H is where I am trying my thing.

Code:


#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "fixedGradientFvPatchField.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"

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

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

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readPISOControls.H"
        #include "CourantNo.H"

    #include "insertModification.H"

        // Pressure-velocity PISO corrector
        {
            // Momentum predictor

            fvVectorMatrix UEqn
            (
                fvm::ddt(U)
              + fvm::div(phi, U)
              + turbulence->divDevReff(U)
            );


gschaider September 23, 2016 07:55

Quote:

Originally Posted by suman91 (Post 619034)
Hi Bernhard,

I have in included in myfile:

#include "fixedGradientFvPatchField.H"


Still am gettinhe the error. Just for your reference an pasting the start part of my file here. insertModification.H is where I am trying my thing.

Fields.H not Field.H. The "s" is important

bigphil September 23, 2016 08:06

Quote:

Originally Posted by ancolli (Post 552323)
I have the following error when i try to compile the code for a scalar T:

DCP.C: In function ‘int main(int, char**)’:
DCP.C:75:66: error: no match for ‘operator[]’ (operand types are ‘Foam::tmp<Foam::Field<double> >’ and ‘Foam::label {aka int}’)
newPatchValues[i] = -T.boundaryField()[patchID].snGrad()[i];
^
make: *** [Make/linux64GccDPOpt/DCP.o] Error 1


any idea?

Hi,

Can you post your code segment?

I think the problem is "snGrad()" returns a "tmp" of a field, so either you modify the code to be:
Code:

      newPatchValues[i] = -T.boundaryField()[patchID].snGrad()()[i];
or a better solution (no need to re-generate the entire snGrad filed every iteration) would be to store the snGrad field:
Code:

const scalarField TsnGrad= T.boundaryField()[patchID].snGrad();
forAll(newPatchValues, i)
{
      newPatchValues[i] = -TsnGrad[i];
}


Philip

suman91 September 23, 2016 13:51

Hi Bernhard,
It is working. Was indeed a spell typo. Thanks a lot for your help.


Hi Philip,
It seems the error was because of a typo. Here is a section of my code:

Code:

label patchID = mesh.boundaryMesh().findPatchID("channelWall");
label nFaces = mesh.boundaryMesh()[patchID].size();
for(int facei = 0 ; facei<nFaces; facei++)
{
  // doing some calculations here ....

  const fixedGradientFvPatchVectorField& wallPatch = refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);
 
 wallPatch.gradient()[facei] == vector(0, dU_dy, 0);

Hoping this will set the gradient calculated for each face on the concerned patch.
Thank you.

suman91 September 26, 2016 09:19

Hi,

I have a doubt related to my last post.

Code:

const fixedGradientFvPatchVectorField& wallPatch =    refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);

for(int facei = 0 ; facei<nFaces; facei++)
{
  /**** 
    calculation of some gradient
  ***/ 

  Info<<" Gradient of face number "<< facei << "  Before changing is  "<<endl;
  Info<<    wallPatch.gradient()[facei]<<endl;

  wallPatch.gradient()[facei] == vector(du_dy, 0, 0);

  Info<<"  After changing is  "<<endl;
  Info<<    wallPatch.gradient()[facei]<<endl;

}

Problem is when I see the values in output, both values for before and after remain same. Here is a small section of the output

Code:

Face number  0
du_dy  -0.00022415
 Gradient of face number 0  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  1
du_dy  -0.00022415
 Gradient of face number 1  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  2
du_dy  -0.00022415
 Gradient of face number 2  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  3
du_dy  -0.00022415
 Gradient of face number 3  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)

Thank you in advance.

bigphil September 26, 2016 10:28

Quote:

Originally Posted by suman91 (Post 619300)
Hi,

I have a doubt related to my last post.

Code:

const fixedGradientFvPatchVectorField& wallPatch =    refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);

for(int facei = 0 ; facei<nFaces; facei++)
{
  /**** 
    calculation of some gradient
  ***/ 

  Info<<" Gradient of face number "<< facei << "  Before changing is  "<<endl;
  Info<<    wallPatch.gradient()[facei]<<endl;

  wallPatch.gradient()[facei] == vector(du_dy, 0, 0);

  Info<<"  After changing is  "<<endl;
  Info<<    wallPatch.gradient()[facei]<<endl;

}

Problem is when I see the values in output, both values for before and after remain same. Here is a small section of the output

Code:

Face number  0
du_dy  -0.00022415
 Gradient of face number 0  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  1
du_dy  -0.00022415
 Gradient of face number 1  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  2
du_dy  -0.00022415
 Gradient of face number 2  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)
Face number  3
du_dy  -0.00022415
 Gradient of face number 3  Before changing is 
(0 0 0)
  After changing is 
(0 0 0)

Thank you in advance.

You should not use the "==" operator; you should use the "=" operator to assign the face value.

The "==" operator is just necessary when assigning the value to an entire fixedValue-type boundary condition.

Philip

suman91 September 27, 2016 02:56

Hi Philip,

Thank you so much for the reply. I changed '==' to '=', then it gave me the following error:

Code:

error: passing ‘const Foam::Vector<double>’ as ‘this’ argument discards qualifiers [-fpermissive]
  wallPatch.gradient()[facei] = vector(du_dy, 0, 0);

A bit of searching and I found that declaring wallPatch as a 'const' was the problem. So I changed it to:

Code:

fixedGradientFvPatchVectorField& wallPatch = refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);
Now on printing, I see that the values are changing. Thanks a lot again for the reply.

Suman

suman91 October 12, 2016 13:17

Modified pisoFoam for pitzDaily - solver blowing up !!
 
2 Attachment(s)
Hi there,

Am having some problem in solving flow over backward facing step in OpenFOAM. I am using LES and Smagorinsky Model. At start it all goes well but after some time utau values at the wall starts increasing and the deltaT of simulation goes on decreasing (till if falls to 10^-12 or so order). I cannot detect what is the problem. I have modified the pisoFoam solver so impose stress boundary condition using the 'fixedGradient' feature of OpenFOAM.

Am attaching my required files here. Please have a look and suggest where am I going wrong. Will be waiting for a reply.

Attachment 51068

Attachment 51069

suman91 October 13, 2016 02:42

Hi there,

I think the problems are coming with the pressure solver settings. At t=0.01 pressure values are about -120 at the inlet and then increases linearly to close to zero at the outlet. But at T=0.05, it falls down to -4000 approx. Obviously in an incompressible flow, the absolute value of pressure does not matter, it is the gradient that matters. But what I am not understanding are such low values. Considering inlet velocity has a magnitude of 10 in the x direction (which is max amongst all 3 directions), a value of -50 to -70 may still be understandable but -4000, am not getting at all. Not understanding why the pressure field is behaving such. Maybe that is triggering the problem.

suman91 October 26, 2016 09:36

Hi,

Back again with an update on this same issue. So to start with, setting the gradient at the wall was not so much of a good idea it seems. Somehow the wall normal velocity was being calculated to be a non zero value at the wall. This might be one of the reasons for the solver to blow up in that manner. The modified code I am using now, uses fixedValue (0 0 0) as a boundary condition at the wall and in the modified solver, I calculate the velocity at the wall using the gradient and set the velocity. The solver now runs properly and does not blow up.

The only issue am facing with this is, the reynolds stress when calculated comes out to be very less; flow is not becoming turbulent only. Am wondering if it's an issue with the gradient schemes being used or not. If anyone has any views on this, it will be really helpful.

So the question is, what is the best gradient scheme to be used for solving incompressible turbulent flow (using LES) in openFOAM?

Sorry for maintaining two threads for the same topic. I will from now on post the updates on one particular thread itself:

Link: http://www.cfd-online.com/Forums/ope...tml#post623038

S.Colucci May 29, 2017 11:50

Dear Suman,
I have a similar problem, that is my piece of code
Quote:

PtrList<PtrList<volScalarField> > Index(yindex.size());

Index.set(0, new PtrList<volScalarField>(2));

Index[0].set
(
0,
new volScalarField
(
IOobject
(
"Index.pressure" ,
this->pair_.phase1().time().timeName(),
this->pair_.phase1().mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
Np_[0]*(p - one_p*pmin_[0])/( one_p*pmax_[0] - one_p*pmin_[0] )
)
);
// correct internalField
forAll(Index[0][0],celli)
{
Index[0][0][celli] = floor(Index[0][0][celli]);
}
// correct BC
forAll(Index[0][0].boundaryField(), patchi)
{

scalarField newPatchValues(Index[0][0].boundaryField()[patchi].size(),scalar(0));

forAll(Index[0][0].boundaryField()[patchi],i)
{
newPatchValues[i] = floor(Index[0][0].boundaryField()[patchi][i]);
}

Index[0][0].boundaryField()[patchi] == newPatchValues;
}
Commenting the last line
Quote:

Index[0][0].boundaryField()[patchi] == newPatchValues
it compiles and the values stored in newPatchValues are correct, so it seems that the problem is in assigning a new bc value to Index. That is the error:

error: no match for ‘operator==’ (operand types are ‘const Foam::fvPatchField<double>’ and ‘Foam::scalarField {aka Foam::Field<double>}’)
Index[0][0].boundaryField()[patchi] == newPatchValues;

If I substitute
Quote:

Index[0][0].boundaryField()[patchi] == newPatchValues
with
Quote:

Index[0][0].boundaryField()[patchi] = newPatchValues
I have the following error
error: passing ‘const Foam::fvPatchField<double>’ as ‘this’ argument of ‘void Foam::fvPatchField<Type>::operator=(const Foam::UList<T>&) [with Type = double]’ discards qualifiers [-fpermissive]
Index[0][0].boundaryField()[patchi] = newPatchValues;

It seems to be the same problem as you, but I do not have any const!
Any idea?
Simone

S.Colucci May 29, 2017 12:14

Quote:

Originally Posted by bigphil (Post 381098)
Yes, sorry that was a typo.

Best regards,
Philip

Dear Philip,
I post here the same problem that I posted in an other thread since here it is older but seems more appropriate

...

I have a similar problem, that is my piece of code
Quote:

PtrList<PtrList<volScalarField> > Index(yindex.size());

Index.set(0, new PtrList<volScalarField>(2));

Index[0].set
(
0,
new volScalarField
(
IOobject
(
"Index.pressure" ,
this->pair_.phase1().time().timeName(),
this->pair_.phase1().mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
Np_[0]*(p - one_p*pmin_[0])/( one_p*pmax_[0] - one_p*pmin_[0] )
)
);
// correct internalField
forAll(Index[0][0],celli)
{
Index[0][0][celli] = floor(Index[0][0][celli]);
}
// correct BC
forAll(Index[0][0].boundaryField(), patchi)
{

scalarField newPatchValues(Index[0][0].boundaryField()[patchi].size(),scalar(0));

forAll(Index[0][0].boundaryField()[patchi],i)
{
newPatchValues[i] = floor(Index[0][0].boundaryField()[patchi][i]);
}

Index[0][0].boundaryField()[patchi] == newPatchValues;
}
Commenting the last line
Quote:

Index[0][0].boundaryField()[patchi] == newPatchValues
it compiles and the values stored in newPatchValues are correct, so it seems that the problem is in assigning a new bc value to Index. That is the error:

error: no match for ‘operator==’ (operand types are ‘const Foam::fvPatchField<double>’ and ‘Foam::scalarField {aka Foam::Field<double>}’)
Index[0][0].boundaryField()[patchi] == newPatchValues;

If I substitute
Quote:

Index[0][0].boundaryField()[patchi] == newPatchValues
with
Quote:

Index[0][0].boundaryField()[patchi] = newPatchValues
I have the following error
error: passing ‘const Foam::fvPatchField<double>’ as ‘this’ argument of ‘void Foam::fvPatchField<Type>::operator=(const Foam::UList<T>&) [with Type = double]’ discards qualifiers [-fpermissive]
Index[0][0].boundaryField()[patchi] = newPatchValues;

It seems to be the same problem as you, but I do not have any const!
Any idea?
Simone

bigphil May 29, 2017 14:10

Hi Simone,

This error:
Code:

error: passing ‘const Foam::fvPatchField<double>’ as ‘this’ argument of ‘void Foam::fvPatchField<Type>::operator=(const Foam::UList<T>&) [with Type = double]’ discards qualifiers [-fpermissive]
is saying that this 'Index[0][0].boundaryField()[patchi]' is const.

From this code:
Code:

PtrList<PtrList<volScalarField> > Index(yindex.size());

Index.set(0, new PtrList<volScalarField>(2));

Index[0].set
(
0,
new volScalarField
(
IOobject
(
"Index.pressure" ,
this->pair_.phase1().time().timeName(),
this->pair_.phase1().mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
Np_[0]*(p - one_p*pmin_[0])/( one_p*pmax_[0] - one_p*pmin_[0] )
)
);
// correct internalField
forAll(Index[0][0],celli)
{
Index[0][0][celli] = floor(Index[0][0][celli]);
}
// correct BC
forAll(Index[0][0].boundaryField(), patchi)
{

scalarField newPatchValues(Index[0][0].boundaryField()[patchi].size(),scalar(0));

forAll(Index[0][0].boundaryField()[patchi],i)
{
newPatchValues[i] = floor(Index[0][0].boundaryField()[patchi][i]);
}

Index[0][0].boundaryField()[patchi] == newPatchValues;
}

it should not be const, and the code should compile; it compiles for me in foam-extend-32 (I would expect the same in other version too); of course I had to make some unrelated changes where I replace "yindex.size" with "1", "Np_[0]*(p - one_p*pmin_[0])/( one_p*pmax_[0] - one_p*pmin_[0] )" with "mag(U)", "this->pair_.phase1().time()" with "runTime" and "this->pair_.phase1().mesh()" with "mesh".

Are you sure this is the code you are compiling? did you, for example, post "simplified" code? If possible I would suggest making a very simple OpenFOAM utility that shows this compilation problem and then upload this for others to try.

Philip


All times are GMT -4. The time now is 13:33.