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-2022 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
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
30
31
32// Oomph-lib includes
33//#include "generic.h"
34
35namespace 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),
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
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
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
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) *
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) +
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
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
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
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 =
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 =
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) *
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) *
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 =
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 =
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 *
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 /
1207 sign * second_invariant_of_rate_of_strain_tensor) *
1208 (*Alpha_pt)) *
1209 1.0 /
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.
void calculate_fitting_parameters_of_cubic(double &a, double &b)
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...
GeneralisedNewtonianConstitutiveEquation()
Empty constructor.
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 ...
virtual ~GeneralisedNewtonianConstitutiveEquation()
Empty virtual destructor.
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 * Regularisation_parameter_pt
regularisation parameter e << 1
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 * Zero_shear_viscosity_pt
the viscosity at zero shear rate
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...
double * Exponential_parameter_pt
Regularisation parameter m >> 1.
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 Viscosity_ratio
Viscosity ratio.
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 ...
double * Mu_0_pt
zero shear rate viscosity
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
double * Mu_inf_pt
high shear rate viscosity
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)
double * Regularisation_parameter_pt
regularisation parameter e << 1
A GeneralisedNewtonianConstitutiveEquation class defining a Sisko fluid using Tanner and Milthorpe's ...
void calculate_fitting_parameters_of_cubic(double &a, double &b)
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 * Mu_inf_pt
high shear rate viscosity
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 * Mu_0_pt
zero shear rate viscosity
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...