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