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-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
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
41namespace 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
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
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.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...