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

Custom viscosity model

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

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   November 2, 2016, 08:14
Default Custom viscosity model
  #1
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Sponsored Links
Hi, I'm trying to implement a precise viscosity model from experimental measurements.

I've created a folder "Spline" and two files "Spline.C" and "Spline.H".
I've redefined calcNu() method in the following was:


Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    volScalarField strain = strainRate();
     volScalarField nu = strain;
     
//for each element of a field
     for (unsigned i=0; i<strain.size(); i++){
         double n = strain[i]; //get a strain-rate
    
        const double dn = 250/sizeof(visc);
        n /= dn; //scale it to match an index of array
    
/*visc array includes 2000 values of measured viscosity*/
        unsigned x1 = floor(n);
        if (x1 < 1) x1 = 0; //check minimum value
        if (x1 >= sizeof(visc)) x1 = sizeof(visc)-2; //check maximum value
        
        unsigned x2 = x1 + 1;
    
        if (x1 != x2){
            double y1 = visc[x1];
            double y2 = visc[x2];
        
            double k = y2 - y1; //linear interpolation
            double b = x2*y1 - x1*y2;
        
            nu[i] = k*n + b;
        }
        else{
            nu[i] = visc[x1];
        }
     }
    
    return nu;
}
It compiles with no error, however after start of pimpleFoam it exits with the following log:
Code:
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::viscosityModels::Spline::calcNu() const at ??:?
#4  Foam::viscosityModels::Spline::Spline(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#5  Foam::viscosityModel::adddictionaryConstructorToTable<Foam::viscosityModels::Spline>::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#6  Foam::viscosityModel::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#7  Foam::singlePhaseTransportModel::singlePhaseTransportModel(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#8  ? at ??:?
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? at ??:?
What can be wrong with my code ?
Svensen is offline   Reply With Quote
Sponsored Links

Old   November 2, 2016, 20:01
Default
  #2
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,978
Blog Entries: 39
Rep Power: 108
wyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of light
Greetings Svensen,

I've come to this thread, after looking for it due to the following bug report: http://bugs.openfoam.org/view.php?id=2316
In that bug report you have a slightly different source code:
Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
        volScalarField strain(strainRate());
     
        //for all elements in the field
     for (unsigned i=0; i<strain.size(); i++){
         double n = strain[i];
    
        const double dn = 250/sizeof(visc);
        n /= dn; //calc an index in array
    
        unsigned x1 = floor(n);
        if (x1 < 1) x1 = 0; //correct minimum value
        if (x1 >= sizeof(visc)) x1 = sizeof(visc)-2; //correct maximum value
        
        unsigned x2 = x1 + 1;
    
        if (x1 != x2){
            double y1 = visc[x1]; //visc array contains the measured viscosity values 
            double y2 = visc[x2];
        
            double k = y2 - y1; //because x2-x1 = 1
            double b = x2*y1 - x1*y2; //because x2-x1 = 1
        
            strain[i] = k*n + b;
        }
        else{
            strain[i] = visc[x1];
        }
     }
        
    return dimensionedScalar("one", dimLength*dimLength, 1.0)*strain;
}
Furthermore you commented that:
Quote:
After some experiments I've just found that problem is in line "strain[i] = k*n + b;". However I still don't understand why the problem happens here. OpenFOAM documentation states that volScalarFiels is based on List, therefore a value can be assigned to individual element. Is it some issue in OpenFOAM sources ?
It's very possible that you've isolated a symptom, but not the cause of the problem.
If we take a look at the printed stack trace, the first tell-tale sign is this:
Code:
#1  Foam::sigFpe::sigHandler(int) at ??:?
This means that there was a https://en.wikipedia.org/wiki/SIGFPE - which quoting from there:
Quote:
The SIGFPE signal is sent to a process when it executes an erroneous arithmetic operation, such as division by zero (the name "FPE", standing for floating-point exception, is a misnomer as the signal covers integer-arithmetic errors as well).
Other possible reasons is if there is a bad logarithm calculation or a division by infinite and so on.

When we take into account this signal for a bad mathematical operation and that the stack trace indicates that this occurs when the transport model is being created for the first time, I'm guessing that something wrong is going on with the initialization of the other fields:
  • It is not clear what "visc" is and how it was defined.
  • These two operations are the most likely suspects for a bad mathematical operation:
    Code:
             const double dn = 250/sizeof(visc);
            n /= dn; //calc an index in array
  • And as I mentioned in the bug report, I suspect that "sizeof(visc)" is simply being used incorrectly, because "sizeof" gives the number of bytes the variable needs on RAM, not the number of items on a list. Which means that if "visc" is in fact a list, then this means that any accesses with "visc[x1]" will result in illegal memory accesses.
When in doubt, use old-school debugging, i.e. use Info to output to screen the values that are being used, e.g.:
Code:
Info << "size visc: " << sizeof(visc) << endl;
Info << "dn: " << dn << endl;
Beyond this, the way that OpenFOAM goes through lists is by using forAll, e.g.:
Code:
forAll(strain, itemi)
{
   Info << strain[itemi] << endl;

}
Best regards,
Bruno
__________________
wyldckat is offline   Reply With Quote

Old   November 3, 2016, 08:12
Default
  #3
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Thanks, wyldcat.

The first error was with sizeof. I forgot that sizeof returns a number of bytes instead of number of elements.

However, the main problem was in the way of how "visc" array was initialized. I defined it as a private member of a class and initialized it in the body of constructor.

BUT OpenFOAM uses a constructor with initialization list, where nu_ object initializes before body of constructor is executed. That's why calcNu() method executes first and only after this my constructor will initialize a "visc" array. This produces "strange" values in "visc" array during execution of calcNu() method.

I've tried to move an initialization of nu_ object from initialization list to body of constructor, but compilation failed.
To overcome this, I've simply moved an initialization of "visc" array to the code of calcNu() method. It compiles without errors and executes first step without error, but then crashes.

I need some time to check my code to find out why it crashes.

But nevertheless, can I initialize a "visc" array outside from calcNu() method ? Because it seems weird to initialize it every time the calcNu() method executes.
I've attached the sources below.
Attached Files
File Type: txt Spline_C.txt (25.4 KB, 2 views)
File Type: txt Spline_H.txt (25.2 KB, 2 views)
Svensen is offline   Reply With Quote

Old   November 3, 2016, 10:04
Default
  #4
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
About crash...

First of all, I've started with Newtonian viscosity, i.e. I've set all elements of visc array to 4.1 mPa*s, which corresponds to kinematic viscosity of 3.9e-6 m2/s. I've checked that every element of "strain" object is initialized by this value during calcNu() method. However the simulation crashes again.

The return command looks like this:
Code:
return dimensionedScalar("one", pow(dimLength, 2), 1.0)*strain;
Error log is the following:
Code:
Courant Number mean: 0.970093 max: 303.708
Time = 0.002

PIMPLE: iteration 1
GAMG:  Solving for Ux, Initial residual = 0.931379, Final residual = 5.80411e-08, No Iterations 1000
GAMG:  Solving for Uy, Initial residual = 0.95176, Final residual = 2.56661e-08, No Iterations 1000
GAMG:  Solving for Uz, Initial residual = 0.924978, Final residual = 1.00825e-07, No Iterations 1000
GAMG:  Solving for p, Initial residual = 0.663798, Final residual = 8.97244e-08, No Iterations 48
time step continuity errors : sum local = 1.64355e-05, global = 1.05359e-05, cumulative = 1.05436e-05
GAMG:  Solving for p, Initial residual = 0.00216456, Final residual = 5.99289e-08, No Iterations 11
time step continuity errors : sum local = 0.00179684, global = -0.00101075, cumulative = -0.00100021
PIMPLE: iteration 2
[9] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[9] #1  Foam::sigFpe::sigHandler(int) at ??:?
[9] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
[9] #3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
[9] #4  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[9] #5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
[9] #6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
[9] #7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[9] #8  ? at ??:?
[9] #9  ? at ??:?
[9] #10  ? at ??:?
[9] #11  ? at ??:?
[9] #12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[9] #13  ? at ??:?
[sergey-HP:109750] *** Process received signal ***
[sergey-HP:109750] Signal: Floating point exception (8)
[sergey-HP:109750] Signal code:  (-6)
[sergey-HP:109750] Failing at address: 0x3e80001acb6
[sergey-HP:109750] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0)[0x7fc9651ffcb0]
[sergey-HP:109750] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37)[0x7fc9651ffc37]
[sergey-HP:109750] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0)[0x7fc9651ffcb0]
[sergey-HP:109750] [ 3] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam7sumProdIdEEdRKNS_5UListIT_EES5_+0x39)[0x7fc966566689]
[sergey-HP:109750] [ 4] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam5PBiCG5solveERNS_5FieldIdEERKS2_h+0x7d5)[0x7fc9663b93d5]
[sergey-HP:109750] [ 5] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver18solveCoarsestLevelERNS_5FieldIdEERKS2_+0x207)[0x7fc9663d46a7]
[sergey-HP:109750] [ 6] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver6VcycleERKNS_7PtrListINS_9lduMatrix8smootherEEERNS_5FieldIdEERKS8_S9_S9_S9_S9_S9_RNS1_IS8_EESD_h+0xbd8)[0x7fc9663d63f8]
[sergey-HP:109750] [ 7] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5solveERNS_5FieldIdEERKS2_h+0x49b)[0x7fc9663d818b]
[sergey-HP:109750] [ 8] pimpleFoam[0x459a28]
[sergey-HP:109750] [ 9] pimpleFoam[0x462b1b]
[sergey-HP:109750] [10] pimpleFoam[0x462dfb]
[sergey-HP:109750] [11] pimpleFoam[0x424978]
[sergey-HP:109750] [12] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7fc9651eaf45]
[sergey-HP:109750] [13] pimpleFoam[0x4253e8]
[sergey-HP:109750] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 9 with PID 109750 on node sergey-HP exited on signal 8 (Floating point exception).
It seems to me that something wrong happens when I return the "strain" object from calcNu() method. Maybe some pointer violation or type mismatch...
Svensen is offline   Reply With Quote

Old   November 4, 2016, 18:35
Default
  #5
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
...
    volScalarField strain(strainRate());
...
    return dimensionedScalar("one", dimLength*dimLength, 1.0)*strain;
}
->
Code:
Foam::viscosityModels::Spline::Spline
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    ...
    nu_(calcNu())
{}

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    const volScalarField& strainRate = strainRate()();
    tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strainRate.time().timeName(),
                strainRate.mesh(),
                strainRate.readOpt(),
                strainRate.writeOpt()
            ),
            strainRate.mesh(),
            sqr(dimLength)*strainRate.dimensions()
            strainRate.boundaryField().types()
        )
    );

    volScalarField& nu = tNu.ref();
    //volScalarField& nu = tNu(); //- for OF < 4.0
    forAll(nu, i)
    {
        ...
        nu[i] = ....
        ....
    }

    //- optionally you can also alter your boundary values:
    volScalarField::GeometricBoundaryField& bNu = nu.boundaryField()
    //volScalarField::Boundary& bNu = nu.boundaryFieldRef() //- for OF < 4.0
    forAll(bNu, i)
    {
        scalarField& pNu = bNu[i]; // patch field
        forAll(pNu, j)
        {
            pNu[j] = ... // patch cell face
        }
    }

    return tNu;
}
Zeppo is online now   Reply With Quote

Old   November 4, 2016, 19:53
Default
  #6
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Dear Zeppo,

thank you very much for help !

I've inserted your code in Spline.C and fixed some syntax issues in it. It compiles fine, however during execution the error with reference occurred. I've tried to comment every line of code from end to beginning and it seems that something wrong with this line:
Code:
tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strain_Rate.time().timeName(),
                strain_Rate.mesh(),
                strain_Rate.readOpt(),
                strain_Rate.writeOpt()
            ),
            strain_Rate.mesh(),
            sqr(dimLength)*strain_Rate.dimensions(),
            strain_Rate.boundaryField().types()
        )
    );
error log states that:
Code:
Selecting incompressible transport model Spline
[3] #0  Foam::error::printStack(Foam::Ostream&)[4] 
[4] 
[4] --> FOAM FATAL ERROR: 
[4] hanging pointer at index 0 (size 6), cannot dereference
[4] 
[4]     From function const T& Foam::UPtrList<T>::operator[](Foam::label) const [with T = Foam::fvPatchField<double>; Foam::label = int]
[4]     in file /home/sergey/OpenFOAM-dev/OpenFOAM-dev/src/OpenFOAM/lnInclude/UPtrListI.H at line 107.
[4] 
FOAM parallel run aborting
I've attached modified sources as well as full log below
Attached Files
File Type: c Spline.C (26.1 KB, 2 views)
File Type: h Spline.H (25.2 KB, 1 views)
File Type: txt log.txt (35.1 KB, 1 views)
Svensen is offline   Reply With Quote

Old   November 5, 2016, 06:11
Default
  #7
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
Ok. Then replace your calcNu() with the code following below:
Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    const volScalarField& strainRate = strainRate()();
    tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strainRate.time().timeName(),
                strainRate.mesh(),
                strainRate.readOpt(),
                strainRate.writeOpt()
            ),
            strainRate.mesh(),
            dimensionedScalar
            (   
                "oneNu", 
                sqr(dimLength)*strainRate.dimensions(), 
                1.0
            )
        )
    );
    return tNu;
}
Zeppo is online now   Reply With Quote

Old   November 5, 2016, 06:48
Default
  #8
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Dear Zeppo, thank you very much for help !

You suggestion is right and I was able to implement just a simple Newtonian model:
Code:
forAll(nu, i){
        nu[i] = visc[0];
    }
This code works fine when nu is not depend on strainRate. However, if I need a viscosity as a function of strain rate:
Code:
forAll(nu, i){
        nu[i] = function(strain_Rate[i]);
    }
it crashes.

It seems to me that i-index for nu is not the same as i-index for strain_Rate.

Maybe you know how the correct index for strain_Rate can be obtained ?
The sources as well as log file are attached below
Attached Files
File Type: txt log.txt (10.0 KB, 4 views)
File Type: txt Spline_C.txt (26.6 KB, 5 views)
File Type: txt Spline_H.txt (25.2 KB, 2 views)
Svensen is offline   Reply With Quote

Old   November 5, 2016, 07:17
Default
  #9
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
Both strain_Rate and nu have the same number of elements, so the code like this
Code:
forAll(nu, i)
{
    Info << nu[i] << " " << strain_Rate[i] <<endl;   
}
should be unconditionally valid. Commenting out the code you should isolate the problem to the single line which is the cause of the error
Zeppo is online now   Reply With Quote

Old   November 5, 2016, 07:55
Default
  #10
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Finally I've figured out a problem. I' found that this code works properly on a normal execution like "pimpleFoam", but code crashes on parallel run "pimpleFoam -parallel".

Maybe you know if there are some additional changes are required to source code to support parallel execution ?

Thanks !
Svensen is offline   Reply With Quote

Old   November 5, 2016, 10:40
Default
  #11
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
Which line triggers the error? This one?
Code:
double n = strain_Rate[i] / dn;
Zeppo is online now   Reply With Quote

Old   November 5, 2016, 11:43
Default
  #12
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Yes, you are right. Something wrong here.
Svensen is offline   Reply With Quote

Old   November 5, 2016, 13:11
Default
  #13
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
run this
Code:
forAll(nu, i)
{
    Info << i << ": " << nu[i] << " " << strain_Rate[i] <<endl;   
}
Zeppo is online now   Reply With Quote

Old   November 5, 2016, 13:55
Default
  #14
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
In the first step for every nu it prints i: 1 1 like
Code:
10380: 1 1
and so on.

At the beginning of second step it crashes before any output starts.
The full log is attached
Attached Files
File Type: gz log-parallel.txt.gz (51.8 KB, 2 views)
Svensen is offline   Reply With Quote

Old   November 5, 2016, 15:21
Default
  #15
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 192
Rep Power: 10
Zeppo is on a distinguished road
I am not an expert but what I see is that the error message never mentions your library file (.so-file). You built your custom viscosity class into a user library, right? It is libOpenFOAM.so that is mensioned in the error message as the sourse where the error comes from.
You said that the very same code runs flawlessly in serial, correct?
wyldckat likes this.
Zeppo is online now   Reply With Quote

Old   November 5, 2016, 15:26
Default
  #16
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Yes, serial execution works fine.
I will post this issue on issue tracker.

Thanks a lot for help !
Svensen is offline   Reply With Quote

Old   November 5, 2016, 15:49
Default
  #17
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,978
Blog Entries: 39
Rep Power: 108
wyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of light
Quick answer:
Quote:
Originally Posted by Svensen View Post
Yes, serial execution works fine.
I will post this issue on issue tracker.
Not a bug in OpenFOAM, it's a limitation of the implementation you've done

I need to understand what is the "visc" vector meant to be. Is it a lookup table?
wyldckat is offline   Reply With Quote

Old   November 5, 2016, 16:17
Default
  #18
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Yes, it is just a viscosity data. visc[0] contains the kinematic viscosity for shearRate=0; visc[1] contains the kinematic viscosity for shearRate=0.125; visc[2] contains viscosity for shearRate=0.250 and so on.

for shearRate between points, for example for shearRate=0.2 the linear interpolation will be used to compute nu for this point.

The sources are attached above.
Svensen is offline   Reply With Quote

Old   November 6, 2016, 04:42
Default
  #19
Senior Member
 
Join Date: Jan 2015
Posts: 142
Rep Power: 4
Svensen is on a distinguished road
Can this problem be treated ?
Svensen is offline   Reply With Quote

Old   November 6, 2016, 05:21
Default
  #20
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,978
Blog Entries: 39
Rep Power: 108
wyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of light
Greetings Svensen,

There are a few boundary conditions that demonstrate how lookup tables can be used. Using lookup tables as OpenFOAM already does, would ensure that the bug is not coming from a possibly flawed implementation of the lookup mechanism you currently have implemented.

The boundary "src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue" has such an example:
  1. The lookup table object is defined like this in the header file ".H":
    Code:
    autoPtr<Function1<Type>> uniformValue_;
    • In your case:
      Code:
      autoPtr<Function1<scalar>> visc;
  2. Constructed like this:
    Code:
    uniformValue_(Function1<Type>::New("uniformValue", dict))
    • In your case:
      Code:
      visc(Function1<scalar>::New("viscTable", dict))
      Then the entry "viscTable" needs to be in "transportProperties", within the respective "*Coefs" section.
  3. Then used like this:
    Code:
        const scalar t = this->db().time().timeOutputValue();
        fvPatchField<Type>::operator==(uniformValue_->value(t));
    • Although in your case, it would be something like:
      Code:
      visc->value(strainRate[celli])
  4. Don't forget to include the header file:
    Code:
    #include "Function1.H"
  5. In the "transportProperties" file, you will need something like what's described here: http://openfoam.org/release/2-1-0/bo...ime-dependent/ - possibly this:
    Code:
    viscTable     table
    (
        (0.0     0.014)
        (0.125   0.0123715)
        (0.25    0.011782)
    ...
        (????  0.00409976)
    );


If you still have trouble implementing it, please do the following steps, so it's easier to help you:
  1. Provide the source code by following these steps:
    1. Run inside the source code folder of your custom library:
      Code:
      wclean all
    2. Then go to the parent folder and package the source code folder, e.g. if the folder is named "splineTransport":
      Code:
      tar -czf splineLib.tar.gz splineTransport
    3. Then attach the file "splineLib.tar.gz" to your next post.
  2. Then we either need a test case, or instructions on how to change a tutorial in OpenFOAM to use your transport model.
Best regards,
Bruno
__________________
wyldckat is offline   Reply With Quote

Reply

Tags
openfoam-dev, viscosity model

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Time constant in Herschel-Bulkley viscosity model Mikel6 Main CFD Forum 0 October 17, 2016 04:52
New Viscosity Model - Best Practice for Temporary Values MaSch OpenFOAM Programming & Development 4 August 25, 2016 09:18
Problem with divergence TDK FLUENT 11 July 31, 2016 06:03
Viscosity ratio in gamma-theta transition model based on k-w sst turb model Qiaol618 Main CFD Forum 8 June 9, 2012 06:43
How to modify the viscosity model mpml OpenFOAM Running, Solving & CFD 4 October 13, 2010 07:44

Sponsored Links


All times are GMT -4. The time now is 19:07.