
[Sponsors] 
April 1, 2017, 13:08 
"Double slope" Bingham viscosity model

#1 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

April 1, 2017, 14:15 

#2 
Senior Member
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 588
Rep Power: 10 
As far as I understand your code:
you have changed the original Code:
volScalarField plugVisc Code:
IOobject
__________________
Uwe Pilz  Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950) 

April 2, 2017, 02:38 

#3 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

April 2, 2017, 08:38 

#4 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18 
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; } Code:
Foam::viscosityModels::myBingham::myBingham ( ... ) : ... nu_ ( IOobject ( name, // "nu"? ... ), calcNu() ) {} 

April 2, 2017, 12:01 

#5 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

April 2, 2017, 12:27 

#6 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18 
Yes, I meant
Code:
nu[j] = nu0_; 

April 2, 2017, 13:16 

#7 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
Here it is...


April 2, 2017, 17:12 

#8 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 255
Rep Power: 18 
Code:
nu()[j] = nu0_; 

April 3, 2017, 03:11 

#9 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

April 3, 2017, 03:33 

#10 
Senior Member

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() ) Code:
return ( volScalarField ( plugVisc("plugVisc", dimensionSet(0,2,1,0,0,0,0), nu0_.value()) ); ); Code:
return nu1*pos(sr()[j]  criticalShearRate) + nu0 Code:
k_ + tau0_*(1.0  k_/nu0_)/(sr() + SMALL)  nu0 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_; 

April 3, 2017, 04:54 

#11 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

April 17, 2017, 06:28 

#12 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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 HerschelBulkley 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. 

April 18, 2017, 02:17 

#13 
Senior Member

Hi,
If you target OpenFOAM, there are no rootfinding methods. In foamextend there is findRoot (https://sourceforge.net/p/foamexten.../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 = 1e8, 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); 

April 18, 2017, 04:20 

#14 
Member
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 6 
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. 

Thread Tools  Search this Thread 
Display Modes  


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 HerschelBulkley 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 CrossArrhenius and CrossWLF viscosity model  awacs  OpenFOAM Running, Solving & CFD  4  August 13, 2009 06:56 