Snippet for redefining fixedGradient boundary condition for a patch inside solver

bigphil · May 27, 2017, 11:08

Quote:

Originally Posted by babakflame

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,

You can manually apply under-relaxation when setting the boundary conditions e.g.

Code:

// No under-relaxation
//gradcMinusField = cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;

// With under-relaxation
const scalarField newGrad = cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;
const scalar relaxFactor = 0.1;
gradcMinusField = relaxFactor*newGrad + (1.0 - relaxFactor)* gradcMinusField;

I have hard-coded the under-relaxation factor above but you could read this from a dict too, if desired.

Philip

babakflame · May 27, 2017, 16:58

Thanks Philip

The idea of coupled solver (coupling these 3 equations before pimple loop) is interesting for me too. I found out that Foam-Extend has some coupled solver examples. I will proceed through BC under-relaxation and coupled solver. Lets see how they perform.

Regards

babakflame · June 2, 2017, 21:51

Dear Philip
The idea of under-relaxing boundary conditions helped in the context of stabilization. Although, still the solution is unstable in the range of Ue greater than 10. But at least helped me to increase the boundary values of Ue up to 10. This is the snippet I am using to under-relax a fixedValue boundary condition:

Code:

                     const scalar& relaxFactor = 0.00001; 
                const scalarField& Ueold = Ue.boundaryField()[patchI]; 
                Ue.boundaryField()[patchI] == (1 - relaxFactor) * Ueold + relaxFactor * 90;

Even I put all those three equations within the pimple loop, so as to be updated in each iteration of pimple loop (to increase stability) I mean:

Code:

       #include "alphaEqnSubCycle.H" 
        // --- Pressure-velocity PIMPLE corrector loop 
        while (pimple.loop()) 
        { 
         
            #include "UEqn.H" 
            #include "ElectricEqn.H" 
            #include "SourceTerm.H" 
 
            // --- Pressure corrector loop 
            while (pimple.correct()) 
            { 
                #include "pEqn.H" 
            } 
 
            if (pimple.turbCorr()) 
            { 
                turbulence->correct(); 
            } 
        } 
 
        runTime.write();

I am gonna make the first two equations to be solved in a blockcoupled manner:

Code:

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

        surfaceScalarField cMinusFlux =  DeKbT*fvc::interpolate(cMinus)*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();

Then, the third one in a segregated manner (as before):

Code:

  fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -(cPlus - cMinus) * NA * elem
            );
        
            UeEqn.relax();
            UeEqn.solve();

Do you think its a reasonable approach?

bigphil · June 6, 2017, 12:53

Hi Bobi,

The block coupled method will help the convergence of the cPlus and cMinus equations; however, it seems that currently you are limited by the convergence of the Ue equation (because of the Ue boundary conditions), so I'm not sure of coupling cMinus and cPlus will have any effect.

I am not familiar with the equations/physics being solved, but personally I would try find how others have solved them (do they use coupled or staggered method, or maybe they use explicit methods or solve in the frequency domain...).

Depending on the strength of the coupling between the different equations, it may be beneficial to resolve the coupling between two of the equations first before solving the third e.g. loop over cPlus and cMinus until converged and then solve Ue in an outer loop.

Philip

babakflame · June 6, 2017, 16:43

Dear Philip

Many thanks for your reply. The idea of looping over cPlus and cMinus equations and then get to Ue equation looks interesting. However, I have a question here

, maybe I am missing sth. These are my equations (all three):
$\frac{\partial c^+}{\partial t} + \nabla \cdot \left( c^{+} V \right) = - \nabla \cdot (c^+ \nabla U_e) + D \nabla^2 c^{+}$

$\frac{\partial c^-}{\partial t} + \nabla \cdot \left( c^{-} V \right) = + \nabla \cdot (c^- \nabla U_e) + D \nabla^2 c^{-}$

$\nabla^2 U_e = - (c^{+} - c^{-})$

As you can see, i have put the header file including these equations inside pimple loop:

Code:

 while (pimple.loop()) 
        { 
         
            #include "UEqn.H" 
             
            #include "ElectricEqn.H" 
            #include "SourceTerm.H" 
            // --- Pressure corrector loop 
            while (pimple.correct()) 
            { 
                #include "pEqn.H" 
            } 
 
            if (pimple.turbCorr()) 
            { 
                turbulence->correct(); 
            } 
        }

Since $c^{+}$ and $c^{-}$ has a convective term, i.e. : $c^{+} V$ .

My Header file (ElectricEqn.H) is:

Code:

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

        surfaceScalarField cMinusFlux =  DeKbT*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) * NA * elem;
            
              // solve
            fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE
            );
            UeEqn.relax();
            UeEqn.solve();

As far as I know, when each equation is solved, it reaches to the tolerance mentioned in fvSolutions file. So, If I want to put for example, cPlus and cMinus in an inner loop in above header file, what would be the convergence criterion?

That's the point I don't get.

I mean even in solving the equations once, we have satisfied the tolerance speciified in fvSolutions file through the linear matrix solver.

Unfortunately, no body was interested in charge equations with super high potential values.

babakflame · June 6, 2017, 18:52

Dear Philip,

This is the snippet that I can think of:

Code:

const scalar& convergenceTolerance = 0.01;
scalar& initialResidualA = 0.0;
scalar& initialResidualB = 0.0;
scalar& initialResidual = 0.0;


do 
    {
        surfaceScalarField cPlusFlux  = -DeKbT*fvc::interpolate(cPlus)*mesh.magSf()*fvc::snGrad(Ue);

        surfaceScalarField cMinusFlux =  DeKbT*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)    
            );
            cPlus.relax();
//            cPlusEqn.solve();
            initialResidualA = cPlusEqn.solve().initialResidual();


            fvScalarMatrix cMinusEqn
            (
                fvm::ddt(cMinus)
          + fvm::div(phi, cMinus)
              + fvc::div(cMinusFlux)
           - fvm::laplacian(kappa, cMinus)    
            );
            cMinus.relax();
//            cMinusEqn.solve();
          initialResidualB = cMinusEqn.solve().initialResidual();

            initialResidual = (initialResidualB > initialResidualA) ? initialResidualB : initialResidualA;


} while (initialResidual > convergenceTolerance);

babakflame · June 6, 2017, 20:35

Unfortunately, above snippet did not help.
It seems that although it gaurantees that in each time step initial residual is less than a specific value, however, after couple of iterations, it reaches to a point that can not reduce initial residual and again divergence happens.

arjun · June 6, 2017, 23:03

Quote:

Originally Posted by babakflame

Unfortunately, above snippet did not help.
It seems that although it gaurantees that in each time step initial residual is less than a specific value, however, after couple of iterations, it reaches to a point that can not reduce initial residual and again divergence happens.

How are you under relaxing the last Poisson equation.

babakflame · June 6, 2017, 23:22

Hi Arjun

At first, I applied the under-relaxation to the matrix of coefficients (under-relaxing equation) as follows:

Code:

 fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE
            );
            UeEqn.relax();
            UeEqn.solve();

Then I applied the under-relaxation to the field as:

Code:

 fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE
            );
            Ue.relax();
            UeEqn.solve();

For the field, I am using lower under-relaxation factors (0.01) compared to 0.1 for the equation.

These are the values for the fields:

Code:

relaxationFactors
{
    fields
    {
        Ue         0.01;
        cMinus     0.001;
        cPlus      0.001;
    }
    equations
    {
     }
}

Moreover, there is a snippet that applies under-relaxation to the specified patch as follows:

Code:

const polyPatchList& patches = mesh.boundaryMesh(); 
           forAll (patches,patchI){   
 
        if(patches[patchI].name()=="midup"){     

                const scalar& relaxFactor = 0.0001; 
                const scalarField& Ueold = Ue.boundaryField()[patchI]; 
                Ue.boundaryField()[patchI] == (1 - relaxFactor) * Ueold + relaxFactor * 2; 
           }

It makes the patch to reach to a high value (here 200) by under-relaxation ( adding a value through geometric progression in each time step).

arjun · June 7, 2017, 00:20

The reason why I asked this question is that you can not apply implicit URF to Laplacian. You can only have explicit urf. By applying implicit URF to last equation it is not really converging for you.

What you should do is to use urf = 1 for this equation.

Personally

I would not use coupled solver for it, I would first solve cplus and cminus (may be coupled or segregated).
Then I would solve the Laplacian with explicit urf. (Not implicit).

Note that the situation here is quite similar to flow solver where momentum takes implicit urf but pressure takes explicit urf (because it is Poisson).

In your situation even urf = 0.95 does not work (only way is to use urf = 1)..

PS (edited to add): In case of Stokes flow where momentum is also Poisson then implicit urf is not applied and there we again apply explicit urf and it is difficult to solve.

babakflame · June 7, 2017, 10:29

Dear Arjun

Thanks for your comment. So, basically, what you recommend is:

Code:

fvScalarMatrix UeEqn
            (
                fvc::laplacian(eps,Ue) == -rhoE
            );
            UeEqn.solve();

with no no urf (urf=1). I am gonna try it. But I think previously, just with Poisson (no cPlus and cMinus) explicit discretizing was giving me lower stability.

I am gonna try it and give the feedback.

Regards

arjun · June 7, 2017, 11:37

Quote:

Originally Posted by babakflame

Dear Arjun

Thanks for your comment. So, basically, what you recommend is:

Code:

fvScalarMatrix UeEqn
            (
                fvc::laplacian(eps,Ue) == -rhoE
            );
            UeEqn.solve();

with no no urf (urf=1). I am gonna try it. But I think previously, just with Poisson (no cPlus and cMinus) explicit discretizing was giving me lower stability.

I am gonna try it and give the feedback.

Regards

Yes, we do not apply implicit urf to Poisson part.

We use explicit urf to make it stable. This is why it is better to take them out of the loop or if it is coupled then use no urf.

And yes, it will be less stable with coupled (contrary to what most believe).

babakflame · June 8, 2017, 21:42

Arjun,

I realized that explicit discretizing of Poisson equation in my case leads to immediate divergence.

Regards

arjun · June 9, 2017, 00:58

Quote:

Originally Posted by babakflame

Arjun,

I realized that explicit discretizing of Poisson equation in my case leads to immediate divergence.

Regards

Let me come upto speed on this one.

option (a) you are solving all three equations in one coupled system?

option (b) only 2 equations are solved in coupled manner but poisson equation was solved separately after cplus and cminus were solved.

So which one it is? For option (a) I expect immediate divergence. For option (b) the divergence depend on explicit under relaxation. (this you seem not to be using).

PS: Is it possible for you to show me what exactly we are solving here because you are solving flow, turbulence etc. So how exactly cplus and cminus are used along with flow? Or are they part of some other model.

babakflame · June 10, 2017, 20:01

Arjun,

Thanks for your reply. Here, I am interested in tracking ions in laminar fluids. Two transport equations are solved (3 unkowns incluing anions, cations and electric potential) plus a Poisson equation (the 3rd required equation).

The fluid flow is laminar. I have mentioned equations in detail in this post:
Snippet for redefining fixedGradient boundary condition for a patch inside solver

What I am doing is :

solving all above 3 equations in a loop (in segregated manner) to assure convergence in each time step:
(shows more stability, however, still diverges)

Code:

const scalar& convergenceTolerance = 0.01;
double initialResidual = 0.0;
double initialResidualA = 0.0;
double initialResidualB = 0.0;
double initialResidualC = 0.0;

do 
    {
        surfaceScalarField cPlusFlux  = -DeKbT*fvc::interpolate(cPlus)*mesh.magSf()*fvc::snGrad(Ue);

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

            fvScalarMatrix cPlusEqn
            (
                fvm::ddt(cPlus)
          + fvm::div(phi, cPlus)
              + fvc::div(cPlusFlux)
           - fvm::laplacian(kappa, cPlus)    
            );

            initialResidualA = cPlusEqn.solve().initialResidual();


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

            initialResidualB = cMinusEqn.solve().initialResidual();

            initialResidual = (initialResidualB > initialResidualA) ? initialResidualB : initialResidualA;



        rhoE = (cPlus - cMinus) * NA * elem;
            
              // solve
            fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE
            );


         initialResidualC = UeEqn.solve().initialResidual();

            initialResidual = (initialResidualC > initialResidual) ? initialResidualC : initialResidual;

} while (initialResidual > convergenceTolerance);

Still, I have not gone toward coupled solver.

anperri · August 9, 2017, 11:49

Hey I just wanted to attach a nice thread if you want to use the refCast snippet in a newer version of OF:

Quote:

Originally Posted by mboesi

since it took me some time to cast away the const, I thought it would be a good idea to post the solution here. In my case this worked to remove the constness:

Code:

fixedValueFvPatchVectorField& Upatch = 
    refCast<fixedValueFvPatchVectorField>(U.boundaryFieldRef()[patchi]);

calculatedFvPatchVectorField& Uhpatch = 
    refCast<calculatedFvPatchVectorField>(Uh.boundaryFieldRef()[patchi])

Thread:
of4 refCast complains

AnthonyP · August 9, 2017, 11:58

Hey I just wanted to attach a nice thread if you want to use the refCast snippet in a newer version of OF:

Quote:

Originally Posted by mboesi

since it took me some time to cast away the const, I thought it would be a good idea to post the solution here. In my case this worked to remove the constness:

Code:

fixedValueFvPatchVectorField& Upatch = 
    refCast<fixedValueFvPatchVectorField>(U.boundaryFieldRef()[patchi]);

calculatedFvPatchVectorField& Uhpatch = 
    refCast<calculatedFvPatchVectorField>(Uh.boundaryFieldRef()[patchi])

Thread:
of4 refCast complains

babakflame · March 30, 2018, 21:28

Dear Fellows

I have a question with regard to this custom boundary conditions as attached, any hint is appreciated.

babakflame · April 1, 2018, 22:51

Dear Fellows

I have written the following boundary condition:

$\frac{d c^+}{dy} = - j^+- c^+ \frac{d \varphi}{dy}$

in OpenFOAM as:

Code:



if(cPlus.boundaryField()[ID].type()=="fixedGradient"){
  fixedGradientFvPatchScalarField& gradcPlusPatch =
  refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[ID]);
  scalarField& gradcPlusField = gradcPlusPatch.gradient();
  gradcPlusField = -cPlus.boundaryField()[patchI]* Ue.boundaryField()[ID].snGrad()+jPlus;
  }

My coordinate system (y) is upwards and perpendicular to the boundary patch. As shown in the attached figure.

However, this boundary condition is based on normal vector (j+) and normal gradients
$\frac{d c^+}{dn}$ & $c^+ \frac{d \varphi}{dn}$

Now, my question is when y is upwards and also inwards; I know that OpenFOAM calculates the gradients outwards (i.e.

in my figure). What happens to the normal vector (here j+) which I know that has an outwards direction.
or in simple words, here in the figure y is upwards, and j+ outwards, should I put a positive value or negative value here in this boundary condition in OpenFOAM (the last line) when my boundary condition is $\frac{d c^+}{dy} = - j^+- c^+ \frac{d \varphi}{dy}$
?

Code:


gradcPlusField= -cPlus.boundaryField()[patchI]*Ue.boundaryField()[ID].snGrad()+jplus;

babakflame · April 5, 2018, 15:03

Dear Fellows

I want to make a slight change to this boundary to not be applied when the value at the patch becomes negative (go back to zero gradient).

Before this change, the working boundary was:

Code:

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

  if(patches[patchI].name()=="midup"){    
       if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){
  fixedGradientFvPatchScalarField& gradcMinusPatch =refCast<fixedGradientFvPatchScalarField>
(cMinus.boundaryField()[patchI][i]);
   scalarField& gradcMinusField = gradcMinusPatch.gradient();
   gradcMinusField = eKbT * cMinus.boundaryField()[patchI][i] * Ue.boundaryField()[patchI][i].snGrad()+Dinv*currentj; 
   }
  }

Now, I want to add if the variable is less than 0, recast zero gradient, otherwise do as before.

I can think of sth like this:

Code:

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

    if(patches[patchI].name()=="midup"){  
        forAll (patchI, facei){                                           
                        if (cMinus.boundaryField()[patchI][facei] < 0)
                        {
                        fixedGradientFvPatchScalarField& gradcMinusPatch =
                                refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI][facei]);
                        scalarField& gradcMinusField = gradcMinusPatch.gradient();
                        gradcMinusField = 0;
                        }
                        else if(cMinus.boundaryField()[patchI][facei] > 0)
                        {
              if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){
                             
            fixedGradientFvPatchScalarField& gradcMinusPatch =
                refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI][facei]);
            scalarField& gradcMinusField = gradcMinusPatch.gradient();
            gradcMinusField = eKbT * cMinus.boundaryField()[patchI][facei] * Ue.boundaryField()[patchI][facei].snGrad() + Dinv *  currentj; //3.13e13;   //midup
               }
               }}        
     }     
}

However, it gives this error:

Code:

/home/babak/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/typeInfo.H:110:35: error: cannot dynamic_cast ‘r’ (of type ‘double’) to type ‘class Foam::fixedGradientFvPatchField<double>&’ (source is not of class type)         return dynamic_cast<To&>(r);                                   ^/home/babak/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/typeInfo.H:115:40: error: request for member ‘type’ in ‘r’, which is of non-class type ‘double’             << "Attempt to cast type " << r.type()                                        ^/home/babak/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/typeInfo.H:119:35: error: cannot dynamic_cast ‘r’ (of type ‘double’) to type ‘class Foam::fixedGradientFvPatchField<double>&’ (source is not of class type)         return dynamic_cast<To&>(r);                                   ^

can anybody hint me how can I add this new condition based on patch value to the previously written snippet?

Regards

June 6, 2017, 12:53		#24
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Hi Bobi, The block coupled method will help the convergence of the cPlus and cMinus equations; however, it seems that currently you are limited by the convergence of the Ue equation (because of the Ue boundary conditions), so I'm not sure of coupling cMinus and cPlus will have any effect. I am not familiar with the equations/physics being solved, but personally I would try find how others have solved them (do they use coupled or staggered method, or maybe they use explicit methods or solve in the frequency domain...). Depending on the strength of the coupling between the different equations, it may be beneficial to resolve the coupling between two of the equations first before solving the third e.g. loop over cPlus and cMinus until converged and then solve Ue in an outer loop. Philip babakflame likes this.

June 6, 2017, 16:43		#25
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Philip Many thanks for your reply. The idea of looping over cPlus and cMinus equations and then get to Ue equation looks interesting. However, I have a question here, maybe I am missing sth. These are my equations (all three): $\frac{\partial c^+}{\partial t} + \nabla \cdot \left( c^{+} V \right) = - \nabla \cdot (c^+ \nabla U_e) + D \nabla^2 c^{+}$ $\frac{\partial c^-}{\partial t} + \nabla \cdot \left( c^{-} V \right) = + \nabla \cdot (c^- \nabla U_e) + D \nabla^2 c^{-}$ $\nabla^2 U_e = - (c^{+} - c^{-})$ As you can see, i have put the header file including these equations inside pimple loop: Code: while (pimple.loop()) { #include "UEqn.H" #include "ElectricEqn.H" #include "SourceTerm.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence->correct(); } } Since $c^{+}$ and $c^{-}$ has a convective term, i.e. : $c^{+} V$ . My Header file (ElectricEqn.H) is: Code: surfaceScalarField cPlusFlux = -DeKbTfvc::interpolate(cPlus)mesh.magSf()fvc::snGrad(Ue); surfaceScalarField cMinusFlux = DeKbTfvc::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) * NA * elem; // solve fvScalarMatrix UeEqn ( fvm::laplacian(eps,Ue) == -rhoE ); UeEqn.relax(); UeEqn.solve(); As far as I know, when each equation is solved, it reaches to the tolerance mentioned in fvSolutions file. So, If I want to put for example, cPlus and cMinus in an inner loop in above header file, what would be the convergence criterion? That's the point I don't get. I mean even in solving the equations once, we have satisfied the tolerance speciified in fvSolutions file through the linear matrix solver. Unfortunately, no body was interested in charge equations with super high potential values.

June 6, 2017, 23:22		#29
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Hi Arjun At first, I applied the under-relaxation to the matrix of coefficients (under-relaxing equation) as follows: Code: fvScalarMatrix UeEqn ( fvm::laplacian(eps,Ue) == -rhoE ); UeEqn.relax(); UeEqn.solve(); Then I applied the under-relaxation to the field as: Code: fvScalarMatrix UeEqn ( fvm::laplacian(eps,Ue) == -rhoE ); Ue.relax(); UeEqn.solve(); For the field, I am using lower under-relaxation factors (0.01) compared to 0.1 for the equation. These are the values for the fields: Code: relaxationFactors { fields { Ue 0.01; cMinus 0.001; cPlus 0.001; } equations { } } Moreover, there is a snippet that applies under-relaxation to the specified patch as follows: Code: const polyPatchList& patches = mesh.boundaryMesh(); forAll (patches,patchI){ if(patches[patchI].name()=="midup"){ const scalar& relaxFactor = 0.0001; const scalarField& Ueold = Ue.boundaryField()[patchI]; Ue.boundaryField()[patchI] == (1 - relaxFactor) * Ueold + relaxFactor * 2; } It makes the patch to reach to a high value (here 200) by under-relaxation ( adding a value through geometric progression in each time step).

June 7, 2017, 00:20		#30
arjun Senior Member Arjun Join Date: Mar 2009 Location: Nurenberg, Germany Posts: 1,274 Rep Power: 34	The reason why I asked this question is that you can not apply implicit URF to Laplacian. You can only have explicit urf. By applying implicit URF to last equation it is not really converging for you. What you should do is to use urf = 1 for this equation. Personally I would not use coupled solver for it, I would first solve cplus and cminus (may be coupled or segregated). Then I would solve the Laplacian with explicit urf. (Not implicit). Note that the situation here is quite similar to flow solver where momentum takes implicit urf but pressure takes explicit urf (because it is Poisson). In your situation even urf = 0.95 does not work (only way is to use urf = 1).. PS (edited to add): In case of Stokes flow where momentum is also Poisson then implicit urf is not applied and there we again apply explicit urf and it is difficult to solve. babakflame likes this.

June 7, 2017, 10:29		#31
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Arjun Thanks for your comment. So, basically, what you recommend is: Code: fvScalarMatrix UeEqn ( fvc::laplacian(eps,Ue) == -rhoE ); UeEqn.solve(); with no no urf (urf=1). I am gonna try it. But I think previously, just with Poisson (no cPlus and cMinus) explicit discretizing was giving me lower stability. I am gonna try it and give the feedback. Regards

May 27, 2017, 16:58		#22
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Thanks Philip The idea of coupled solver (coupling these 3 equations before pimple loop) is interesting for me too. I found out that Foam-Extend has some coupled solver examples. I will proceed through BC under-relaxation and coupled solver. Lets see how they perform. Regards

June 6, 2017, 20:35		#27
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Unfortunately, above snippet did not help. It seems that although it gaurantees that in each time step initial residual is less than a specific value, however, after couple of iterations, it reaches to a point that can not reduce initial residual and again divergence happens.

June 8, 2017, 21:42		#33
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Arjun, I realized that explicit discretizing of Poisson equation in my case leads to immediate divergence. Regards

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