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