generalised_newtonian_constitutive_models.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 
32 #ifndef OOMPH_GENERALISED_NEWTONIAN_CONSTITUTIVE_MODELS_HEADER
33 #define OOMPH_GENERALISED_NEWTONIAN_CONSTITUTIVE_MODELS_HEADER
34 
35 
36 //Oomph-lib includes
37 //#include "../generic.h"
38 
39 namespace oomph
40 {
41 
42 
43 //======================================================================
44 /// A Base class defining the generalise Newtonian constitutive relation
45 //======================================================================
46  template<unsigned DIM>
48  {
49 
50  public:
51 
52  /// Empty constructor
54 
55 
56  /// Empty virtual destructor
58 
59  /// Function implementing the constitutive model
60  /// Input: second invariant of the rate of strain
61  /// Output: the viscosity
62  /// For Newtonian behaviour this returns 1
63  virtual double viscosity(const double&
64  second_invariant_of_rate_of_strain_tensor)=0;
65 
66  /// Function returning the derivative of the viscosity w.r.t.
67  /// the second invariant of the rate of strain tensor
68  /// For Newtonian behaviour this returns 0.0
69  virtual double dviscosity_dinvariant
70  (const double& second_invariant_of_rate_of_strain_tensor)=0;
71 
72  };
73 
74 //===================================================================
75 /// A GeneralisedNewtonianConstitutiveEquation class
76 /// defining a Newtonian fluid
77 //===================================================================
78  template<unsigned DIM>
81  {
82 
83  public:
84 
85  /// Constructor: specify viscosity ratio (defaults to one)
86  NewtonianConstitutiveEquation(const double& viscosity_ratio=1.0) :
87  Viscosity_ratio(viscosity_ratio) {}
88 
89  /// in the Newtonian case the viscosity is constant
90  double viscosity(const double&
91  second_invariant_of_rate_of_strain_tensor)
92  {
93  return Viscosity_ratio;
94  }
95 
96  /// the derivative w.r.t. I2 is zero
98  (const double& second_invariant_of_rate_of_strain_tensor)
99  {
100  return 0.0;
101  }
102 
103  private:
104 
105  /// Viscosity ratio
107 
108  };
109 
110 
111 
112 //===================================================================
113 /// A GeneralisedNewtonianConstitutiveEquation class defining a power-law fluid
114 /// regularised according to Bercovier and Engelman (1980) to allow for n < 1
115 //==================================================================
116 template<unsigned DIM>
119 {
120  private:
121 
122  /// power law index n
123  double* Power_pt;
124 
125  /// regularisation parameter e << 1
127 
128 
129  public:
130 
131  PowerLawBerEngRegConstitutiveEquation(double* power_pt, double* reg_par_pt) :
133  Power_pt(power_pt), Regularisation_parameter_pt(reg_par_pt){}
134 
135  double viscosity(const double&
136  second_invariant_of_rate_of_strain_tensor)
137  {
138  // Pre-multiply the second invariant with +/-1 depending on whether it's
139  // positive or not
140  double sign=-1.0;
141  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
142  {
143  sign=1.0;
144  }
145 
146  // Calculate the square root of the absolute value of the
147  // second invariant of the rate of strain tensor
148  double measure_of_rate_of_strain=
149  sqrt(sign*second_invariant_of_rate_of_strain_tensor);
150 
151  return pow((2.0*measure_of_rate_of_strain + *Regularisation_parameter_pt),
152  *Power_pt-1.0);
153  }
154 
155 
156 };
157 
158 //===================================================================
159 /// A GeneralisedNewtonianConstitutiveEquation class
160 /// defining a Herschel-Bulkley fluid
161 /// using Bercovier and Engelman's (1980) regularisation
162 //==================================================================
163 template<unsigned DIM>
166 {
167  private:
168 
169  /// yield stress tau_y
171 
172  /// power law index n
173  double* Flow_index_pt;
174 
175  /// regularisation parameter e << 1
177 
178 
179  public:
180 
182  (double* yield_stress_pt, double* flow_index_pt,
183  double* regularisation_parameter_pt) :
185  Yield_stress_pt(yield_stress_pt), Flow_index_pt(flow_index_pt),
186  Regularisation_parameter_pt(regularisation_parameter_pt){}
187 
188  double viscosity(const double&
189  second_invariant_of_rate_of_strain_tensor)
190  {
191  // Pre-multiply the second invariant with +/-1 depending on whether it's
192  // positive or not
193  double sign=-1.0;
194  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
195  {
196  sign=1.0;
197  }
198 
199  // Calculate the square root of the absolute value of the
200  // second invariant of the rate of strain tensor
201  double measure_of_rate_of_strain=
202  sqrt(sign*second_invariant_of_rate_of_strain_tensor+
204 
205  return (*Yield_stress_pt)/(2.0*measure_of_rate_of_strain)+
206  pow(2.0*measure_of_rate_of_strain,*Flow_index_pt-1.0);
207  }
208 
209  // not implemented yet
210  double dviscosity_dinvariant
211  (const double& second_invariant_of_rate_of_strain_tensor)
212  {
213  // Pre-multiply the second invariant with +/-1 depending on whether it's
214  // positive or not
215  double sign=-1.0;
216  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
217  {
218  sign=1.0;
219  }
220 
221  // std::ostringstream error_stream;
222  // error_stream << "This has not been implemented yet!";
223  // throw OomphLibError(
224  // error_stream.str(),
225  // OOMPH_CURRENT_FUNCTION,
226  // OOMPH_EXCEPTION_LOCATION);
227 
228  return sign*pow(2.0,(*Flow_index_pt)-2.0)*((*Flow_index_pt)-1.0)*
229  pow(sign*second_invariant_of_rate_of_strain_tensor+
231  ((*Flow_index_pt)-1.0)/2.0-1.0)-sign*(*Yield_stress_pt)/
232  (4.0*pow(sign*second_invariant_of_rate_of_strain_tensor+
233  (*Regularisation_parameter_pt),3.0/2.0));
234  }
235 
236 };
237 
238 //===================================================================
239 /// A GeneralisedNewtonianConstitutiveEquation class
240 /// defining a Herschel-Bulkley fluid
241 /// using Tanner and Milthorpe's (1983) regularisation
242 //==================================================================
243 template<unsigned DIM>
246 {
247  private:
248 
249  /// yield stress tau_y
251 
252  /// power law index n
253  double* Flow_index_pt;
254 
255  /// value of the second invariant below which we have
256  /// constant (Newtonian) viscosity -- assumed to be always positive
258 
259  public:
260 
261 
262  /// "Cutoff regularised" Herschel Bulkley constitutive equation
264  (double* yield_stress_pt, double* flow_index_pt,
265  double* critical_second_invariant_pt) :
267  Yield_stress_pt(yield_stress_pt), Flow_index_pt(flow_index_pt),
268  Critical_second_invariant_pt(critical_second_invariant_pt)
269  {
270 
271  // Calculate the Newtonian cutoff viscosity from the constitutive
272  // equation and the cutoff value of the second invariant
273  double cut_off_viscosity=calculate_cut_off_viscosity();
274 
275  oomph_info << "HerschelBulkleyTanMilRegConstitutiveEquation: "
276  << " cutoff viscosity = " << cut_off_viscosity
277  << std::endl;
278 
279  oomph_info << " "
280  << " cutoff invariant = "
282  << std::endl;
283  }
284 
285  /// Function that calculates the cut off viscosity
287  {
288  return (*Yield_stress_pt)/
289  (2.0*sqrt((*Critical_second_invariant_pt)))+
290  pow((2.0*sqrt((*Critical_second_invariant_pt))),
291  *Flow_index_pt-1.0);
292  }
293 
294  /// Report cutoff values
295  void report_cut_off_values(double& cut_off_invariant,
296  double& cut_off_viscosity)
297  {
298  cut_off_invariant=*Critical_second_invariant_pt;
299  cut_off_viscosity=calculate_cut_off_viscosity();
300  }
301 
302 
303  /// Viscosity ratio as a fct of strain rate invariant
304  double viscosity(const double&
305  second_invariant_of_rate_of_strain_tensor)
306  {
307  // Pre-multiply the second invariant with +/-1 depending on whether it's
308  // positive or not
309  double sign=-1.0;
310  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
311  {
312  sign=1.0;
313  }
314 
315  //std::cout<<"I2 "<<second_invariant_of_rate_of_strain_tensor<<std::endl;
316 
317  // if the second invariant is below the cutoff we have a constant, Newtonian
318  // viscosity
319  if(sign*second_invariant_of_rate_of_strain_tensor <
321  {
323  }
324  else
325  {
326  // Calculate the square root of the absolute value of the
327  // second invariant of the rate of strain tensor
328  double measure_of_rate_of_strain=
329  sqrt(sign*second_invariant_of_rate_of_strain_tensor);
330 
331  return (*Yield_stress_pt)/(2.0*measure_of_rate_of_strain)+
332  pow((2.0*measure_of_rate_of_strain),*Flow_index_pt-1.0);
333  }
334  }
335 
336  /// Deriv of viscosity w.r.t. strain rate invariant
337  double dviscosity_dinvariant
338  (const double& second_invariant_of_rate_of_strain_tensor)
339  {
340  // Pre-multiply the second invariant with +/-1 depending on whether it's
341  // positive or not
342  double sign=-1.0;
343  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
344  {
345  sign=1.0;
346  }
347 
348  if(sign*second_invariant_of_rate_of_strain_tensor <
350  {
351  return 0.0;
352  }
353  else
354  {
355  return sign*pow(2.0,(*Flow_index_pt)-2.0)*((*Flow_index_pt)-1.0)*
356  pow(sign*second_invariant_of_rate_of_strain_tensor,
357  ((*Flow_index_pt)-1.0)/2.0-1.0)-sign*(*Yield_stress_pt)/
358  (4.0*pow(sign*second_invariant_of_rate_of_strain_tensor,3.0/2.0));
359  }
360  }
361 
362 };
363 
364 //===================================================================
365 /// A GeneralisedNewtonianConstitutiveEquation class
366 /// defining a Herschel-Bulkley fluid
367 /// using Tanner and Milthorpe's (1983) regularisation
368 /// with a smooth transition using a cubic
369 //==================================================================
370 template<unsigned DIM>
373 {
374  private:
375 
376  /// yield stress tau_y
378 
379  /// power law index n
380  double* Flow_index_pt;
381 
382  /// value of the second invariant below which we have
383  /// constant (Newtonian) viscosity -- assumed to be always positive
385 
386 
387  public:
388 
389 
390  /// "Cutoff regularised" Herschel Bulkley constitutive equation
392  (double* yield_stress_pt, double* flow_index_pt,
393  double* critical_second_invariant_pt) :
395  Yield_stress_pt(yield_stress_pt), Flow_index_pt(flow_index_pt),
396  Critical_second_invariant_pt(critical_second_invariant_pt)
397  {
398 
399  /// get the cutoff viscosity
400  double cut_off_viscosity=calculate_cutoff_viscosity();
401 
402  /// get the zero shear viscosity
403  double zero_shear_viscosity=calculate_zero_shear_viscosity();
404 
405  oomph_info << "HerschelBulkleyTanMilRegWithBlendingConstitutiveEquation: "
406  << " zero shear viscosity = " << zero_shear_viscosity
407  << std::endl;
408 
409  oomph_info << "HerschelBulkleyTanMilRegWithBlendingConstitutiveEquation: "
410  << " cut off viscosity = " << cut_off_viscosity
411  << std::endl;
412 
413  oomph_info << " "
414  << " cutoff invariant = "
416  << std::endl;
417 
418 
419  }
420 
421  /// Function that calculates the viscosity at the cutoff invariant
422  /// Note: this is NOT the viscosity at zero I2
424  {
425  // Calculate the Newtonian cutoff viscosity from the constitutive
426  // equation and the cutoff value of the second invariant
427  return (*Yield_stress_pt)/
428  (2.0*sqrt((*Critical_second_invariant_pt)))+
429  pow((2.0*sqrt((*Critical_second_invariant_pt))),
430  *Flow_index_pt-1.0);
431  }
432 
433  /// Offset by how much the zero shear rate viscosity lies
434  /// above the viscosity at I2_cutoff
435  /// Hard-coded to a value that ensures a smooth transition
436  double calculate_viscosity_offset_at_zero_shear(double& cut_off_viscosity)
437  {
438  return cut_off_viscosity/5.0;
439  }
440 
441  /// Function that calculates the viscosity at zero I2
443  {
444  // get the viscosity at the cutoff point
445  double cut_off_viscosity=calculate_cutoff_viscosity();
446 
447  /// Offset by how much the zero shear rate viscosity lies
448  /// above the viscosity at I2_cutoff
449  double epsilon=calculate_viscosity_offset_at_zero_shear(cut_off_viscosity);
450 
451  return cut_off_viscosity+epsilon;
452  }
453 
454  /// Report cutoff values
455  void report_cut_off_values(double& cut_off_invariant,
456  double& cut_off_viscosity,
457  double& zero_shear_viscosity)
458  {
459  cut_off_invariant=*Critical_second_invariant_pt;
460  cut_off_viscosity=calculate_cutoff_viscosity();
461  zero_shear_viscosity=calculate_zero_shear_viscosity();
462  }
463 
464  // Calculate the fitting parameters for the cubic blending
465  void calculate_fitting_parameters_of_cubic(double& a, double& b)
466  {
467  // get the viscosity at the cutoff invariant
468  double Cut_off_viscosity=calculate_cutoff_viscosity();
469 
470  // calculate the offset at zero shear
471  double epsilon=calculate_viscosity_offset_at_zero_shear(Cut_off_viscosity);
472 
473  a=-1.0/pow((*Critical_second_invariant_pt),4.0)*
474  (pow((*Critical_second_invariant_pt),2.0)*
475  (-pow(2.0,((*Flow_index_pt)-2.0))*((*Flow_index_pt)-1.0)*
477  (((*Flow_index_pt)-1.0)/2.0-1.0))+
478  (*Yield_stress_pt)/(4.0*pow((*Critical_second_invariant_pt),
479  (3.0/2.0))))-
480  2.0*(*Critical_second_invariant_pt)*
481  (Cut_off_viscosity+epsilon-pow(2.0,((*Flow_index_pt)-1.0))*
482  pow((*Critical_second_invariant_pt),(((*Flow_index_pt)-1.0)/2.0))-
483  (*Yield_stress_pt)/(2.0*sqrt((*Critical_second_invariant_pt)))));
484 
485  b=-1.0/(8.0*pow((*Critical_second_invariant_pt),5.0)*
486  pow((*Critical_second_invariant_pt),(7.0/2.0)))*
487  (24.0*epsilon*pow((*Critical_second_invariant_pt),3.0)*
488  pow((*Critical_second_invariant_pt),(7.0/2.0))+
489  24.0*Cut_off_viscosity*pow((*Critical_second_invariant_pt),3.0)*
490  pow((*Critical_second_invariant_pt),(7.0/2.0))-
491  pow(2.0,((*Flow_index_pt)+1.0))*
492  pow((*Critical_second_invariant_pt),4.0)*
493  pow((*Critical_second_invariant_pt),(((*Flow_index_pt)+4.0)/2.0))+
494  pow(2.0,((*Flow_index_pt)+1.0))*(*Flow_index_pt)*
495  pow((*Critical_second_invariant_pt),4.0)*
496  pow((*Critical_second_invariant_pt),(((*Flow_index_pt)+4.0)/2.0))-
497  12.0*pow(2.0,(*Flow_index_pt))*pow((*Critical_second_invariant_pt),3.0)*
498  pow((*Critical_second_invariant_pt),(((*Flow_index_pt)+6.0)/2.0))-
499  2.0*(*Yield_stress_pt)*pow((*Critical_second_invariant_pt),4.0)*
500  pow((*Critical_second_invariant_pt),2.0)-12.0*(*Yield_stress_pt)*
501  pow((*Critical_second_invariant_pt),3.0)*
502  pow((*Critical_second_invariant_pt),3.0));
503  }
504 
505 
506  /// Viscosity ratio as a fct of strain rate invariant
507  double viscosity(const double&
508  second_invariant_of_rate_of_strain_tensor)
509  {
510  // Get the parameters of the cubic
511  double a;
512  double b;
513 
515 
516  double zero_shear_viscosity=calculate_zero_shear_viscosity();
517 
518  // Pre-multiply the second invariant with +/-1 depending on whether it's
519  // positive or not
520  double sign=-1.0;
521  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
522  {
523  sign=1.0;
524  }
525 
526  // if the second invariant is below the cutoff we have a constant, Newtonian
527  // viscosity
528  if(sign*second_invariant_of_rate_of_strain_tensor <
530  {
531  return a*pow(sign*second_invariant_of_rate_of_strain_tensor,3.0)+
532  b*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
533  zero_shear_viscosity;
534  }
535  else
536  {
537  // Calculate the square root of the absolute value of the
538  // second invariant of the rate of strain tensor
539  double measure_of_rate_of_strain=
540  sqrt(sign*second_invariant_of_rate_of_strain_tensor);
541 
542  return (*Yield_stress_pt)/(2.0*measure_of_rate_of_strain)+
543  pow((2.0*measure_of_rate_of_strain),*Flow_index_pt-1.0);
544  }
545  }
546 
547  /// Deriv of viscosity w.r.t. strain rate invariant
548  double dviscosity_dinvariant
549  (const double& second_invariant_of_rate_of_strain_tensor)
550  {
551  // Get the parameters of the cubic
552  double a;
553  double b;
554 
556 
557  // Pre-multiply the second invariant with +/-1 depending on whether it's
558  // positive or not
559  double sign=-1.0;
560  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
561  {
562  sign=1.0;
563  }
564 
565  if(sign*second_invariant_of_rate_of_strain_tensor <
567  {
568  return sign*3.0*a*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
569  2.0*b*second_invariant_of_rate_of_strain_tensor;
570  }
571  else
572  {
573  return pow(2.0,(*Flow_index_pt)-2.0)*((*Flow_index_pt)-1.0)*
574  second_invariant_of_rate_of_strain_tensor*
575  pow(sign*second_invariant_of_rate_of_strain_tensor,
576  ((*Flow_index_pt)-1.0)/2.0-2.0)-
577  (*Yield_stress_pt)*second_invariant_of_rate_of_strain_tensor/
578  (4.0*pow(sign*second_invariant_of_rate_of_strain_tensor,5.0/2.0));
579  }
580  }
581 
582 };
583 
584 
585 //===================================================================
586 /// A GeneralisedNewtonianConstitutiveEquation class
587 /// defining a Herschel-Bulkley fluid
588 /// using Papanastasiou's (1987) regularisation
589 //==================================================================
590 template<unsigned DIM>
593 {
594  private:
595 
596  /// Yield stress tau_y
598 
599  /// Power law index n
600  double* Flow_index_pt;
601 
602  /// Regularisation parameter m >> 1
604 
605 
606  public:
607 
609  (double* yield_stress_pt, double* flow_index_pt,
610  double* exponential_parameter_pt) :
612  Yield_stress_pt(yield_stress_pt), Flow_index_pt(flow_index_pt),
613  Exponential_parameter_pt(exponential_parameter_pt){}
614 
615  double viscosity(const double&
616  second_invariant_of_rate_of_strain_tensor)
617  {
618  // Calculate magnitude of the rate of strain tensor
619  double measure_of_rate_of_strain=
620  sqrt(std::fabs(second_invariant_of_rate_of_strain_tensor));
621 
622  return (1.0-exp(-2.0*(*Exponential_parameter_pt)*
623  measure_of_rate_of_strain))/(2.0*measure_of_rate_of_strain)*
624  (*Yield_stress_pt)+
625  (1.0-exp(-2.0*(*Exponential_parameter_pt)*
626  measure_of_rate_of_strain))/(2.0*measure_of_rate_of_strain)*
627  pow((2.0*measure_of_rate_of_strain),*Flow_index_pt);
628  }
629 
630 };
631 
632 //===================================================================
633 /// A GeneralisedNewtonianConstitutiveEquation class
634 /// defining a Herschel-Bulkley fluid
635 /// using Mendes and Dutra's (2004) regularisation
636 //==================================================================
637 template<unsigned DIM>
640 {
641  private:
642 
643  /// yield stress tau_y
645 
646  /// power law index n
647  double* Flow_index_pt;
648 
649  /// the viscosity at zero shear rate
651 
652  public:
653 
654 
655  /// "Exponentially regularised" Herschel Bulkley constitutive equation
657  (double* yield_stress_pt, double* flow_index_pt,
658  double* zero_shear_viscosity_pt) :
660  Yield_stress_pt(yield_stress_pt), Flow_index_pt(flow_index_pt),
661  Zero_shear_viscosity_pt(zero_shear_viscosity_pt){}
662 
663  /// Viscosity ratio as a fct of strain rate invariant
664  double viscosity(const double&
665  second_invariant_of_rate_of_strain_tensor)
666  {
667  // Pre-multiply the second invariant with +/-1 depending on whether it's
668  // positive or not
669  // Also, because the viscosity is exactly zero at zero invariant,
670  // we have to add a small value to it
671  double sign=-1.0;
672  double eps=0.0;
673  if(second_invariant_of_rate_of_strain_tensor == 0.0)
674  {
675  eps=1.0e-30;
676  sign=1.0;
677  }
678  else if(second_invariant_of_rate_of_strain_tensor > 0.0)
679  {
680  sign=1.0;
681  }
682 
683  // Calculate the square root of the absolute value of the
684  // second invariant of the rate of strain tensor
685  double measure_of_rate_of_strain=
686  sqrt(sign*(second_invariant_of_rate_of_strain_tensor+eps));
687 
688  return (1.0-exp(-(*Zero_shear_viscosity_pt)*
689  2.0*measure_of_rate_of_strain/(*Yield_stress_pt)))*
690  ((*Yield_stress_pt)/(2.0*measure_of_rate_of_strain)+
691  pow((2.0*measure_of_rate_of_strain),*Flow_index_pt-1.0));
692 
693  }
694 
695  /// Deriv of viscosity w.r.t. strain rate invariant
696  double dviscosity_dinvariant
697  (const double& second_invariant_of_rate_of_strain_tensor)
698  {
699  // Pre-multiply the second invariant with +/-1 depending on whether it's
700  // positive or not
701  // Also, because the viscosity is exactly zero at zero invariant,
702  // we have to add a small value to it
703  double sign=-1.0;
704  double eps=0.0;
705  if(second_invariant_of_rate_of_strain_tensor == 0.0)
706  {
707  eps=1.0e-30;
708  sign=1.0;
709  }
710  else if(second_invariant_of_rate_of_strain_tensor > 0.0)
711  {
712  sign=1.0;
713  }
714 
715  // Calculate the square root of the absolute value of the
716  // second invariant of the rate of strain tensor
717  double measure_of_rate_of_strain=
718  sqrt(sign*(second_invariant_of_rate_of_strain_tensor+eps));
719 
720  return (1.0-exp(-(*Zero_shear_viscosity_pt)*
721  2.0*measure_of_rate_of_strain/(*Yield_stress_pt)))*
722  (sign*pow(2.0,(*Flow_index_pt)-2.0)*((*Flow_index_pt)-1.0)*
723  pow(sign*(second_invariant_of_rate_of_strain_tensor+eps),
724  ((*Flow_index_pt)-1.0)/2.0-1.0)-sign*(*Yield_stress_pt)/
725  (4.0*pow(sign*(second_invariant_of_rate_of_strain_tensor+eps),3.0/2.0)))+
726  sign*(((*Zero_shear_viscosity_pt)*
727  exp(-(*Zero_shear_viscosity_pt)*
728  2.0*measure_of_rate_of_strain/(*Yield_stress_pt)))/
729  ((*Yield_stress_pt)*measure_of_rate_of_strain))*
730  (pow(2.0,(*Flow_index_pt)-1.0)*
731  pow(sign*(second_invariant_of_rate_of_strain_tensor+eps),
732  ((*Flow_index_pt)-1.0)/2.0)+
733  (*Yield_stress_pt)/(2.0*measure_of_rate_of_strain));
734 
735  }
736 
737 };
738 
739 //===================================================================
740 /// A GeneralisedNewtonianConstitutiveEquation class
741 /// defining a Sisko fluid
742 /// using Tanner and Milthorpe's (1983) regularisation
743 /// with a smooth transition using a cubic (for n < 1)
744 //==================================================================
745 template<unsigned DIM>
748 {
749  private:
750 
751  /// pre-factor alpha
752  double* Alpha_pt;
753 
754  /// power law index n
755  double* Flow_index_pt;
756 
757  /// value of the second invariant below which we have
758  /// constant (Newtonian) viscosity -- assumed to be always positive
760 
761 
762  public:
763 
764 
765  /// "Cutoff regularised" Sisko constitutive equation
767  (double* alpha_pt, double* flow_index_pt,
768  double* critical_second_invariant_pt) :
770  Alpha_pt(alpha_pt), Flow_index_pt(flow_index_pt),
771  Critical_second_invariant_pt(critical_second_invariant_pt)
772  {
773 
774  /// get the cutoff viscosity
775  double cut_off_viscosity=calculate_cutoff_viscosity();
776 
777  /// get the zero shear viscosity
778  double zero_shear_viscosity=calculate_zero_shear_viscosity();
779 
780  oomph_info << "SiskoTanMilRegWithBlendingConstitutiveEquation: "
781  << " zero shear viscosity = " << zero_shear_viscosity
782  << std::endl;
783 
784  oomph_info << "SiskoTanMilRegWithBlendingConstitutiveEquation: "
785  << " cut off viscosity = " << cut_off_viscosity
786  << std::endl;
787 
788  oomph_info << " "
789  << " cutoff invariant = "
791  << std::endl;
792 
793 
794  }
795 
796  /// Function that calculates the viscosity at the cutoff invariant
797  /// Note: this is NOT the viscosity at zero I2
799  {
800  // Calculate the Newtonian cutoff viscosity from the constitutive
801  // equation and the cutoff value of the second invariant
802  return 1.0+(*Alpha_pt)*pow((2.0*sqrt((*Critical_second_invariant_pt))),
803  *Flow_index_pt-1.0);
804  }
805 
806  /// Offset by how much the zero shear rate viscosity lies
807  /// above the viscosity at I2_cutoff
808  /// Hard-coded to a value that ensures a smooth transition
809  double calculate_viscosity_offset_at_zero_shear(double& cut_off_viscosity)
810  {
811  return cut_off_viscosity/5.0;
812  }
813 
814  /// Function that calculates the viscosity at zero I2
816  {
817  // get the viscosity at the cutoff point
818  double cut_off_viscosity=calculate_cutoff_viscosity();
819 
820  /// Offset by how much the zero shear rate viscosity lies
821  /// above the viscosity at I2_cutoff
822  double epsilon=calculate_viscosity_offset_at_zero_shear(cut_off_viscosity);
823 
824  return cut_off_viscosity+epsilon;
825  }
826 
827  /// Report cutoff values
828  void report_cut_off_values(double& cut_off_invariant,
829  double& cut_off_viscosity,
830  double& zero_shear_viscosity)
831  {
832  cut_off_invariant=*Critical_second_invariant_pt;
833  cut_off_viscosity=calculate_cutoff_viscosity();
834  zero_shear_viscosity=calculate_zero_shear_viscosity();
835  }
836 
837  // Calculate the fitting parameters for the cubic blending
838  void calculate_fitting_parameters_of_cubic(double& a, double& b)
839  {
840  // get the viscosity at the cutoff invariant
841  double Cut_off_viscosity=calculate_cutoff_viscosity();
842 
843  // calculate the offset at zero shear
844  double epsilon=calculate_viscosity_offset_at_zero_shear(Cut_off_viscosity);
845 
846  a=1.0/(4.0*pow((*Critical_second_invariant_pt),3.0)*
847  pow((*Critical_second_invariant_pt),5.0/2.0))*
848  (8.0*(Cut_off_viscosity+epsilon-1.0)*
849  pow((*Critical_second_invariant_pt),5.0/2.0)+
850  pow(2.0,(*Flow_index_pt))*((*Flow_index_pt)-1.0)*
851  (*Alpha_pt)*pow((*Critical_second_invariant_pt),2.0)*
853  pow(2.0,(*Flow_index_pt)+2.0)*(*Alpha_pt)*
854  pow((*Critical_second_invariant_pt),(*Flow_index_pt)/2.0+2.0));
855 
856  b=1.0/(4.0*pow((*Critical_second_invariant_pt),2.0)*
857  pow((*Critical_second_invariant_pt),5.0/2.0))*
858  (-12.0*(Cut_off_viscosity+epsilon-1.0)*
859  pow((*Critical_second_invariant_pt),5.0/2.0)-
860  pow(2.0,(*Flow_index_pt))*((*Flow_index_pt)-1.0)*
861  (*Alpha_pt)*pow((*Critical_second_invariant_pt),2.0)*
863  3.0*pow(2.0,(*Flow_index_pt)+1.0)*(*Alpha_pt)*
864  pow((*Critical_second_invariant_pt),(*Flow_index_pt)/2.0+2.0));
865  }
866 
867 
868  /// Viscosity ratio as a fct of strain rate invariant
869  double viscosity(const double&
870  second_invariant_of_rate_of_strain_tensor)
871  {
872  // Get the parameters of the cubic
873  double a;
874  double b;
875 
877 
878  double zero_shear_viscosity=calculate_zero_shear_viscosity();
879 
880  // Pre-multiply the second invariant with +/-1 depending on whether it's
881  // positive or not
882  double sign=-1.0;
883  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
884  {
885  sign=1.0;
886  }
887 
888  // if the second invariant is below the cutoff we have a constant, Newtonian
889  // viscosity
890  if(sign*second_invariant_of_rate_of_strain_tensor <
892  {
893  return a*pow(sign*second_invariant_of_rate_of_strain_tensor,3.0)+
894  b*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
895  zero_shear_viscosity;
896  }
897  else
898  {
899  // Calculate the square root of the absolute value of the
900  // second invariant of the rate of strain tensor
901  double measure_of_rate_of_strain=
902  sqrt(sign*second_invariant_of_rate_of_strain_tensor);
903 
904  return 1.0+(*Alpha_pt)*pow((2.0*measure_of_rate_of_strain),
905  *Flow_index_pt-1.0);
906  }
907  }
908 
909  /// Deriv of viscosity w.r.t. strain rate invariant
910  double dviscosity_dinvariant
911  (const double& second_invariant_of_rate_of_strain_tensor)
912  {
913  // Get the parameters of the cubic
914  double a;
915  double b;
916 
918 
919  // Pre-multiply the second invariant with +/-1 depending on whether it's
920  // positive or not
921  double sign=-1.0;
922  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
923  {
924  sign=1.0;
925  }
926 
927  if(sign*second_invariant_of_rate_of_strain_tensor <
929  {
930  return sign*3.0*a*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
931  2.0*b*second_invariant_of_rate_of_strain_tensor;
932  }
933  else
934  {
935  return (*Alpha_pt)*pow(2.0,(*Flow_index_pt)-2.0)*((*Flow_index_pt)-1.0)*
936  second_invariant_of_rate_of_strain_tensor*
937  pow(sign*second_invariant_of_rate_of_strain_tensor,
938  ((*Flow_index_pt)-1.0)/2.0-2.0);
939  }
940  }
941 
942 };
943 
944 //===================================================================
945 /// A GeneralisedNewtonianConstitutiveEquation class
946 /// defining a Casson model fluid
947 /// using Tanner and Milthorpe's (1983) regularisation
948 /// with a smooth transition using a cubic
949 //==================================================================
950 template<unsigned DIM>
953 {
954  private:
955 
956  /// Yield stress
958 
959  /// value of the second invariant below which we have
960  /// constant (Newtonian) viscosity -- assumed to be always positive
962 
963 
964  public:
965 
966 
967  /// "Cutoff regularised" Casson constitutive equation
969  (double* yield_stress_pt, double* critical_second_invariant_pt) :
971  Yield_stress_pt(yield_stress_pt),
972  Critical_second_invariant_pt(critical_second_invariant_pt)
973  {
974 
975  /// get the cutoff viscosity
976  double cut_off_viscosity=calculate_cutoff_viscosity();
977 
978  /// get the zero shear viscosity
979  double zero_shear_viscosity=calculate_zero_shear_viscosity();
980 
981  oomph_info << "CassonTanMilRegWithBlendingConstitutiveEquation: "
982  << " zero shear viscosity = " << zero_shear_viscosity
983  << std::endl;
984 
985  oomph_info << "CassonTanMilRegWithBlendingConstitutiveEquation: "
986  << " cut off viscosity = " << cut_off_viscosity
987  << std::endl;
988 
989  oomph_info << " "
990  << " cutoff invariant = "
992  << std::endl;
993 
994 
995  }
996 
997  /// Function that calculates the viscosity at the cutoff invariant
998  /// Note: this is NOT the viscosity at zero I2
1000  {
1001  // Calculate the Newtonian cutoff viscosity from the constitutive
1002  // equation and the cutoff value of the second invariant
1003  return (*Yield_stress_pt)/(2.0*sqrt(*Critical_second_invariant_pt))+
1004  2.0*sqrt(*Yield_stress_pt/(2.0*sqrt(*Critical_second_invariant_pt)))+1.0;
1005  }
1006 
1007  /// Offset by how much the zero shear rate viscosity lies
1008  /// above the viscosity at I2_cutoff
1009  /// Hard-coded to a value that ensures a smooth transition
1010  double calculate_viscosity_offset_at_zero_shear(double& cut_off_viscosity)
1011  {
1012  return cut_off_viscosity/5.0;
1013  }
1014 
1015  /// Function that calculates the viscosity at zero I2
1017  {
1018  // get the viscosity at the cutoff point
1019  double cut_off_viscosity=calculate_cutoff_viscosity();
1020 
1021  /// Offset by how much the zero shear rate viscosity lies
1022  /// above the viscosity at I2_cutoff
1023  double epsilon=calculate_viscosity_offset_at_zero_shear(cut_off_viscosity);
1024 
1025  return cut_off_viscosity+epsilon;
1026  }
1027 
1028  /// Report cutoff values
1029  void report_cut_off_values(double& cut_off_invariant,
1030  double& cut_off_viscosity,
1031  double& zero_shear_viscosity)
1032  {
1033  cut_off_invariant=*Critical_second_invariant_pt;
1034  cut_off_viscosity=calculate_cutoff_viscosity();
1035  zero_shear_viscosity=calculate_zero_shear_viscosity();
1036  }
1037 
1038  // Calculate the fitting parameters for the cubic blending
1039  void calculate_fitting_parameters_of_cubic(double& a, double& b)
1040  {
1041  // get the viscosity at the cutoff invariant
1042  double Cut_off_viscosity=calculate_cutoff_viscosity();
1043 
1044  // calculate the offset at zero shear
1045  double epsilon=calculate_viscosity_offset_at_zero_shear(Cut_off_viscosity);
1046 
1047  a=1.0/pow(*Critical_second_invariant_pt,39.0/4.0)*
1048  (pow(*Critical_second_invariant_pt,27.0/4.0)*
1049  (2.0*(Cut_off_viscosity+epsilon)-2.0*sqrt(2.0)*
1050  sqrt(*Yield_stress_pt/sqrt(*Critical_second_invariant_pt))-2.0)-
1051  5.0/4.0*(*Yield_stress_pt)*pow(*Critical_second_invariant_pt,25.0/4.0)-
1052  1.0/(2.0*sqrt(2.0))*sqrt(*Yield_stress_pt)*
1053  pow(*Critical_second_invariant_pt,13.0/2.0));
1054 
1055  b=1.0/pow(*Critical_second_invariant_pt,27.0/4.0)*
1056  (pow(*Critical_second_invariant_pt,19.0/4.0)*
1057  (-3.0*(Cut_off_viscosity+epsilon)+3.0*sqrt(2.0)*
1058  sqrt(*Yield_stress_pt/sqrt(*Critical_second_invariant_pt))+3.0)+
1059  7.0/4.0*(*Yield_stress_pt)*pow(*Critical_second_invariant_pt,17.0/4.0)+
1060  1.0/(2.0*sqrt(2.0))*sqrt(*Yield_stress_pt)*
1061  pow(*Critical_second_invariant_pt,9.0/2.0));
1062  }
1063 
1064 
1065  /// Viscosity ratio as a fct of strain rate invariant
1066  double viscosity(const double&
1067  second_invariant_of_rate_of_strain_tensor)
1068  {
1069  // Get the parameters of the cubic
1070  double a;
1071  double b;
1072 
1074 
1075  double zero_shear_viscosity=calculate_zero_shear_viscosity();
1076 
1077  // Pre-multiply the second invariant with +/-1 depending on whether it's
1078  // positive or not
1079  double sign=-1.0;
1080  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
1081  {
1082  sign=1.0;
1083  }
1084 
1085  // if the second invariant is below the cutoff we have a constant, Newtonian
1086  // viscosity
1087  if(sign*second_invariant_of_rate_of_strain_tensor <
1089  {
1090  return a*pow(sign*second_invariant_of_rate_of_strain_tensor,3.0)+
1091  b*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
1092  zero_shear_viscosity;
1093  }
1094  else
1095  {
1096  // Calculate the square root of the absolute value of the
1097  // second invariant of the rate of strain tensor
1098  double measure_of_rate_of_strain=
1099  sqrt(sign*second_invariant_of_rate_of_strain_tensor);
1100 
1101  return *Yield_stress_pt/(2.0*measure_of_rate_of_strain)+
1102  2.0*sqrt(*Yield_stress_pt/(2.0*measure_of_rate_of_strain))+1.0;
1103  }
1104  }
1105 
1106  /// Deriv of viscosity w.r.t. strain rate invariant
1107  double dviscosity_dinvariant
1108  (const double& second_invariant_of_rate_of_strain_tensor)
1109  {
1110  // Get the parameters of the cubic
1111  double a;
1112  double b;
1113 
1115 
1116  // Pre-multiply the second invariant with +/-1 depending on whether it's
1117  // positive or not
1118  double sign=-1.0;
1119  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
1120  {
1121  sign=1.0;
1122  }
1123 
1124  if(sign*second_invariant_of_rate_of_strain_tensor <
1126  {
1127  return sign*3.0*a*pow(sign*second_invariant_of_rate_of_strain_tensor,2.0)+
1128  2.0*b*second_invariant_of_rate_of_strain_tensor;
1129  }
1130  else
1131  {
1132  return -sqrt(*Yield_stress_pt)*
1133  second_invariant_of_rate_of_strain_tensor/
1134  (2.0*sqrt(2.0)*
1135  pow(sign*second_invariant_of_rate_of_strain_tensor,9.0/4.0))-
1136  (*Yield_stress_pt)*second_invariant_of_rate_of_strain_tensor/
1137  (4.0*pow(sign*second_invariant_of_rate_of_strain_tensor,5.0/2.0));
1138  }
1139  }
1140 
1141 };
1142 
1143 //===================================================================
1144 /// A GeneralisedNewtonianConstitutiveEquation class
1145 /// defining an arbitrary shear-thinning fluid
1146 //==================================================================
1147 template<unsigned DIM>
1150 {
1151  private:
1152 
1153  /// high shear rate viscosity
1154  double* Mu_inf_pt;
1155 
1156  /// zero shear rate viscosity
1157  double* Mu_0_pt;
1158 
1159  /// parameter that controls the steepness of the curve
1160  double* Alpha_pt;
1161 
1162 
1163  public:
1164 
1166  (double* mu_inf_pt, double* mu_0_pt, double* alpha_pt) :
1168  Mu_inf_pt(mu_inf_pt), Mu_0_pt(mu_0_pt), Alpha_pt(alpha_pt){}
1169 
1170  double viscosity(const double&
1171  second_invariant_of_rate_of_strain_tensor)
1172  {
1173  return (*Mu_inf_pt)+((*Mu_0_pt)-(*Mu_inf_pt))*
1174  exp(-(*Alpha_pt)*second_invariant_of_rate_of_strain_tensor);
1175  }
1176 
1177  double dviscosity_dinvariant
1178  (const double& second_invariant_of_rate_of_strain_tensor)
1179  {
1180  //std::ostringstream error_stream;
1181  //error_stream << "This has not been implemented yet!";
1182  //throw OomphLibError(
1183  // error_stream.str(),
1184  // OOMPH_CURRENT_FUNCTION,
1185  // OOMPH_EXCEPTION_LOCATION);
1186 
1187  //return 0.0;
1188 
1189  return (*Alpha_pt)*((*Mu_inf_pt)-(*Mu_0_pt))*
1190  exp(-(*Alpha_pt)*second_invariant_of_rate_of_strain_tensor);
1191  }
1192 
1193 };
1194 
1195 //===================================================================
1196 /// A GeneralisedNewtonianConstitutiveEquation class
1197 /// defining a fluid following a tanh-profile
1198 //==================================================================
1199 template<unsigned DIM>
1202 {
1203  private:
1204 
1205  /// high shear rate viscosity
1206  double* Mu_inf_pt;
1207 
1208  /// zero shear rate viscosity
1209  double* Mu_0_pt;
1210 
1211  /// parameter controlling the steepness of the step
1212  /// nb -- I used 10.0/(*Critical_second_invariant_pt)
1213  double* Alpha_pt;
1214 
1215  /// parameter that controls the location of the step -- assumed
1216  /// to be always positive
1218 
1219  public:
1220 
1222  (double* mu_inf_pt, double* mu_0_pt, double* alpha_pt,
1223  double* critical_second_invariant_pt) :
1225  Mu_inf_pt(mu_inf_pt), Mu_0_pt(mu_0_pt), Alpha_pt(alpha_pt),
1226  Critical_second_invariant_pt(critical_second_invariant_pt){}
1227 
1228  double viscosity(const double&
1229  second_invariant_of_rate_of_strain_tensor)
1230  {
1231  // Pre-multiply the second invariant with +/-1 depending on whether it's
1232  // positive or not
1233  double sign=-1.0;
1234  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
1235  {
1236  sign=1.0;
1237  }
1238 
1239  return ((*Mu_0_pt)-(*Mu_inf_pt))/2.0*
1240  tanh(((*Critical_second_invariant_pt)-sign*
1241  second_invariant_of_rate_of_strain_tensor)*(*Alpha_pt))-
1242  ((*Mu_0_pt)-(*Mu_inf_pt))/2.0+(*Mu_0_pt);
1243  }
1244 
1245  double dviscosity_dinvariant
1246  (const double& second_invariant_of_rate_of_strain_tensor)
1247  {
1248  // Pre-multiply the second invariant with +/-1 depending on whether it's
1249  // positive or not
1250  double sign=-1.0;
1251  if(second_invariant_of_rate_of_strain_tensor >= 0.0)
1252  {
1253  sign=1.0;
1254  }
1255 
1256  return -sign*((*Mu_0_pt)-(*Mu_inf_pt))*10.0/
1258  1.0/cosh(((*Critical_second_invariant_pt)-sign*
1259  second_invariant_of_rate_of_strain_tensor)*(*Alpha_pt))*
1260  1.0/cosh(((*Critical_second_invariant_pt)-sign*
1261  second_invariant_of_rate_of_strain_tensor)*(*Alpha_pt));
1262  }
1263 
1264 };
1265 
1266 
1267 /////////////////////////////////////////////////////////////////////
1268 /////////////////////////////////////////////////////////////////////
1269 /////////////////////////////////////////////////////////////////////
1270 
1271 }
1272 
1273 
1274 #endif
A Base class defining the generalise Newtonian constitutive relation.
PowerLawBerEngRegConstitutiveEquation(double *power_pt, double *reg_par_pt)
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
Viscosity ratio as a fct of strain rate invariant.
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
double calculate_zero_shear_viscosity()
Function that calculates the viscosity at zero I2.
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
Deriv of viscosity w.r.t. strain rate invariant.
double * Alpha_pt
parameter that controls the steepness of the curve
CassonTanMilRegWithBlendingConstitutiveEquation(double *yield_stress_pt, double *critical_second_invariant_pt)
"Cutoff regularised" Casson constitutive equation
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
Viscosity ratio as a fct of strain rate invariant.
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
in the Newtonian case the viscosity is constant
double * Regularisation_parameter_pt
regularisation parameter e << 1
virtual double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)=0
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
Deriv of viscosity w.r.t. strain rate invariant.
TanhProfileConstitutiveEquation(double *mu_inf_pt, double *mu_0_pt, double *alpha_pt, double *critical_second_invariant_pt)
void report_cut_off_values(double &cut_off_invariant, double &cut_off_viscosity, double &zero_shear_viscosity)
Report cutoff values.
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
Viscosity ratio as a fct of strain rate invariant.
HerschelBulkleyTanMilRegConstitutiveEquation(double *yield_stress_pt, double *flow_index_pt, double *critical_second_invariant_pt)
"Cutoff regularised" Herschel Bulkley constitutive equation
OomphInfo oomph_info
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
Viscosity ratio as a fct of strain rate invariant.
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
NicosConstitutiveEquation(double *mu_inf_pt, double *mu_0_pt, double *alpha_pt)
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
Deriv of viscosity w.r.t. strain rate invariant.
void report_cut_off_values(double &cut_off_invariant, double &cut_off_viscosity, double &zero_shear_viscosity)
Report cutoff values.
double calculate_cut_off_viscosity()
Function that calculates the cut off viscosity.
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
Deriv of viscosity w.r.t. strain rate invariant.
double calculate_zero_shear_viscosity()
Function that calculates the viscosity at zero I2.
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
Viscosity ratio as a fct of strain rate invariant.
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
Deriv of viscosity w.r.t. strain rate invariant.
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
the derivative w.r.t. I2 is zero
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
SiskoTanMilRegWithBlendingConstitutiveEquation(double *alpha_pt, double *flow_index_pt, double *critical_second_invariant_pt)
"Cutoff regularised" Sisko constitutive equation
HerschelBulkleyBerEngRegConstitutiveEquation(double *yield_stress_pt, double *flow_index_pt, double *regularisation_parameter_pt)
NewtonianConstitutiveEquation(const double &viscosity_ratio=1.0)
Constructor: specify viscosity ratio (defaults to one)
HerschelBulkleyPapRegConstitutiveEquation(double *yield_stress_pt, double *flow_index_pt, double *exponential_parameter_pt)
double calculate_zero_shear_viscosity()
Function that calculates the viscosity at zero I2.
void report_cut_off_values(double &cut_off_invariant, double &cut_off_viscosity, double &zero_shear_viscosity)
Report cutoff values.
HerschelBulkleyMenDutRegConstitutiveEquation(double *yield_stress_pt, double *flow_index_pt, double *zero_shear_viscosity_pt)
"Exponentially regularised" Herschel Bulkley constitutive equation
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
double viscosity(const double &second_invariant_of_rate_of_strain_tensor)
HerschelBulkleyTanMilRegWithBlendingConstitutiveEquation(double *yield_stress_pt, double *flow_index_pt, double *critical_second_invariant_pt)
"Cutoff regularised" Herschel Bulkley constitutive equation
virtual double viscosity(const double &second_invariant_of_rate_of_strain_tensor)=0
double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)
const double eps
Definition: orthpoly.h:57
void report_cut_off_values(double &cut_off_invariant, double &cut_off_viscosity)
Report cutoff values.