[TemC]

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.

[piu58]

As far as I understand your code:

you have changed the original

volScalarField plugVisc

to

IOobject

. If that is really necessary, you have to cast your result.

[TemC]

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.

[Zeppo]

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

->

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

Foam::viscosityModels::myBingham::myBingham
(
    ...
)
:
    ...
    nu_
    (
        IOobject
        (
            name, //- "nu"?
            ...
        ),
        calcNu()
    )   
{}

[TemC]

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.

[Zeppo]

Yes, I meant

nu[j] = nu0_;

What is the error you get compiling this?

[TemC]

Here it is...

[Zeppo]

nu()[j] = nu0_;

[TemC]

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.

[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

               volScalarField 
               (
                   "plugVisc", dimensionSet(0,2,-1,0,0,0,0), nu0_.value()
               )

and not

            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

return nu1*pos(sr()[j] - criticalShearRate) + nu0

where nu1 is

k_ + tau0_*(1.0 - k_/nu0_)/(sr() + SMALL) - nu0

And so your calcNu should be just something like:

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

[TemC]

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.

[TemC]

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.

[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):

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

[TemC]

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.