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

"Double slope" Bingham viscosity model

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

Like Tree2Likes
  • 1 Post By alexeym
  • 1 Post By TemC

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 1, 2017, 13:08
Default "Double slope" Bingham viscosity model
  #1
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Hi foamers,

I'am working with simpleFoam, and I'am trying to implement a "Double slope" Bingham viscosity model.
The 'classical' Bingham model has only two parameters: the yield stress "tau0", and the bingham viscosity "k".
The "Double slope" Bingham model suggest to add a new parameter "nu0", which is the viscosity of the central plugged zone.

With :
"muApp" = the apparent viscoity, and "sr()" = the shearRate ;

This model states:

if ( sr() >= criticalShearRate)
muApp = k_ + tau0_*(1 - k_/nu0_)/sr()
else
muApp = nu0_

Knowing that the parameter 'nu0_' is specified in constant/transportProperties, what I want is to create a volScalarField "plugVisc", which has uniform values of 'nu0_'. So that I can write the previous condition as :

forAll( sr(), j)
{
if ( sr()[j] >= criticalShearRate)
return
k_ + tau0_*(1 - k_/nu0_)/sr()
else
return
plugVisc()
}

The attachment "Bingham.C" shows how I'am currently trying to do the task. (I haven't change anything in the file "Bingham.H" ). When I try to compile with 'wmake libso', I got the message in attachment.

Does somebody have some advice about how I can create this volScalarField 'plugVisc' so that the compliation works?

Thanks in advance, and have a nice week.

Regards.
Attached Images
File Type: png Compilation_Message.PNG (121.1 KB, 15 views)
Attached Files
File Type: c Bingham.C (4.4 KB, 11 views)
File Type: h Bingham.H (3.3 KB, 8 views)
TemC is offline   Reply With Quote

Old   April 1, 2017, 14:15
Default
  #2
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 588
Rep Power: 10
piu58 is on a distinguished road
As far as I understand your code:

you have changed the original
Code:
volScalarField plugVisc
to
Code:
IOobject
. If that is really necessary, you have to cast your result.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   April 2, 2017, 02:38
Default
  #3
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Hi, and thanks for your reply.

Actually (for a Bingham fluid), the orgininal formulation of the model just states :

muApp = min (nu0_, k_ + tau0_/strainRate)

As you can see, it doesn't consider any condition on the strain rate. I have used this original formulation extensively, and It doesn't give me the expected theoretical results.

This is why now I want to implement the most commun formulation which takes into account a condition on the srain rate, as described in the first post.

And to do that, I need to create a uniform volScalarField with the value of the parameter 'nu0_'. This field I called it 'plugVisc', but it is just an arbitrary name and of course it can be changed.

My question is how do I create such a field, so that I can use it in my loop :

forAll( sr(), j)
{
if ( sr()[j] >= criticalShearRate)
return
k_ + tau0_*(1 - k_/nu0_)/sr()
else
return
plugVisc()
}

For the moment I'm trying to do that just by modifying the file Bingham.C. And apparently what I do doesn't work. I'am wondering how can I change my current declaration of this volScalarField 'plugVisc'? Do I have to declare it also into the file 'Bingham.H'? If yes, how?

At the end of the field Bingham.C, in comment, it is another way of declaring this volScalarField. I also tried this declaration, but the compilation also led to an error message.

Hope that more people will have a glance at it and share their point of view... Thanks in advance.

Regards.
TemC is offline   Reply With Quote

Old   April 2, 2017, 08:38
Default
  #4
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18
Zeppo will become famous soon enough
Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::myBingham::calcNu() const
{
    dimensionedScalar tone("tone", dimTime, 1.0);
    dimensionedScalar rtone("rtone", dimless/dimTime, 1.0);
    
    tmp<volScalarField> sr(strainRate());
    
    
    scalar criticalShearRate(tau0_.value()/nu0_.value());
    
    forAll(sr(), j)
    {      
        if ( sr()[j] >= criticalShearRate ) 
        {
            return
            (
               k_ + tau0_*(1.0 - k_/nu0_)/max(sr(), dimensionedScalar("VSMALL", dimless/dimTime, VSMALL))
            );  
        }
        else 
        {
            return
            (
               volScalarField 
               (     
                   plugVisc("plugVisc", dimensionSet(0,2,-1,0,0,0,0), nu0_.value()))
               );
            );
        } 
    } 
}
->
Code:
Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::myBingham::calcNu() const
{
    tmp<volScalarField> srtmp(strainRate());
    const volScalarField& sr = srtmp();

    scalar criticalShearRate(tau0_.value()/nu0_.value());
    
    tmp<volScalarField> nu
    (
        //- You should check the dimensions by yourself.
        k_
      + tau0_*(1.0 - k_/nu0_)
      / max
        (
            sr,
            dimensionedScalar("VSMALL", sr.dimensions(), VSMALL)
        )
    );
	
    forAll(sr, j)
    {      
        if ( sr[j] < criticalShearRate ) 
	{
		nu = nu0_;
	}
    }

    nu.correctBoundaryConditions();

    return nu;

}
Btw you might accidentally replaced "nu" with name
Code:
Foam::viscosityModels::myBingham::myBingham
(
    ...
)
:
    ...
    nu_
    (
        IOobject
        (
            name, //- "nu"?
            ...
        ),
        calcNu()
    )   
{}
Zeppo is offline   Reply With Quote

Old   April 2, 2017, 12:01
Default
  #5
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Thanks for your reply.

I understand your approach. Just something that puzzles me:

Into the loop forAll(sr, j), 'nu' is a volScalarField and 'nu0_' is a parameter. Of course they share the same unit, but I think we can't write "nu = nu0_"...

Actually I tried it and I got an error mesage while compiling...

I also tried something like : ' nu[j] = nu0.value() '

It also didn't work and now I realize that it is not appropriate, since the loop is done on the values of the field 'sr', and not 'nu'.

Any idea about how to adjust the line "nu = nu0_" into the loop? ...


Concerning the last point of your post, it was originally 'name', and I didn't change this part of the code.

Regards.
TemC is offline   Reply With Quote

Old   April 2, 2017, 12:27
Default
  #6
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18
Zeppo will become famous soon enough
Yes, I meant
Code:
nu[j] = nu0_;
What is the error you get compiling this?
Zeppo is offline   Reply With Quote

Old   April 2, 2017, 13:16
Default
  #7
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Here it is...
Attached Images
File Type: jpg Comp_Message.jpg (83.9 KB, 11 views)
TemC is offline   Reply With Quote

Old   April 2, 2017, 17:12
Default
  #8
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18
Zeppo will become famous soon enough
Code:
nu()[j] = nu0_;
Zeppo is offline   Reply With Quote

Old   April 3, 2017, 03:11
Default
  #9
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Morning,

Thanks for your reply. I tried what you suggested, now I just got what seems to be an error linked to the type of the scalars. Any idea how to fix it?..
Thanks again for your time.

Regards.
Attached Images
File Type: jpg Compilation.jpg (124.2 KB, 8 views)
TemC is offline   Reply With Quote

Old   April 3, 2017, 03:33
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,926
Rep Power: 35
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi all,

@TemC

1. The original error was caused by extra ) on line 76, a little bit strange usage of contructor, and extra ; here and there, i.e. in general you write

Code:
               volScalarField 
               (
                   "plugVisc", dimensionSet(0,2,-1,0,0,0,0), nu0_.value()
               )
and not

Code:
            return
            (
               volScalarField 
               (     
                   plugVisc("plugVisc", dimensionSet(0,2,-1,0,0,0,0), nu0_.value())
               );
            );
2. Since you are dealing with scalars, there is pos function (https://cpp.openfoam.org/v4/a09265_source.html#l00130), so your code has basically to do

Code:
return nu1*pos(sr()[j] - criticalShearRate) + nu0
where nu1 is

Code:
k_ + tau0_*(1.0 - k_/nu0_)/(sr() + SMALL) - nu0
And so your calcNu should be just something like:

Code:
    tmp<volScalarField> sr(strainRate());
    scalar sr_c(tau0_.value()/nu0_.value());
    volScalarField nu1(k_ + tau0_*(1.0 - k_/nu0_)/max(sr(), SMALL)) - nu0_;
    return nu1*pos(sr - sr_c) + nu0_;
tooran likes this.
alexeym is offline   Reply With Quote

Old   April 3, 2017, 04:54
Default
  #11
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Sir,

Thank you very much for sharing your knowledge with us. That's a sharp and very effective way to do the task.
After making a few adjustments to have the right units, like

return nu1*pos(sr - tone*sr_c) + nu0_; // where 'tone' is a dimensionedScalar with thhe dimension [1/s]

the compilation works fine now.

Again thank you all for your time and contributions. I really appreciate it.

Have a nive week.

Regards.
tooran likes this.
TemC is offline   Reply With Quote

Old   April 17, 2017, 06:28
Default
  #12
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Hi foamers,

I'am currently trying to make the viscosity model previously discussed more complex.
Up to this day, I was working with Bingham fluids. Which means that the parameter 'n' was always fixed to 1. Therefore the critical shear rate was:

sr_c = tau0.value()/nu0.value()

Now with the same approach, I want to handle the case of Herschel-Bulkley fluids. Meaning that 'n' can be 0.27 , 0.36 , 0.5 , 0.73, 1.18, ...

And this involves to solve an equation in order to find the critical shear rate 'sr_c'.

This equation is the following:

tau0 + k*(sr_c^n) = nu0*sr_c.

For the case n = 0.5, This equation is just transformed into an second order equation, and I solve it in the very classical way to have an analytical expression of 'sr_c'.

Now my challenge is when 'n' is different than 0.5 (let's say 0.36), I don't know how to deal with it...

Is there any way that I can to write a line of code like :

solve ( tau0 + k*(sr_c^n) == nu0*sr_c , sr_c) ?

So that right after this line, I can use the solution 'sr_c' in the expression of my viscosity model. Or should I write a code into a separate file to solve this equation?

If you have any idea or suggestions, please don't hesitate to share them. Of course any comment/suggestion will be much of help.

Thanks in advance, and have a good week.

Regards.
TemC is offline   Reply With Quote

Old   April 18, 2017, 02:17
Default
  #13
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,926
Rep Power: 35
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

If you target OpenFOAM, there are no root-finding methods.

In foam-extend there is findRoot (https://sourceforge.net/p/foam-exten.../ODE/findRoot/) and certain methods are implemented. Though it seems API provokes more coding than solving equations yourself.

If you really want to have proposed syntax, then you can write something like (you need to include header <functional> and to sort out things with units):

Code:
auto solve = [](std::function<scalar(scalar)> fn,
                std::function<scalar(scalar)> deriv,
                scalar x0,
                scalar tol = 1e-8,
                uLabel max_iter = 1024)
{
  uLabel iter = 1;
  scalar x = x0;
  for (;;) {
    scalar xn = x - fn(x)/deriv(x);
    if (std::abs(x - xn) <= tol or iter > max_iter) {
      return xn;
    }
    iter += 1;
    x = xn;
  }
};

sr_c = solve(
  [&tau0, &n, &nu0, &k](scalar x) -> scalar { return tau0 + k*pow(x, n) - nu0*x; },
  [&tau0, &n, &nu0, &k](scalar x) -> scalar { return k*n*pow(x, n - 1) - nu0; },
  tau0/nu0);
alexeym is offline   Reply With Quote

Old   April 18, 2017, 04:20
Default
  #14
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6
TemC is on a distinguished road
Sir,

Thank you very much for your reply, for the advice and for the code. I will take the time to digest it

Have a good week.

Regards.
TemC 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
Using a new implemented viscosity model with simpleFoam TemC OpenFOAM Running, Solving & CFD 6 March 8, 2017 03:07
Wrong multiphase flow at rotating interface Sanyo CFX 14 February 7, 2017 17:19
Time constant in Herschel-Bulkley viscosity model Mikel6 Main CFD Forum 0 October 17, 2016 04:52
modelling solids viscosity in eulerian multiphase model derkaiser FLUENT 1 December 5, 2011 03:42
Questions about Cross-Arrhenius and Cross-WLF viscosity model awacs OpenFOAM Running, Solving & CFD 4 August 13, 2009 06:56


All times are GMT -4. The time now is 00:26.