CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Snippet for redefining fixedGradient boundary condition for a patch inside solver

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

Like Tree12Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 3, 2017, 19:36
Default Snippet for redefining fixedGradient boundary condition for a patch inside solver
  #1
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Fellows

I want to redifine fixedGradient boundary condition for a specified patch inside my source code. i.e. whenever the code sees "fixedGradient" in 0/cPlus file (cPlus is a volScalarField and a dependent variable)

My snippet now looks like:

Code:
const polyPatchList& patches = mesh.boundaryMesh();
   forAll (patches,patchI){
     
      if(cPlus.boundaryField()[patchI].type() == 'fixedGradient'){


           fixedGradientPatchField& fgPatchField
         (
          refCast<fixedGradientPatchField>(cPlus.boundaryField()[patchI])
         );
          fgPatchField.gradient() = fvc::snGrad(cPlus) + cPlus * fvc::snGrad(Ue);
       }
    }
In above snippet, I want to tell the code that whenever it saw "fixedGradient" for a patch inside 0/cPlus file calculate the following relation at the boundary surface and assign it as boundary condition:

\frac{d c^+}{d n} + c^+ \frac{d \varphi}{d n} in which c^+ & \varphi are dependent scalar variables that should be calculated at the boundary patch from previous time step.

However the code compilation returns this error:

Code:
createFields.H:501:3: error: ‘fgPatchField’ was not declared in this scope
In file included from interFoamEHDIon.C:79:0:
My first question is How should I declare fgPatchField ? since I tried volScalarField or volVectorField but did not work. Moreover it does not accept fvc::snGrad. Should I try fvc::reconstruct(snGrad(cPlus))?


The second error is:
Code:
createFields.H:532:12: error: ‘fixedGradientFvPatchVectorField’ was not declared in this scope

createFields.H:533:74: error: no matching function for call to ‘refCast(Foam::fvPatchField<Foam::Vector<double> >&)’
By the way, My OpenFoam version is 2.1.x

Last edited by babakflame; March 4, 2017 at 18:17.
babakflame is offline   Reply With Quote

Old   March 4, 2017, 19:09
Default
  #2
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Fellows

I made my snippet working. This snippet applies above formulation to a specific patch that has been predefined in 0/cPlus file with type fixedGradient. I haven't tested it though.

You need to add this header to your main.C file as well:
Code:
#include "fixedGradientFvPatchFields.H"
The snippet now looks like:

Code:
const polyPatchList& patches = mesh.boundaryMesh();
       forAll (patches,patchI){      

    if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){
        fixedGradientFvPatchScalarField& gradcPlusPatch =
            refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]);
        scalarField& gradcPlusField = gradcPlusPatch.gradient();
        gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();
     }
       }
I will keep ya updated with the results.


Keep Foaming Fellows.
Luttappy likes this.
babakflame is offline   Reply With Quote

Old   March 4, 2017, 19:27
Default
  #3
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 253
Rep Power: 18
Zeppo will become famous soon enough
Code:
fixedGradientPatchField& fgPatchField
(
    refCast<fixedGradientPatchField>(cPlus.boundaryField()[patchI])
);
fgPatchField.gradient() = fvc::snGrad(cPlus) + cPlus * fvc::snGrad(Ue);
->

Code:
fixedGradientFvPatchField<scalar>& fgPatchField
(
    refCast<fixedGradientFvPatchField<scalar> >(cPlus.boundaryField()[patchI])
);
fgPatchField.gradient() = 
    cPlus.boundaryField()[patchI].fvPatchField<scalar>::snGrad()
  + cPlus.boundaryField[patchI].snGrad() * Ue.boundaryField[patchI].snGrad();
Should be compiled but I am not sure the code corresponds to your formula. I think you'll have to develop a custom boundary condition.
babakflame likes this.

Last edited by Zeppo; March 4, 2017 at 19:30. Reason: oh.. you forestalled me
Zeppo is offline   Reply With Quote

Old   March 4, 2017, 19:41
Default
  #4
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Thanks Sergei

Actually, I have doubts on my final line in the snippet addressing formulation i.e.:

Code:
gradcPlusField = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();
I am gonna try the formulation of last line and give the feedbacks.

Regards

Last edited by babakflame; March 4, 2017 at 20:45.
babakflame is offline   Reply With Quote

Old   March 6, 2017, 04:34
Default
  #5
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 884
Rep Power: 28
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by babakflame View Post
Thanks Sergei

Actually, I have doubts on my final line in the snippet addressing formulation i.e.:

Code:
gradcPlusField = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();
I am gonna try the formulation of last line and give the feedbacks.

Regards
Hi babakflame,

What equation are you trying to implement as the boundary condition?

I don't think the following makes sense:
Code:
cPlus.boundaryField()[patchI].snGrad() = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();
which is essentially what you are doing.

This would only converge if "cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad()" is zero.

Maybe you are trying to set something like the following?
Code:
cPlus.boundaryField()[patchI].snGrad() = 
       - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();
Philip
bigphil is offline   Reply With Quote

Old   March 6, 2017, 18:08
Default
  #6
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Hi BigPhil

Many thanks for your reply.

In my domain, I have upper and lower walls which fluid is passing inbetween.

Actually I am trying to set a fixed flux (at lower wall) and zero flux (at upper wall) boundary conditions for scalar cPlus and zero flux (at lower wall) and fixed flux (at upper wall) boundary conditions for scalar cMinus.

My final boundary condition relations are as follows: (A_1 and B_1 are constants.) Actually, B_1 is the value of flux. However, these relations are the boundary conditions I have for cPlus and cMinus.
\left\{
\begin{array}{ll}
\frac{d c^+}{dn} = - A_1 c^+  \frac{d \varphi}{dn}\\
\frac{d c^-}{dn} = -  B_1 + A_1 c^- \frac{d \varphi}{dn} 
\end{array}
\right. for the upper wall

and the following relations:
\frac{d c^+}{dn} = -  B_1 - A_1 c^+  \frac{d \varphi}{dn}

\frac{d c^-}{dn} = - A_1 c^- \frac{d \varphi}{dn} for the lower wall.
n is the direction normal to the wall and outwards.

At the moment, I am trying this snippet in my main.C file:
(\varphi is denoted by Ue in my code which is another scalar (dependent variable)).

Code:
const label& ID= mesh.boundaryMesh().findPatchID("midup");
if(cPlus.boundaryField()[ID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& gradcPlusPatch =
		refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[ID]);
	scalarField& gradcPlusField = gradcPlusPatch.gradient();
	gradcPlusField = - eKbT * gradcPlusPatch * Ue.boundaryField()[ID].snGrad();
}
else if(cMinus.boundaryField()[ID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& gradcMinusPatch =
		refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[ID]);
	scalarField& gradcMinusField = gradcMinusPatch.gradient();
	gradcMinusField =  Dinv * 2.8e-11 + eKbT * gradcMinusPatch * Ue.boundaryField()[ID].snGrad();
}
       
const label& patchID = mesh.boundaryMesh().findPatchID("middown");
if(cPlus.boundaryField()[patchID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& cPlusPatch =
		refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchID]);
	scalarField& gradcPlusField = cPlusPatch.gradient();
	gradcPlusField =  Dinv * 2.8e-11 + eKbT * cPlusPatch * Ue.boundaryField()[patchID].snGrad();
	 }

else if(cMinus.boundaryField()[patchID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& cMinusPatch =
		refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchID]);
	scalarField& gradcMinusField = cMinusPatch.gradient();
		gradcMinusField = - eKbT * cMinusPatch * Ue.boundaryField()[patchID].snGrad();
}
Attached, you can see a picture of my domain with coordinates and n direction.

Another quick question:

I have defined the type of boundary condition as fixedGradient in 0/cPlus and 0/cMinus for upper and lower walls. However, the value of this gradient in 0 folder (I need to assign a value there like 0) still affects my results. When we use this method, the code should use my snippet as a custom boundary condition not the initial value in 0 folder, Isn't it?

Regards
Attached Images
File Type: png Cfd-Online Case.png (5.7 KB, 53 views)

Last edited by babakflame; March 6, 2017 at 21:29.
babakflame is offline   Reply With Quote

Old   March 7, 2017, 12:48
Default
  #7
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 253
Rep Power: 18
Zeppo will become famous soon enough
Quote:
Originally Posted by babakflame View Post
\left\{
\begin{array}{ll}
\frac{d c^+}{dn} = - A_1 c^+  \frac{d \varphi}{dn}\\
\frac{d c^-}{dn} = -  B_1 + A_1 c^- \frac{d \varphi}{dn} 
\end{array}
\right.
Now as we know what Bobi really wanted to get Philip's suggestion looks rather feasible
Quote:
Originally Posted by bigphil View Post
Maybe you are trying to set something like the following?
Code:
cPlus.boundaryField()[patchI].snGrad() = 
       - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();
Zhiheng Wang likes this.
Zeppo is offline   Reply With Quote

Old   March 7, 2017, 16:47
Default
  #8
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 884
Rep Power: 28
bigphil will become famous soon enoughbigphil will become famous soon enough
Hi babakflame,

A couple of comments/questions:
  • in the picture of your domain, the normal vectors "n" point inwards (i.e. in to the domain); however, OpenFOAM convention assumes the normals point outwards;
  • where do you place this code snippet? If you place it in createFields.H then the boundary gradient will only be set once at the start of the simulation e.g. dcPlusdn will be set based on the initial value of cPlus on the boundary, so if cPlus is zero at the start then dcPlusdn will also be set to zero. So maybe you want to put this code at the start of the time-loop so it is called once at the start of every time-step: this would set dcPlusdn based on the previous time-step value of cPlus (an explicit method), or if you need it more implicit then you would need to call this code multiple times per time-step (e.g. put a outer loop around your cPlus equation).

Philip
bigphil is offline   Reply With Quote

Old   March 8, 2017, 00:17
Default
  #9
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Philip

Many thanks for your reply.

As a respond to your comments:

  • I wasn't aware of the direction of n in OpenFoam in the way that you had mentioned, but by resimulating my case couple of times, I figured out that I need to multiply the lower wall boundary conditions by a negative. (This with your description now makes sense to me).
  • I put this snippet inside my time-loop before the pimple algorithm inner-loop. So, I think it updates the boundary conditions once at each time step.


I have a doubt here:
Although I have a gradient of \varphi inside my domain, I think the term \frac{d \varphi}{d n} within the boundary condition relations is not calculated (The code considers it as zero). Basically, its just the value of flux B_1 in my above boundary conditions that affects the gradient of c^+ or c^-.



Just for clarification, this is my final snippet:


Code:
const polyPatchList& patches = mesh.boundaryMesh(); 
    forAll (patches,patchI){   
 
    if(patches[patchI].name()=="midup"){     
 
         if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcPlusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
             scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
             gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();  // midup
         } 
 
         if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
              fixedGradientFvPatchScalarField& gradcMinusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
              scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
              gradcMinusField =  cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + Dinv * 2.8e-11;  //midup 
          } 
      }      
      else if(patches[patchI].name()=="middown"){ 
          if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcPlusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
             scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
             gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + Dinv * 2.8e-11;  // middown 
 
          } 
 
          if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcMinusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
              scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
              gradcMinusField =  cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() ;  //middown 
 
            } 
        } 
   }
A question: Is the way I have implemented c^+ \frac{d \varphi}{d n} or c^+ \frac{d \varphi}{d n} in above snippet correct or should I just use :


Code:
fvc::snGrad(Ue)
for \frac{d \varphi}{d n} term?




Regards
babakflame is offline   Reply With Quote

Old   March 8, 2017, 04:35
Default
  #10
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 884
Rep Power: 28
bigphil will become famous soon enoughbigphil will become famous soon enough
The snGrad at a patch is correctly calculated from:
Code:
Ue.boundaryField()[patchI].snGrad()
If you were to use:
Code:
fvc::snGrad(Ue)
this would give you the snGrad on every internal face and boundary face; you could then take the patch values; of course, in your case you only want the patch values so the first makes sense (note: they will both give the same values).

So, I am not sure if your method is working for you now but the general implementation seems OK; if the values are still not as expected, I would suggest printing them to the screen and figuring out whats wrong.

Philip
babakflame likes this.
bigphil is offline   Reply With Quote

Old   March 8, 2017, 11:21
Default
  #11
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Philip

Many thanks. I am gonna investigate term by term (in the way that you mentioned) to see exactly whats going on.

Regards
babakflame is offline   Reply With Quote

Old   May 19, 2017, 15:22
Default
  #12
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Fellows

As abovementioned, I have successfully implemented fixed flux boundary condition. However, after non-dimensionalizing my solver and boundary conditions, some scales are introduced and boundary conditions (for example for lower wall) becomes as follow


\frac {d c^-}{d n}= - \frac {e \varphi_0}{k_B T} c^- \frac{d \varphi}{d n}
Here, the coefficient \frac {e \varphi_0}{k_B T} has the value 0f 10^4 which causes my solver crashing.

If I ran the non-dimensional code with the non-dimensional boundary condition without this coefficient every thing is fine. I will post the error soon. However, is there any way to improve the solver performance?

What happens is that my courant number of non-dimensional code starts oscilating before divergence in initial time steps.

I am aware of some approaches such as:
  • grid refinement
  • increasing relTol to 0.1

to improve the performance of the solver, However, none of these worked for me.


I tried to increase that coefficient gradually, however, did not work. here is the error:


Code:
Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#4  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#5  Foam::fvMatrix<double>::solve() at ??:?
#6  
 at ??:?
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
Any suggestion or previous experience?

Last edited by babakflame; May 19, 2017 at 22:53.
babakflame is offline   Reply With Quote

Old   May 19, 2017, 22:52
Default
  #13
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Even underrelaxing the related equations as follows did not work.

Code:
relaxationFactors
  {
     fields
      {

      }
     equations
     {
          Ue               0.1;
         "(cPlus|cMinus)*" 0.1;
    }
  }
babakflame is offline   Reply With Quote

Old   May 22, 2017, 10:27
Default
  #14
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 884
Rep Power: 28
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by babakflame View Post
Even underrelaxing the related equations as follows did not work.

Code:
relaxationFactors
  {
     fields
      {

      }
     equations
     {
          Ue               0.1;
         "(cPlus|cMinus)*" 0.1;
    }
  }

Hi babakflame,

Can you post the code you are using?
This might help give an idea of how the stability of your method could be improved (also it will help others who might want to try something similar).

Philip
bigphil is offline   Reply With Quote

Old   May 25, 2017, 15:19
Default
  #15
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Philip

Many thanks for your help.

Here is the main part of code (.C) file:

Code:
int main(int argc, char *argv[]) 
{ 
    #include "setRootCase.H" 
    #include "createTime.H" 
    #include "createMesh.H" 
 
     pimpleControl pimple(mesh); 
 
    #include "initContinuityErrs.H" 
    #include "createFields.H" 
    #include "readTimeControls.H" 
    #include "correctPhi.H" 
    #include "CourantNo.H" 
    #include "setInitialDeltaT.H" 
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
 
    Info<< "\nStarting time loop\n" << endl; 
 
    while (runTime.run()) 
    { 
        #include "readTimeControls.H" 
        #include "CourantNo.H" 
        #include "alphaCourantNo.H" 
        #include "setDeltaT.H" 
 
        runTime++; 
 
        Info<< "Time = " << runTime.timeName() << nl << endl; 
 
        twoPhaseProperties.correct(); 
 
    const polyPatchList& patches = mesh.boundaryMesh(); 
           forAll (patches,patchI){   
 
        if(patches[patchI].name()=="midup"){     
 
            if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcPlusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
                scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
                gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();  // midup 
 
             } 
 
                  if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcMinusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
                scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
                gradcMinusField = cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;  // 1D case 
 
             } 
              }      
        else if(patches[patchI].name()=="middown"){ 
            if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcPlusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
                scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
                gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;  //channelBB case 
                                 
             } 
 
                  if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcMinusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
                scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
                gradcMinusField = - cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() ;  //middown  channelBB case 
                               
             } 
        } 
    } 
 
 
        #include "alphaEqnSubCycle.H" 
        #include "ElectricEqn.H" 
         //Info<< "ElectricEqn over " << nl << endl; 
        // --- Pressure-velocity PIMPLE corrector loop 
        while (pimple.loop()) 
        { 
         
            #include "UEqn.H" 
 
            // --- Pressure corrector loop 
            while (pimple.correct()) 
            { 
                #include "pEqn.H" 
            } 
 
            if (pimple.turbCorr()) 
            { 
                turbulence->correct(); 
            } 
        } 
 
        runTime.write(); 
 
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" 
            << "  ClockTime = " << runTime.elapsedClockTime() << " s" 
            << nl << endl; 
    } 
 
    Info<< "End\n" << endl; 
 
    return 0; 
}
And the ElectricEqn.H is as follows:

Code:
        surfaceScalarField cPlusFlux  = -DLU * fvc::interpolate(cPlus)*mesh.magSf()*fvc::snGrad(Ue);

        surfaceScalarField cMinusFlux =  DLU * fvc::interpolate(cMinus)*mesh.magSf()*fvc::snGrad(Ue);

 //           surfaceScalarField cPlusFlux  = -fvc::interpolate(sgm)*mesh.magSf()*fvc::snGrad(Ue);
            fvScalarMatrix cPlusEqn
            (
                fvm::ddt(cPlus)
          + fvm::div(phi, cPlus)
              + fvc::div(cPlusFlux)
           - fvm::laplacian(kappa, cPlus)    
            );
//            cPlusEqn.relax();
            cPlusEqn.solve();


            fvScalarMatrix cMinusEqn
            (
                fvm::ddt(cMinus)
          + fvm::div(phi, cMinus)
              + fvc::div(cMinusFlux)
           - fvm::laplacian(kappa, cMinus)    
            );

//        cMinusEqn.relax();
            cMinusEqn.solve();


        rhoE = (cPlus - cMinus);
            
              // solve
            fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE * gamma
            );
//            UeEqn.relax();
            UeEqn.solve();
As you can see above, there are three variables : cPlus, cMinus and Ue that are solved once in each pimple loop. There are two transport equations for cPlus and cMinus and one Poisson equation for Ue.
gamma, DLU, Kappa and eta are the non-dimensional groups (constants) appear as non-dimensionalization.
The UEqn.H is as follow:
Code:
 surfaceScalarField muEff
    (
        "muEff",
        twoPhaseProperties.muf()
      + fvc::interpolate(rho*turbulence->nut())
    );




    tmp1 = rhoE * E;

    tmp2 = rhoE * Eext;

       Ux = U.component(vector::X);

       rhoEUx = rhoE * Ux;

 
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muEff, U)
      - (fvc::grad(U) & fvc::grad(muEff))
      - tmp2 * eta
      - tmp1 * eta
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)

                ) * mesh.magSf()
            )

        );
    }
The settings I am using for fvSolutions and fvSchemes are as follows:

Code:
solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       0.001;
            relTol          0;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       0.0001;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.05;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration false;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-08;
            relTol          0;
            nVcycles        2;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-08;
        relTol          0;
        maxIter         20;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-09;
        relTol          0;
        nSweeps         1;
    }

    "(Ue|UeExt)"
    {
//        solver         PCG;
//        preconditioner   DIC;
//        tolerance        1e-14;
//        relTol           0.2;
        solver smoothSolver;
        smoother GaussSeidel;
        preconditioner DIC;
        tolerance 1e-14;
        relTol 0.1;
    };

//    rhoE
    "(cPlus|cMinus)"
    {
        solver         PBiCG;
        preconditioner   DILU;
        tolerance        1e-16;
        relTol           0.000;
    };


    "(k|B|nuTilda)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    nAlphaCorr      1;
    nAlphaSubCycles 3;
    cAlpha          1;

    pRefPoint      (0 5e-5 1e-5);
    pRefValue      0;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
/*         phi        0.1;
         phiFinal     1;
         cPlus      0.1;
         cPlusFinal   1;
         cMinus     0.1;
         cMinusFinal  1;
*/     }
}
Code:
ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default          Gauss linear;
    div(rho*phi,U)   Gauss upwind;
//    div(phi,rhoE)    Gauss upwind;
    div(phi,cMinus)    Gauss upwind;
    div(phi,cPlus)    Gauss upwind;
    div(phi,alpha)   Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
    rhoE;
}
As you can see in the main file (.C) one, for some patches I have defined fixedGradient which will be overwritten by my defined boundary condition.

for the Ue boundary condition on that patches, I want to use Dirichlet boundary condition, with high values, However, the code diverges, just for low values it will continue the simulation.


The code is based on interFoam in OpenFoam version 2.1.x.

Many thanks

Last edited by babakflame; May 25, 2017 at 19:03.
babakflame is offline   Reply With Quote

Old   May 26, 2017, 13:07
Default
  #16
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
As you can see above philip, the Ue variable is found through the Poisson equation (laplace with source term).

My boundary condition for Ue is Dirichlet type. I have tried the followings:

  • smooth solver (Gauss seidel with relTol = 0.1)
  • refining grid close to boundaries
  • relaxing the equation
  • higher order fv schemes
None of these worked. The initial residual of either cPlus or cMinus or Ue starts to grow and reaching to one which throws a floating point exception error.

Last edited by babakflame; May 26, 2017 at 17:14.
babakflame is offline   Reply With Quote

Old   May 26, 2017, 14:03
Default
  #17
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 884
Rep Power: 28
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by babakflame View Post
As you can see above philip, the phi variable is found through the poisson equation (laplacae with source term).

My boundary condition for phi is Dirichlet type. I have tried the followings:

  • smooth solver (Gauss seidel with relTol = 0.1)
  • refining grid close to boundaries
  • relaxing the equation
  • higher order fv schemes
None of these worked. The initial residual of either cPlus or cMinus or Ue starts to grow and reaching to one which throws a floating point exception error.
Hi babakflame,

Can you also post the output from the solver? What exactly is diverging? the linear solver for one of the equations? or the outer pimple loop? or time-steps?

From a quick look through the code you have posted, it seems that you are solving three coupled PDEs in a staggered manner; of course, convergence in these cases may not be easy to achieve using a staggered method (if the PDEs are strongly coupled).
A fruitful direction might be performing outer iterations over the equations with significant under-relaxation of the key variables and/or under-relaxation of the boundary conditions (I would start with the BCs). At least try figure out what exactly is causing the divergence.

Philip
bigphil is offline   Reply With Quote

Old   May 26, 2017, 15:00
Default
  #18
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Dear Philip

many thanks for your help.

First of all, these are boundary files in 0 folder:
cPlus
Code:
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      cPlus;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 -3 1 0 0 1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}
cMinus:

Code:
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      cMinus;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 -3 1 0 0 1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}
Code:
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      Ue;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 2 -3 0 0 -1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedValue;
        value           uniform 20;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedValue;
        value           uniform -20;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}
This is a channel flow case with these additional equations. As you can see in 0/Ue file only I am capable to run the case with values of Ue in order of 10 for midup and middown patches (middle part of channel). Even with this order of values in the middle of simulation, the divergence happens with the following error:

Code:
Courant Number mean: 5.03946621477e+51 max: 7.02462243817e+55
Interface Courant Number mean: 0 max: 0
deltaT = 3.52225240641e-108
Time = 2.16860254237690730505505598558

--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 30 to 31 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 31 to 32 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 32 to 33 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 33 to 34 to distinguish between timeNames at time 2.16860254238
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 2.62480833735e-11, Final residual = 1.57380692088e-25, No Iterations 1
DILUPBiCG:  Solving for cMinus, Initial residual = 1, Final residual = 1.16439623825e-18, No Iterations 3
smoothSolver:  Solving for Ue, Initial residual = 1, Final residual = 0.0999529656162, No Iterations 33
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  void Foam::dot<Foam::Vector<double>, Foam::Vector<double> >(Foam::Field<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type>&, Foam::UList<Foam::Vector<double> > const&, Foam::UList<Foam::Vector<double> > const&) at ??:?
#4  void Foam::dot<Foam::Vector<double>, Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#5  Foam::tmp<Foam::GeometricField<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type, Foam::fvPatchField, Foam::volMesh> > Foam::operator&<Foam::Vector<double>, Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6  
 at ??:?
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8  
 at ??:?
Floating point exception (core dumped)
This error is similar to the one that arises in higher orders of Ue in the very first time steps.

Regards
babakflame is offline   Reply With Quote

Old   May 26, 2017, 15:08
Default
  #19
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
This is the situation of residuals a few time steps before divergence:

Code:
Courant Number mean: 0.0135133107361 max: 0.0199948706422
Interface Courant Number mean: 0 max: 0
deltaT = 9.99750062485e-05
Time = 2.16858285429

MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 2.31041374006e-06, Final residual = 3.23788533007e-17, No Iterations 3
DILUPBiCG:  Solving for cMinus, Initial residual = 0.961364873856, Final residual = 3.50366344753e-18, No Iterations 5
smoothSolver:  Solving for Ue, Initial residual = 0.946454822981, Final residual = 0.0928075889037, No Iterations 14
GAMG:  Solving for p_rgh, Initial residual = 0.99287837977, Final residual = 0.0320233117181, No Iterations 3
time step continuity errors : sum local = 1.7866809096e-06, global = 4.76494332971e-10, cumulative = 4.7773692495e-10
GAMG:  Solving for p_rgh, Initial residual = 0.347493694957, Final residual = 0.0169072858213, No Iterations 2
time step continuity errors : sum local = 1.99215077851e-06, global = 1.3176994971e-09, cumulative = 1.79543642205e-09
GAMGPCG:  Solving for p_rgh, Initial residual = 0.251746573971, Final residual = 6.03393734406e-09, No Iterations 18
time step continuity errors : sum local = 7.01209128599e-13, global = 2.68332995554e-15, cumulative = 1.79543910538e-09
ExecutionTime = 1580.05 s  ClockTime = 1585 s

Courant Number mean: 0.0135540948435 max: 0.507795570766
Interface Courant Number mean: 0 max: 0
deltaT = 1.96873675724e-05
Time = 2.16860254165

MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 4.54849797817e-07, Final residual = 4.1721661111e-17, No Iterations 2
DILUPBiCG:  Solving for cMinus, Initial residual = 0.991818666534, Final residual = 3.3247608014e-20, No Iterations 5
smoothSolver:  Solving for Ue, Initial residual = 0.993844252023, Final residual = 0.0980642096096, No Iterations 10
GAMG:  Solving for p_rgh, Initial residual = 0.999978321915, Final residual = 0.0355731771938, No Iterations 3
time step continuity errors : sum local = 0.0332311172846, global = 6.40661804589e-06, cumulative = 6.40841348499e-06
GAMG:  Solving for p_rgh, Initial residual = 0.251482888814, Final residual = 0.0110229824325, No Iterations 3
time step continuity errors : sum local = 0.0239228239629, global = 2.60788641251e-06, cumulative = 9.0162998975e-06
GAMGPCG:  Solving for p_rgh, Initial residual = 0.145311847528, Final residual = 7.59825317339e-09, No Iterations 20
time step continuity errors : sum local = 1.75845002939e-08, global = -7.27185444404e-11, cumulative = 9.01622717895e-06
ExecutionTime = 1581.15 s  ClockTime = 1586
As you can see above, the intial residuals of p_rgh, cMinus and Ue start to reach 1 just before crashing.

Last edited by babakflame; May 26, 2017 at 16:31.
babakflame is offline   Reply With Quote

Old   May 26, 2017, 17:35
Default
  #20
Senior Member
 
Bobi
Join Date: Oct 2012
Location: Chicago, Illinois
Posts: 452
Rep Power: 12
babakflame is on a distinguished road
Quote:
Originally Posted by bigphil View Post
Hi babakflame,

A fruitful direction might be performing outer iterations over the equations with significant under-relaxation of the key variables and/or under-relaxation of the boundary conditions (I would start with the BCs). At least try figure out what exactly is causing the divergence.

Philip
I have applied under-relaxation for the 3 equations of cPlus, cMinus and Ue. But how can I under-relax boundary conditions??

If you meant Fields, I have applied that too with values around 0.1, no results.

Is there any syntax besides fields and equations (in fvSolutions file) that does under-relaxation in openFoam?
babakflame is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind turbine simulation Saturn CFX 58 July 3, 2020 02:13
Problem with cyclic boundaries in Openfoam 1.5 fs82 OpenFOAM 36 January 7, 2015 01:31
Error finding variable "THERMX" sunilpatil CFX 8 April 26, 2013 08:00
[blockMesh] Cyclic BC's: Possible face ordering problem? (Channel flow) sega OpenFOAM Meshing & Mesh Conversion 3 September 28, 2010 13:46
[Gmsh] Import gmsh msh to Foam adorean OpenFOAM Meshing & Mesh Conversion 24 April 27, 2005 09:19


All times are GMT -4. The time now is 15:27.