constitutive_laws.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 // Header file for ConstitutiveLaw objects that will be used in all
27 // elasticity-type elements
28 
29 #ifndef OOMPH_CONSTITUTIVE_LAWS_HEADER
30 #define OOMPH_CONSTITUTIVE_LAWS_HEADER
31 
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 // OOMPH-LIB includes
38 #include "../generic/oomph_utilities.h"
39 #include "../generic/matrices.h"
40 
41 namespace oomph
42 {
43  //=====================================================================
44  /// Base class for strain energy functions to be used in solid
45  /// mechanics computations.
46  //====================================================================
48  {
49  public:
50  /// Constructor takes no arguments
52 
53  /// Empty virtual destructor
54  virtual ~StrainEnergyFunction() {}
55 
56 
57  /// Return the strain energy in terms of the strain tensor
58  virtual double W(const DenseMatrix<double>& gamma)
59  {
60  std::string error_message =
61  "The strain-energy function as a function of the strain-tensor,\n";
62  error_message +=
63  "gamma, is not implemented for this strain energy function.\n";
64 
65  throw OomphLibError(
66  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
67  return 0.0;
68  }
69 
70 
71  /// Return the strain energy in terms of the strain invariants
72  virtual double W(const Vector<double>& I)
73  {
74  std::string error_message =
75  "The strain-energy function as a function of the strain\n ";
76  error_message +=
77  "invariants, I1, I2, I3, is not implemented for this strain\n ";
78  error_message += "energy function\n";
79 
80  throw OomphLibError(
81  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
82  return 0.0;
83  }
84 
85 
86  /// Return the derivatives of the strain energy function with
87  /// respect to the components of the strain tensor (default is to use
88  /// finite differences).
89  virtual void derivative(const DenseMatrix<double>& gamma,
90  DenseMatrix<double>& dWdgamma)
91  {
92  throw OomphLibError(
93  "Sorry, the FD setup of dW/dgamma hasn't been implemented yet",
94  OOMPH_CURRENT_FUNCTION,
95  OOMPH_EXCEPTION_LOCATION);
96  }
97 
98 
99  /// Return the derivatives of the strain energy function with
100  /// respect to the strain invariants. Default version is to use finite
101  /// differences
103  {
104  // Calculate the derivatives of the strain-energy-function wrt the strain
105  // invariants
106  double FD_Jstep = 1.0e-8; // Usual comments about global stuff
107  double energy = W(I);
108 
109  // Loop over the strain invariants
110  for (unsigned i = 0; i < 3; i++)
111  {
112  // Store old value
113  double I_prev = I[i];
114  // Increase ith strain invariant
115  I[i] += FD_Jstep;
116  // Get the new value of the strain energy
117  double energy_new = W(I);
118  // Calculate the value of the derivative
119  dWdI[i] = (energy_new - energy) / FD_Jstep;
120  // Reset value of ith strain invariant
121  I[i] = I_prev;
122  }
123  }
124 
125  /// Pure virtual function in which the user must declare if the
126  /// constitutive equation requires an incompressible formulation
127  /// in which the volume constraint is enforced explicitly.
128  /// Used as a sanity check in PARANOID mode.
130  };
131 
132 
133  /// /////////////////////////////////////////////////////////////////
134  /// /////////////////////////////////////////////////////////////////
135  /// /////////////////////////////////////////////////////////////////
136 
137 
138  //=====================================================================
139  /// MooneyRivlin strain-energy function.
140  /// with constitutive parameters C1 and C2:
141  /// \f[ W = C_1 (I_0 - 3) + C_2 (I_1 - 3) \f]
142  /// where incompressibility (\f$ I_2 \equiv 1\f$) is assumed.
143  //====================================================================
145  {
146  public:
147  /// Constructor takes the pointer to the value of the constants
148  MooneyRivlin(double* c1_pt, double* c2_pt)
149  : StrainEnergyFunction(), C1_pt(c1_pt), C2_pt(c2_pt)
150  {
151  }
152 
153  /// Empty Virtual destructor
154  virtual ~MooneyRivlin() {}
155 
156  /// Return the strain energy in terms of strain tensor
157  double W(const DenseMatrix<double>& gamma)
158  {
159  return StrainEnergyFunction::W(gamma);
160  }
161 
162  /// Return the strain energy in terms of the strain invariants
163  double W(const Vector<double>& I)
164  {
165  return (*C1_pt) * (I[0] - 3.0) + (*C2_pt) * (I[1] - 3.0);
166  }
167 
168 
169  /// Return the derivatives of the strain energy function with
170  /// respect to the strain invariants
172  {
173  dWdI[0] = (*C1_pt);
174  dWdI[1] = (*C2_pt);
175  dWdI[2] = 0.0;
176  }
177 
178  /// Pure virtual function in which the user must declare if the
179  /// constitutive equation requires an incompressible formulation
180  /// in which the volume constraint is enforced explicitly.
181  /// Used as a sanity check in PARANOID mode. True
183  {
184  return true;
185  }
186 
187 
188  private:
189  /// Pointer to first Mooney Rivlin constant
190  double* C1_pt;
191 
192  /// Pointer to second Mooney Rivlin constant
193  double* C2_pt;
194  };
195 
196 
197  /// /////////////////////////////////////////////////////////////////
198  /// /////////////////////////////////////////////////////////////////
199  /// /////////////////////////////////////////////////////////////////
200 
201 
202  //=====================================================================
203  /// Generalisation of Mooney Rivlin constitutive law to compressible
204  /// media as suggested on p. 553 of Fung, Y.C. & Tong, P. "Classical and
205  /// Computational Solid Mechanics" World Scientific (2001).
206  /// Input parameters are Young's modulus E, Poisson ratio nu and
207  /// the Mooney-Rivlin constant C1. In the small-deformation-limit
208  /// the behaviour becomes equivalent to that of linear elasticity
209  /// with the same E and nu.
210  ///
211  /// Note that there's a factor of 2 difference between C1 and the Mooney
212  /// Rivlin C1!
213  //====================================================================
215  {
216  public:
217  /// Constructor takes the pointers to the constitutive parameters:
218  /// Poisson's ratio, the Mooney-Rivlin parameter. Young's modulus is set
219  /// to 1, implying that it has been used to scale the stresses
220  GeneralisedMooneyRivlin(double* nu_pt, double* c1_pt)
222  Nu_pt(nu_pt),
223  C1_pt(c1_pt),
224  E_pt(new double(1.0)),
225  Must_delete_e(true)
226  {
227  }
228 
229  /// Constructor takes the pointers to the constitutive parameters:
230  /// Poisson's ratio, the Mooney-Rivlin parameter and Young's modulus
231  GeneralisedMooneyRivlin(double* nu_pt, double* c1_pt, double* e_pt)
233  Nu_pt(nu_pt),
234  C1_pt(c1_pt),
235  E_pt(e_pt),
236  Must_delete_e(false)
237  {
238  }
239 
240 
241  /// Virtual destructor
243  {
244  if (Must_delete_e) delete E_pt;
245  }
246 
247  /// Return the strain energy in terms of strain tensor
248  double W(const DenseMatrix<double>& gamma)
249  {
250  return StrainEnergyFunction::W(gamma);
251  }
252 
253 
254  /// Return the strain energy in terms of the strain invariants
255  double W(const Vector<double>& I)
256  {
257  double G = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
258  return 0.5 * ((*C1_pt) * (I[0] - 3.0) + (G - (*C1_pt)) * (I[1] - 3.0) +
259  ((*C1_pt) - 2.0 * G) * (I[2] - 1.0) +
260  (1.0 - (*Nu_pt)) * G * (I[2] - 1.0) * (I[2] - 1.0) /
261  (2.0 * (1.0 - 2.0 * (*Nu_pt))));
262  }
263 
264 
265  /// Return the derivatives of the strain energy function with
266  /// respect to the strain invariants
268  {
269  double G = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
270  dWdI[0] = 0.5 * (*C1_pt);
271  dWdI[1] = 0.5 * (G - (*C1_pt));
272  dWdI[2] = 0.5 * ((*C1_pt) - 2.0 * G +
273  2.0 * (1.0 - (*Nu_pt)) * G * (I[2] - 1.0) /
274  (2.0 * (1.0 - 2.0 * (*Nu_pt))));
275  }
276 
277 
278  /// Pure virtual function in which the user must declare if the
279  /// constitutive equation requires an incompressible formulation
280  /// in which the volume constraint is enforced explicitly.
281  /// Used as a sanity check in PARANOID mode. False.
283  {
284  return false;
285  }
286 
287  private:
288  /// Poisson's ratio
289  double* Nu_pt;
290 
291  /// Mooney-Rivlin parameter
292  double* C1_pt;
293 
294  /// Young's modulus
295  double* E_pt;
296 
297  /// Boolean flag to indicate if storage for elastic modulus
298  /// must be deleted in destructor
300  };
301 
302 
303  /// /////////////////////////////////////////////////////////////////
304  /// /////////////////////////////////////////////////////////////////
305  /// /////////////////////////////////////////////////////////////////
306 
307 
308  //===========================================================================
309  /// A class for constitutive laws for elements that solve
310  /// the equations of solid mechanics based upon the principle of virtual
311  /// displacements. In that formulation, the information required from a
312  /// constitutive law is the (2nd Piola-Kirchhoff) stress tensor
313  /// \f$ \sigma^{ij} \f$ as a function of the (Green) strain
314  /// \f$ \gamma^{ij} \f$:
315  /// \f[ \sigma^{ij} = \sigma^{ij}(\gamma_{ij}). \f]
316  /// The Green strain is defined as
317  /// \f[ \gamma_{ij} = \frac{1}{2} (G_{ij} - g_{ij}), \ \ \ \ \ \ \ \ \ \ \ (1) \f]
318  /// where \f$G_{ij} \f$ and \f$ g_{ij}\f$ are the metric tensors
319  /// in the deformed and undeformed (stress-free) configurations, respectively.
320  /// A specific ConstitutiveLaw needs to be implement the pure
321  /// virtual function
322  /// \code
323  /// ConstitutiveLaw::calculate_second_piola_kirchhoff_stress(...)
324  /// \endcode
325  /// Equation (1) shows that the strain may be calculated from the
326  /// undeformed and deformed metric tensors. Frequently, these tensors are
327  /// also required in the constitutive law itself.
328  /// To avoid unnecessary re-computation of these quantities, we
329  /// pass the deformed and undeformed metric tensor to
330  /// \c calculate_second_piola_kirchhoff_stress(...)
331  /// rather than the strain tensor itself.
332  ///
333  /// The functional form of the constitutive equation is different
334  /// for compressible/incompressible/near-incompressible behaviour
335  /// and we provide interfaces that are appropriate for all of these cases.
336  /// -# \b Compressible \b Behaviour: \n If the material is compressible,
337  /// the stress can be computed from the deformed and undeformed
338  /// metric tensors,
339  ///
340  /// \f[ \sigma^{ij} = \sigma^{ij}(\gamma_{ij}) = \sigma^{ij}\bigg( \frac{1}{2} (G_{ij} - g_{ij})\bigg), \f]
341  /// using the interface
342  /// \code
343  /// // 2nd Piola Kirchhoff stress tensor
344  /// DenseMatrix<double> sigma(DIM,DIM);
345  ///
346  /// // Metric tensor in the undeformed (stress-free) configuration
347  /// DenseMatrix<double> g(DIM,DIM);
348  ///
349  /// // Metric tensor in the deformed configuration
350  /// DenseMatrix<double> G(DIM,DIM);
351  ///
352  /// // Compute stress from the two metric tensors:
353  /// calculate_second_piola_kirchhoff_stress(g,G,sigma);
354  /// \endcode
355  /// \n \n \n
356  /// -# \b Incompressible \b Behaviour: \n If the material is incompressible,
357  /// its deformation is constrained by the condition that
358  ///
359  /// \f[ \det G_{ij} - \det g_{ij}= 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2) \f]
360  /// which ensures that the volume of infinitesimal material
361  /// elements remains constant during the deformation. This
362  /// condition is typically enforced by a Lagrange multiplier which
363  /// plays the role of a pressure. In such cases, the
364  /// stress tensor has form
365  ///
366  /// \f[ \sigma^{ij} = -p G^{ij} + \overline{\sigma}^{ij}\big(\gamma_{kl}\big), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3) \f]
367  /// where only the deviatoric part of the stress tensor,
368  /// \f$ \overline{\sigma}^{ij}, \f$ depends directly on the
369  /// strain. The pressure \f$ p \f$ needs to be determined
370  /// independently from (2).
371  /// Given the deformed and undeformed metric tensors,
372  /// the computation of the stress tensor \f$ \sigma^{ij} \f$
373  /// for an incompressible
374  /// material therefore requires the computation of the following
375  /// quantities:
376  /// - The deviatoric stress \f$ \overline{\sigma}^{ij} \f$
377  /// - The contravariant deformed metric tensor \f$ G^{ij} \f$
378  /// - The determinant of the deformed
379  /// metric tensors, \f$ \det G_{ij}, \f$ which
380  /// is required in equation (2) whose solution determines the pressure.
381  /// .
382  /// \n
383  /// These quantities can be obtained from the following interface \n
384  /// \code
385  /// // Deviatoric part of the 2nd Piola Kirchhoff stress tensor
386  /// DenseMatrix<double> sigma_dev(DIM,DIM);
387  ///
388  /// // Metric tensor in the undeformed (stress-free) configuration
389  /// DenseMatrix<double> g(DIM,DIM);
390  ///
391  /// // Metric tensor in the deformed configuration
392  /// DenseMatrix<double> G(DIM,DIM);
393  ///
394  /// // Determinant of the deformed metric tensor
395  /// double Gdet;
396  ///
397  /// // Contravariant deformed metric tensor
398  /// DenseMatrix<double> G_contra(DIM,DIM);
399  ///
400  /// // Compute stress from the two metric tensors:
401  /// calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,G_contra,Gdet);
402  /// \endcode
403  /// \n \n \n
404  /// -# \b Nearly \b Incompressible \b Behaviour: \n If the material is nearly
405  /// incompressible, it is advantageous to split the stress into
406  /// its deviatoric and hydrostatic parts by writing the
407  /// constitutive law in the form
408  ///
409  /// \f[ \sigma^{ij} = -p G^{ij} + \overline{\sigma}^{ij}\big(\gamma_{kl}\big), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3) \f]
410  /// where the deviatoric part of the stress tensor,
411  /// \f$ \overline{\sigma}^{ij}, \f$ depends on the
412  /// strain. This form of the constitutive
413  /// law is identical to that of the incompressible
414  /// case and it involves a pressure \f$ p \f$ which needs to be
415  /// determined from an additional equation. In the
416  /// incompressible case, this equation was given by the incompressibility
417  /// constraint (2). Here, we need to augment the constitutive law (3) by
418  /// a separate equation for the pressure. Generally this takes the
419  /// form
420  ///
421  /// \f[ p = - \kappa \ d \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4) \f]
422  /// where \f$ \kappa \f$ is the "bulk modulus", a material property
423  /// that needs to be specified by the constitutive law.
424  /// \f$ d \f$ is the (generalised) dilatation, i.e. the relative change
425  /// in the volume of an infinitesimal material element (or some
426  /// suitable generalised quantitiy that is related to it). As the
427  /// material approaches incompressibility, \f$ \kappa \to \infty\f$, so
428  /// that infinitely large pressures would be required to achieve any change
429  /// in volume. To facilitate the implementation of (4) as the equation for
430  /// the pressure, we re-write it in the form
431  /// \f[ p \ \frac{1}{\kappa} + d\big(g_{ij},G_{ij}\big) = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (5) \f]
432  /// which only involves quantities that remain finite
433  /// as we approach true incompressibility.
434  /// \n
435  /// Given the deformed and undeformed metric tensors,
436  /// the computation of the stress tensor \f$ \sigma^{ij} \f$
437  /// for a nearly incompressible
438  /// material therefore requires the computation of the following
439  /// quantities:
440  /// - The deviatoric stress \f$ \overline{\sigma}^{ij} \f$
441  /// - The contravariant deformed metric tensor \f$ G^{ij} \f$
442  /// - The generalised dilatation \f$ d \f$
443  /// - The inverse of the bulk modulus \f$ \kappa \f$
444  /// .
445  /// \n
446  /// These quantities can be obtained from the following interface \n
447  /// \code
448  /// // Deviatoric part of the 2nd Piola Kirchhoff stress tensor
449  /// DenseMatrix<double> sigma_dev(DIM,DIM);
450  ///
451  /// // Metric tensor in the undeformed (stress-free) configuration
452  /// DenseMatrix<double> g(DIM,DIM);
453  ///
454  /// // Metric tensor in the deformed configuration
455  /// DenseMatrix<double> G(DIM,DIM);
456  ///
457  /// // Contravariant deformed metric tensor
458  /// DenseMatrix<double> G_contra(DIM,DIM);
459  ///
460  /// // Inverse of the bulk modulus
461  /// double inv_kappa;
462  ///
463  /// // Generalised dilatation
464  /// double gen_dil;
465  ///
466  /// // Compute stress from the two metric tensors:
467  /// calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,G_contra,inv_kappa,gen_dil);
468  /// \endcode
469  //==========================================================================
471  {
472  protected:
473  /// Test whether a matrix is square
474  bool is_matrix_square(const DenseMatrix<double>& M);
475 
476  /// Test whether two matrices are of equal dimensions
478  const DenseMatrix<double>& M2);
479 
480  /// Check for errors in the input,
481  /// i.e. check that the dimensions of the arrays are all consistent
483  const DenseMatrix<double>& G,
484  DenseMatrix<double>& sigma);
485 
486  /// Calculate a contravariant tensor from a covariant tensor,
487  /// and return the determinant of the covariant tensor.
488  double calculate_contravariant(const DenseMatrix<double>& Gcov,
489  DenseMatrix<double>& Gcontra);
490 
491  /// Calculate the derivatives of the contravariant tensor
492  /// and the derivatives of the determinant of the covariant tensor
493  /// with respect to the components of the covariant tensor
495  RankFourTensor<double>& dGcontra_dG,
496  DenseMatrix<double>& d_detG_dG);
497 
498 
499  public:
500  /// Empty constructor
502 
503 
504  /// Empty virtual destructor
505  virtual ~ConstitutiveLaw() {}
506 
507 
508  /// Calculate the contravariant 2nd Piola Kirchhoff
509  /// stress tensor. Arguments are the
510  /// covariant undeformed and deformed metric tensor and the
511  /// matrix in which to return the stress tensor
513  const DenseMatrix<double>& g,
514  const DenseMatrix<double>& G,
515  DenseMatrix<double>& sigma) = 0;
516 
517  /// Calculate the derivatives of the contravariant
518  /// 2nd Piola Kirchhoff stress tensor with respect to the deformed metric
519  /// tensor. Arguments are the
520  /// covariant undeformed and deformed metric tensor, the current value of
521  /// the stress tensor and the
522  /// rank four tensor in which to return the derivatives of the stress tensor
523  /// The default implementation uses finite differences, but can be
524  /// overloaded for constitutive laws in which an analytic formulation
525  /// is possible.
526  /// If the boolean flag symmetrize_tensor is false, only the
527  /// "upper triangular" entries of the tensor will be filled in. This is
528  /// a useful efficiency when using the derivatives in Jacobian calculations.
530  const DenseMatrix<double>& g,
531  const DenseMatrix<double>& G,
532  const DenseMatrix<double>& sigma,
533  RankFourTensor<double>& d_sigma_dG,
534  const bool& symmetrize_tensor = true);
535 
536 
537  /// Calculate the deviatoric part
538  /// \f$ \overline{ \sigma^{ij}}\f$ of the contravariant
539  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
540  /// Also return the contravariant deformed metric
541  /// tensor and the determinant of the deformed metric tensor.
542  /// This form is appropriate
543  /// for truly-incompressible materials for which
544  /// \f$ \sigma^{ij} = - p G^{ij} +\overline{ \sigma^{ij}} \f$
545  /// where the "pressure" \f$ p \f$ is determined by
546  /// \f$ \det G_{ij} - \det g_{ij} = 0 \f$.
548  const DenseMatrix<double>& g,
549  const DenseMatrix<double>& G,
550  DenseMatrix<double>& sigma_dev,
551  DenseMatrix<double>& G_contra,
552  double& Gdet)
553  {
554  throw OomphLibError(
555  "Incompressible formulation not implemented for this constitutive law",
556  OOMPH_CURRENT_FUNCTION,
557  OOMPH_EXCEPTION_LOCATION);
558  }
559 
560  /// Calculate the derivatives of the contravariant
561  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
562  /// with respect to the deformed metric tensor.
563  /// Also return the derivatives of the determinant of the
564  /// deformed metric tensor with respect to the deformed metric tensor.
565  /// This form is appropriate
566  /// for truly-incompressible materials.
567  /// The default implementation uses finite differences for the
568  /// derivatives that depend on the constitutive law, but not
569  /// for the derivatives of the determinant, which are generic.
570  /// / If the boolean flag symmetrize_tensor is false, only the
571  /// "upper triangular" entries of the tensor will be filled in. This is
572  /// a useful efficiency when using the derivatives in Jacobian calculations.
574  const DenseMatrix<double>& g,
575  const DenseMatrix<double>& G,
576  const DenseMatrix<double>& sigma,
577  const double& detG,
578  const double& interpolated_solid_p,
579  RankFourTensor<double>& d_sigma_dG,
580  DenseMatrix<double>& d_detG_dG,
581  const bool& symmetrize_tensor = true);
582 
583 
584  /// Calculate the deviatoric part of the contravariant
585  /// 2nd Piola Kirchoff stress tensor. Also return the contravariant
586  /// deformed metric tensor, the generalised dilatation, \f$ d, \f$ and
587  /// the inverse of the bulk modulus \f$ \kappa\f$. This form is appropriate
588  /// for near-incompressible materials for which
589  /// \f$ \sigma^{ij} = -p G^{ij} + \overline{ \sigma^{ij}} \f$
590  /// where the "pressure" \f$ p \f$ is determined from
591  /// \f$ p / \kappa - d =0 \f$.
593  const DenseMatrix<double>& g,
594  const DenseMatrix<double>& G,
595  DenseMatrix<double>& sigma_dev,
596  DenseMatrix<double>& Gcontra,
597  double& gen_dil,
598  double& inv_kappa)
599  {
600  throw OomphLibError(
601  "Near-incompressible formulation not implemented for constitutive law",
602  OOMPH_CURRENT_FUNCTION,
603  OOMPH_EXCEPTION_LOCATION);
604  }
605 
606  /// Calculate the derivatives of the contravariant
607  /// 2nd Piola Kirchoff stress tensor with respect to the deformed metric
608  /// tensor. Also return the derivatives of the generalised dilatation,
609  /// \f$ d, \f$ with respect to the deformed metric tensor.
610  /// This form is appropriate
611  /// for near-incompressible materials.
612  /// The default implementation uses finite differences.
613  /// If the boolean flag symmetrize_tensor is false, only the
614  /// "upper triangular" entries of the tensor will be filled in. This is
615  /// a useful efficiency when using the derivatives in Jacobian calculations.
617  const DenseMatrix<double>& g,
618  const DenseMatrix<double>& G,
619  const DenseMatrix<double>& sigma,
620  const double& gen_dil,
621  const double& inv_kappa,
622  const double& interpolated_solid_p,
623  RankFourTensor<double>& d_sigma_dG,
624  DenseMatrix<double>& d_gen_dil_dG,
625  const bool& symmetrize_tensor = true);
626 
627 
628  /// Pure virtual function in which the user must declare if the
629  /// constitutive equation requires an incompressible formulation
630  /// in which the volume constraint is enforced explicitly.
631  /// Used as a sanity check in PARANOID mode.
633  };
634 
635 
636  /// //////////////////////////////////////////////////////////////////////
637  /// //////////////////////////////////////////////////////////////////////
638  /// //////////////////////////////////////////////////////////////////////
639 
640 
641  //========================================================================
642  /// Class for a "non-rational" extension of classical linear elasticity
643  /// to large displacements:
644  /// \f[ \sigma^{ij} = E^{ijkl} \gamma_{kl} \f]
645  /// where
646  /// \f[ E^{ijkl} = \frac{E}{(1+\nu)} \left( \frac{\nu}{(1-2\nu)} G^{ij} G^{kl} + \frac{1}{2} \left( G^{ik} G^{jl} + G^{il} G^{jk} \right) \right) \f]
647  /// For small strains \f$ (| G_{ij} - g_{ij} | \ll 1)\f$ this approaches
648  /// the version appropriate for linear elasticity, obtained
649  /// by replacing \f$ G^{ij}\f$ with \f$ g^{ij}\f$.
650  ///
651  /// We provide three versions of \c calculate_second_piola_kirchhoff_stress():
652  /// -# If \f$ \nu \ne 1/2 \f$ (and not close to \f$ 1/2 \f$), the
653  /// constitutive law can be used directly in the above form, using
654  /// the deformed and undeformed metric tensors as input.
655  /// -# If the material is incompressible (\f$ \nu = 1/2 \f$),
656  /// the first term in the above expression for \f$ E^{ijkl} \f$
657  /// is singular. We re-write the constitutive equation for this
658  /// case as
659  ///
660  /// \f[ \sigma^{ij} = -p G^{ij} + \frac{E}{3} \left( G^{ik} G^{jl} + G^{il} G^{jk} \right) \gamma_{kl} \f]
661  /// where the pressure \f$ p \f$ needs to be determined independently
662  /// via the incompressibility constraint.
663  /// In this case, the stress returned by
664  /// \c calculate_second_piola_kirchhoff_stress()
665  /// contains only the deviatoric part of the 2nd Piola Kirchhoff stress,
666  ///
667  /// \f[ \overline{\sigma}^{ij} = \frac{E}{3} \left( G^{ik} G^{jl} + G^{il} G^{jk} \right) \gamma_{kl}. \f]
668  /// The function also returns the contravariant metric tensor
669  /// \f$ G^{ij}\f$ (since it is needed to form the complete stress
670  /// tensor), and the determinant of the deformed covariant metric
671  /// tensor \f$ {\tt detG} = \det G_{ij} \f$ (since it is needed
672  /// in the equation that enforces the incompressibility).
673  /// -# If \f$ \nu \approx 1/2 \f$, the original form of the
674  /// constitutive equation could be used, but the resulting
675  /// equations tend to be ill-conditioned since they contain
676  /// the product of the large "bulk modulus"
677  ///
678  /// \f[ \kappa = \frac{E\nu}{(1+\nu)(1-2\nu)} \f]
679  /// and the small "generalised dilatation"
680  ///
681  /// \f[ d = \frac{1}{2} G^{ij} (G_{ij}-g_{ij}). \f]
682  /// [\f$ d \f$ represents the actual dilatation in the small
683  /// strain limit; for large deformations it doesn't have
684  /// any sensible interpretation (or does it?). It is simply
685  /// the term that needs to go to zero as \f$ \kappa \to \infty\f$.]
686  /// In this case, the stress returned by
687  /// \c calculate_second_piola_kirchhoff_stress()
688  /// contains only the deviatoric part of the 2nd Piola Kirchhoff stress,
689  ///
690  /// \f[ \overline{\sigma}^{ij} = \frac{E}{3} \left( G^{ik} G^{jl} + G^{il} G^{jk} \right) \gamma_{kl}. \f]
691  /// The function also returns the contravariant metric tensor
692  /// \f$ G^{ij}\f$ (since it is needed to form the complete stress
693  /// tensor), the inverse of the bulk modulus, and the generalised
694  /// dilatation (since they are needed in the equation
695  /// that determines the pressure).
696  ///
697  //=========================================================================
699  {
700  public:
701  /// The constructor takes the pointers to values of material parameters:
702  /// Poisson's ratio and Young's modulus.
703  GeneralisedHookean(double* nu_pt, double* e_pt)
704  : ConstitutiveLaw(), Nu_pt(nu_pt), E_pt(e_pt), Must_delete_e(false)
705  {
706  }
707 
708  /// The constructor takes the pointers to value of
709  /// Poisson's ratio . Young's modulus is set to E=1.0,
710  /// implying that all stresses have been non-dimensionalised
711  /// on on it.
712  GeneralisedHookean(double* nu_pt)
713  : ConstitutiveLaw(),
714  Nu_pt(nu_pt),
715  E_pt(new double(1.0)),
716  Must_delete_e(true)
717  {
718  }
719 
720 
721  /// Virtual destructor
723  {
724  if (Must_delete_e) delete E_pt;
725  }
726 
727  /// Calculate the contravariant 2nd Piola Kirchhoff
728  /// stress tensor. Arguments are the
729  /// covariant undeformed and deformed metric tensor and the
730  /// matrix in which to return the stress tensor
732  const DenseMatrix<double>& G,
733  DenseMatrix<double>& sigma);
734 
735 
736  /// Calculate the deviatoric part
737  /// \f$ \overline{ \sigma^{ij}}\f$ of the contravariant
738  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
739  /// Also return the contravariant deformed metric
740  /// tensor and the determinant of the deformed metric tensor.
741  /// This form is appropriate
742  /// for truly-incompressible materials for which
743  /// \f$ \sigma^{ij} = - p G^{ij} +\overline{ \sigma^{ij}} \f$
744  /// where the "pressure" \f$ p \f$ is determined by
745  /// \f$ \det G_{ij} - \det g_{ij} = 0 \f$.
747  const DenseMatrix<double>& G,
748  DenseMatrix<double>& sigma_dev,
749  DenseMatrix<double>& G_contra,
750  double& Gdet);
751 
752 
753  /// Calculate the deviatoric part of the contravariant
754  /// 2nd Piola Kirchoff stress tensor. Also return the contravariant
755  /// deformed metric tensor, the generalised dilatation, \f$ d, \f$ and
756  /// the inverse of the bulk modulus \f$ \kappa\f$. This form is appropriate
757  /// for near-incompressible materials for which
758  /// \f$ \sigma^{ij} = -p G^{ij} + \overline{ \sigma^{ij}} \f$
759  /// where the "pressure" \f$ p \f$ is determined from
760  /// \f$ p / \kappa - d =0 \f$.
762  const DenseMatrix<double>& G,
763  DenseMatrix<double>& sigma_dev,
764  DenseMatrix<double>& Gcontra,
765  double& gen_dil,
766  double& inv_kappa);
767 
768 
769  /// Pure virtual function in which the writer must declare if the
770  /// constitutive equation requires an incompressible formulation
771  /// in which the volume constraint is enforced explicitly.
772  /// Used as a sanity check in PARANOID mode. False.
774  {
775  return false;
776  }
777 
778  private:
779  /// Poisson ratio
780  double* Nu_pt;
781 
782  /// Young's modulus
783  double* E_pt;
784 
785  /// Boolean flag to indicate if storage for elastic modulus
786  /// must be deleted in destructor
788  };
789 
790 
791  /// //////////////////////////////////////////////////////////////////////
792  /// //////////////////////////////////////////////////////////////////////
793  /// //////////////////////////////////////////////////////////////////////
794 
795 
796  //=====================================================================
797  /// A class for constitutive laws derived from strain-energy functions.
798  /// Theory is in Green and Zerna.
799  //=====================================================================
801  {
802  private:
803  /// Pointer to the strain energy function
805 
806  public:
807  /// Constructor takes a pointer to the strain energy function
809  StrainEnergyFunction* const& strain_energy_function_pt)
810  : ConstitutiveLaw(), Strain_energy_function_pt(strain_energy_function_pt)
811  {
812  }
813 
814  /// Calculate the contravariant 2nd Piola Kirchhoff
815  /// stress tensor. Arguments are the
816  /// covariant undeformed and deformed metric tensor and the
817  /// matrix in which to return the stress tensor.
818  /// Uses correct 3D invariants for 2D (plane strain) problems.
820  const DenseMatrix<double>& G,
821  DenseMatrix<double>& sigma);
822 
823 
824  /// Calculate the deviatoric part
825  /// \f$ \overline{ \sigma^{ij}}\f$ of the contravariant
826  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
827  /// Also return the contravariant deformed metric
828  /// tensor and the determinant of the deformed metric tensor.
829  /// This form is appropriate
830  /// for truly-incompressible materials for which
831  /// \f$ \sigma^{ij} = - p G^{ij} +\overline{ \sigma^{ij}} \f$
832  /// where the "pressure" \f$ p \f$ is determined by
833  /// \f$ \det G_{ij} - \det g_{ij} = 0 \f$.
835  const DenseMatrix<double>& G,
836  DenseMatrix<double>& sigma_dev,
837  DenseMatrix<double>& G_contra,
838  double& Gdet);
839 
840 
841  /// Calculate the deviatoric part of the contravariant
842  /// 2nd Piola Kirchoff stress tensor. Also return the contravariant
843  /// deformed metric tensor, the generalised dilatation, \f$ d, \f$ and
844  /// the inverse of the bulk modulus \f$ \kappa\f$. This form is appropriate
845  /// for near-incompressible materials for which
846  /// \f$ \sigma^{ij} = -p G^{ij} + \overline{ \sigma^{ij}} \f$
847  /// where the "pressure" \f$ p \f$ is determined from
848  /// \f$ p / \kappa - d =0 \f$.
850  const DenseMatrix<double>& G,
851  DenseMatrix<double>& sigma_dev,
852  DenseMatrix<double>& Gcontra,
853  double& gen_dil,
854  double& inv_kappa);
855 
856 
857  /// State if the constitutive equation requires an incompressible
858  /// formulation in which the volume constraint is enforced explicitly.
859  /// Used as a sanity check in PARANOID mode. This is determined
860  /// by interrogating the associated strain energy function.
862  {
864  }
865  };
866 
867 } // namespace oomph
868 
869 #endif
cstr elem_len * i
Definition: cfortran.h:603
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void error_checking_in_input(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Check for errors in the input, i.e. check that the dimensions of the arrays are all consistent.
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &gen_dil, double &inv_kappa)
Calculate the deviatoric part of the contravariant 2nd Piola Kirchoff stress tensor....
virtual bool requires_incompressibility_constraint()=0
Pure virtual function in which the user must declare if the constitutive equation requires an incompr...
double calculate_contravariant(const DenseMatrix< double > &Gcov, DenseMatrix< double > &Gcontra)
Calculate a contravariant tensor from a covariant tensor, and return the determinant of the covariant...
virtual void calculate_d_second_piola_kirchhoff_stress_dG(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG, const bool &symmetrize_tensor=true)
Calculate the derivatives of the contravariant 2nd Piola Kirchhoff stress tensor with respect to the ...
void calculate_d_contravariant_dG(const DenseMatrix< double > &Gcov, RankFourTensor< double > &dGcontra_dG, DenseMatrix< double > &d_detG_dG)
Calculate the derivatives of the contravariant tensor and the derivatives of the determinant of the c...
bool is_matrix_square(const DenseMatrix< double > &M)
Test whether a matrix is square.
ConstitutiveLaw()
Empty constructor.
bool are_matrices_of_equal_dimensions(const DenseMatrix< double > &M1, const DenseMatrix< double > &M2)
Test whether two matrices are of equal dimensions.
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)=0
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
virtual ~ConstitutiveLaw()
Empty virtual destructor.
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &G_contra, double &Gdet)
Calculate the deviatoric part of the contravariant 2nd Piola Kirchhoff stress tensor ....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
virtual ~GeneralisedHookean()
Virtual destructor.
bool Must_delete_e
Boolean flag to indicate if storage for elastic modulus must be deleted in destructor.
GeneralisedHookean(double *nu_pt)
The constructor takes the pointers to value of Poisson's ratio . Young's modulus is set to E=1....
void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
double * E_pt
Young's modulus.
bool requires_incompressibility_constraint()
Pure virtual function in which the writer must declare if the constitutive equation requires an incom...
double * Nu_pt
Poisson ratio.
GeneralisedHookean(double *nu_pt, double *e_pt)
The constructor takes the pointers to values of material parameters: Poisson's ratio and Young's modu...
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
GeneralisedMooneyRivlin(double *nu_pt, double *c1_pt)
Constructor takes the pointers to the constitutive parameters: Poisson's ratio, the Mooney-Rivlin par...
bool requires_incompressibility_constraint()
Pure virtual function in which the user must declare if the constitutive equation requires an incompr...
double W(const Vector< double > &I)
Return the strain energy in terms of the strain invariants.
bool Must_delete_e
Boolean flag to indicate if storage for elastic modulus must be deleted in destructor.
double * E_pt
Young's modulus.
GeneralisedMooneyRivlin(double *nu_pt, double *c1_pt, double *e_pt)
Constructor takes the pointers to the constitutive parameters: Poisson's ratio, the Mooney-Rivlin par...
virtual ~GeneralisedMooneyRivlin()
Virtual destructor.
double W(const DenseMatrix< double > &gamma)
Return the strain energy in terms of strain tensor.
void derivatives(Vector< double > &I, Vector< double > &dWdI)
Return the derivatives of the strain energy function with respect to the strain invariants.
double * Nu_pt
Poisson's ratio.
double * C1_pt
Mooney-Rivlin parameter.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
StrainEnergyFunction * Strain_energy_function_pt
Pointer to the strain energy function.
bool requires_incompressibility_constraint()
State if the constitutive equation requires an incompressible formulation in which the volume constra...
IsotropicStrainEnergyFunctionConstitutiveLaw(StrainEnergyFunction *const &strain_energy_function_pt)
Constructor takes a pointer to the strain energy function.
void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void derivatives(Vector< double > &I, Vector< double > &dWdI)
Return the derivatives of the strain energy function with respect to the strain invariants.
double * C2_pt
Pointer to second Mooney Rivlin constant.
virtual ~MooneyRivlin()
Empty Virtual destructor.
double W(const DenseMatrix< double > &gamma)
Return the strain energy in terms of strain tensor.
double W(const Vector< double > &I)
Return the strain energy in terms of the strain invariants.
double * C1_pt
Pointer to first Mooney Rivlin constant.
MooneyRivlin(double *c1_pt, double *c2_pt)
Constructor takes the pointer to the value of the constants.
bool requires_incompressibility_constraint()
Pure virtual function in which the user must declare if the constitutive equation requires an incompr...
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
Base class for strain energy functions to be used in solid mechanics computations.
StrainEnergyFunction()
Constructor takes no arguments.
virtual bool requires_incompressibility_constraint()=0
Pure virtual function in which the user must declare if the constitutive equation requires an incompr...
virtual void derivative(const DenseMatrix< double > &gamma, DenseMatrix< double > &dWdgamma)
Return the derivatives of the strain energy function with respect to the components of the strain ten...
virtual ~StrainEnergyFunction()
Empty virtual destructor.
virtual double W(const DenseMatrix< double > &gamma)
Return the strain energy in terms of the strain tensor.
virtual void derivatives(Vector< double > &I, Vector< double > &dWdI)
Return the derivatives of the strain energy function with respect to the strain invariants....
virtual double W(const Vector< double > &I)
Return the strain energy in terms of the strain invariants.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...