[babakflame]

Dear Fellows

I am solving a 1D case with two transport equations for ions and a Poisson equation for additional parameter relating concentrations to each other. The simulation is stable and my initial residual keeps reducing.

However, an artificial oscillations occurs with my concentrations. I mean although the values seem fine , they alternatively flip sign. This is a sample of my cPlus concentration file written during simulation (first lines of concentration file)

internalField   nonuniform List<scalar> 
20000
(
-42240801498
-42240801498
-42240801497.7
-42240801497.9
-42240801497.9
-42240801497.9
-42240801498.1
-42240801498.1
-42240801498.1
-42240801498.1
6571685.77352
6571685.77353
6571685.77345
6571685.77353
6571685.7735
6571685.77351
6571685.77359
6571685.7736
6571685.77359
6571685.77358
-676.881846854
-676.881846855
-676.881846844
-676.881846857
-676.881846851
-676.881846853
-676.881846865
-676.881846866
-676.881846864
-676.881846863

I want to underline that values are sensible for me. However, not the sign change. I think sth is missing in my schemes. This is my fvSchems file:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
    default          Gauss linear;
    div(rho*phi,U)   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;
    cMinus;
    cPlus;
}

Any hint fellows? My concentration variables are cPlus and cMinus. There is no velocity in this simulation. Just Diffusion and time derivative plus electro-migrative terms.

So, I have to come up with other schemes for gradient and interpolation schemes.
I tried to employ NVD schemes for interpolationSchemes. I mean changing fvSchemes to:

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
    default          Gauss linear;
    div(rho*phi,U)   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 limited 1.0;
}

interpolationSchemes
{
    default         SFCD;
}

snGradSchemes
{
    default         limited 0.5;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
    cMinus;
    cPlus;
}

--> FOAM FATAL IO ERROR: 
attempt to read beyond EOF

file: /home/babak/Desktop/NoFlow/1D/Steric/a_3nm/Ue_500/Ue_100/system/fvSchemes::interpolationSchemes::default at line 45.

    From function ITstream::read(token&)
    in file db/IOstreams/Tstreams/ITstream.C at line 83.

FOAM exiting

Can anybody hint me how can I employ NVD schemes correctly or how can I remove this artificial oscillation?

Regards

[alexeym]

Hi,

(N|T)VD schemes need flux for limiters. In case of 'div(phi,?)' schemes can find flux, in other cases, they can not. So to use SFCD for interpolation (have doubts it helps with oscillations) you have to put 'SFCD phi'. Error message is awful but it is a result of "universality".

[babakflame]

Dear Alexey

Many thanks for your reply. Unfortunately, I am solving the equations without convective term.

My equations are as follows:

 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);



           surfaceScalarField cPlusArtDiff  = 
        -acubkappa*fvc::interpolate(cPlus)*mesh.magSf()*fvc::snGrad(cAgg)/(1-acub*fvc::interpolate(cPlus)-acub*fvc::interpolate(cMinus));
const double convergenceTolerance = 1e-02;
double initialResidualA = 0.0;
double initialResidualB = 0.0;
double initialResidualC = 0.0;
double initialResidual = 0.0;


do {
    cMinus.storePrevIter();
    cPlus.storePrevIter();
    Ue.storePrevIter();

       surfaceScalarField cMinusArtDiff  = 
        -acubkappa*fvc::interpolate(cMinus)*mesh.magSf()*fvc::snGrad(cAgg)/(1-acub*fvc::interpolate(cMinus)-acub*fvc::interpolate(cPlus));
           
           
            fvScalarMatrix cPlusEqn
            (
                fvm::ddt(cPlus)
              + fvc::div(cPlusFlux)
          - fvc::div(cPlusArtDiff)
           - fvm::laplacian(kappa, cPlus)    
            );
            initialResidualA = cPlusEqn.solve().initialResidual();


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

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

       
        cAgg = (cPlus + cMinus);

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


       initialResidualC = UeEqn.solve().initialResidual();
          initialResidual = (initialResidualC > initialResidual) ? initialResidualC : initialResidual;

} while (initialResidual > convergenceTolerance);

//          Convective = fvm::div(phi, rhoE);
    
           E = -fvc::grad(Ue);

    volVectorField gradUe = fvc::grad(Ue);

        gradUex = gradUe.component(vector::X);

        gradUey = gradUe.component(vector::Y);

As you can see, there is no convective term (phi in OpenFOAM notation). I tried NVD schemes through the following updated fvSchemes file:

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
  
}

laplacianSchemes
{
    default         Gauss linear limited 1.0;
}

interpolationSchemes
{
    default         linear;
cPlusFlux    Gamma 0.5;
cMinusFlux    Gamma 0.5;
cMinusArtDiff    Gamma 0.5;
cPlusArtDiff    Gamma 0.5;


}

snGradSchemes
{
    default         limited 0.5;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
    cMinus;
    cPlus;
}

However, I still face the problem of artificial oscillations in concentration values. I mean flipping signs.

Regards