v0.13.1
SmallStrainPlasticity.hpp
Go to the documentation of this file.
1/** \file SmallStrainPlasticity.hpp
2 * \ingroup small_strain_plasticity
3 * \brief Operators and data structures for small strain plasticity
4 * \example SmallStrainPlasticity.hpp
5
6 \section problem_formulation Problem formulation
7
8 \subsection theory_small_strain Theory
9
10 \subsubsection introduction Introduction
11
12 In the following document, we define functions,
13 \f$\psi(\varepsilon^e,\alpha)\f$, \f$f(\sigma,\mathbf{q})\f$ and
14 \f$\Psi(\sigma,\mathbf{q})\f$ which are free Helmholtz energy, yield function
15 and flow potential, respectively. Note that free Helmholtz energy is defined
16 by stains and kinematic internal state variables, whereas yield function and
17 flow potential is defined by fluxes (stress tensor and fluxes work conjugated
18 to internal kinematic state variables). In this implantation, we assume that all
19 functions are convex. Also, is assumed that energy potential and flow
20 potential is positive and takes zero value at zero strain and zero internal
21 variable states. It can be shown that for such defined problem, the first law
22 of thermodynamics and non-negative dissipation of energy is automatically
23 satisfied, see \cite de2011computational \cite simo2006computational.
24 Moreover, in the following implementation is assumed that all functions are
25 smooth enough that appropriate derivatives could be calculated.
26
27 In the following implementation, the backward-Euler difference scheme is
28 applied to solve set nonlinear system of equations with unilateral constrains.
29 Admissible stress state is enforced by Lagrange multiplier method and
30 (Khun-Tucker) loading/unloading conditions, see \cite simo2006computational.
31 That lead to well known \em Closet \em point \em projection algorithm.
32
33 Following \cite simo2006computational can be shown that plastic loading
34 unloading conditions can be exclusively characterised in terms of trial
35 state
36 \f[
37 \varepsilon^{e,\textrm{trial}}_{n+1} := \varepsilon_{n+1} - \varepsilon^p_n
38 \f]
39 \f[
40 \alpha^\textrm{trial}_{n+1} := \alpha_n
41 \f]
42 \f[
43 \sigma^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n)
44 \f]
45 \f[
46 \mathbf{q}^\textrm{trial}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^{e,\textrm{trial}}_{n+1},\alpha_n)
47 \f]
48 \f[
49 f^\textrm{trial}_{n+1} := f(\sigma^\textrm{trial}_{n+1},\mathbf{q}_n)
50 \f]
51 then loading/unloading conditions are given
52 \f[
53 \begin{split}
54 f^\textrm{trial}_{n+1} < 0,\quad\textrm{elastic step}\,\to\,\Delta\gamma_{n+1} = 0 \\
55 f^\textrm{trial}_{n+1} \geq 0,\quad\textrm{plastic step}\,\to\,\Delta\gamma_{n+1} > 0
56 \end{split}
57 \f]
58 If first condition is true \f$f^\textrm{trial}_{n+1} < 0\f$, trial stress is solution
59 at integration point. If second condition is true, plastic strains and values
60 of internal state variables are calculated using \em Closet \em point \em
61 projection algorithm.
62
63 \subsection cloaset_point_projection Closet point projection algorithm
64
65 \subsubsection elastic_strain Additive strain decomposition
66
67 \f[
68 \varepsilon^e_{n+1} := \varepsilon_{n+1} - \varepsilon^p_{n+1}
69 \f]
70
71 \subsubsection constitutive_equation Constitutive relations
72
73 \f[
74 \sigma_{n+1} := \frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e_{n+1},\alpha_{n+1})
75 \f]
76
77 \f[
78 \mathbf{q}_{n+1} := \frac{\partial \psi}{\partial \alpha}(\varepsilon^e_{n+1},\alpha_{n+1})
79 \f]
80
81 \subsubsection constrains Constrains
82
83 \f[
84 f_{n+1}:=f(\sigma_{n+1},\mathbf{q}_{n+1})
85 \f]
86
87 \subsubsection evolution Evolution functions
88
89 Flow vector
90 \f[
91 \mathbf{n}_{n+1}:=\frac{\partial \Psi}{\partial \sigma}(\sigma_{n+1},\mathbf{q}_{n+1})
92 \f]
93
94 Hardening law:
95 \f[
96 \mathbf{h}_{n+1}:=-\frac{\partial \Psi}{\partial \mathbf{q}}(\sigma_{n+1},\mathbf{q}_{n+1})
97 \f]
98
99 \subsubsection residual Residual
100
101 \f[
102 \mathbf{R}_{n+1}:=
103 \left\{
104 \begin{array}{c}
105 -\varepsilon^p_{n+1}+\varepsilon^p_n \\
106 -\alpha_{n+1}+\alpha_n \\
107 f_{n+1}
108 \end{array}
109 \right\}
110 +
111 \Delta\gamma_{n+1}
112 \left\{
113 \begin{array}{c}
114 \mathbf{n}_{n+1} \\
115 \mathbf{h}_{n+1} \\
116 0
117 \end{array}
118 \right\}
119 \f]
120
121 \f[
122 \left[
123 \begin{array}{ccc}
124 -1+\gamma
125 \left(
126 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
127 +
128 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\varepsilon^p}
129 \right)
130 &
131 \gamma
132 \left(
133 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
134 +
135 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha}
136 \right)
137 &
138 \mathbf{n}
139 \\
140 \gamma
141 \left(
142 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
143 +
144 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial\alpha\partial\varepsilon^p}
145 \right)
146 &
147 -1+\gamma
148 \left(
149 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
150 +
151 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \alpha}
152 \right)
153 &
154 \mathbf{h}\\
155 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
156 +
157 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha\partial\varepsilon^p}
158 &
159 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
160 +
161 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial\alpha}
162 & 0
163 \end{array}
164 \right]_n
165 \left\{
166 \begin{array}{c}
167 \delta \varepsilon^p \\
168 \delta \alpha \\
169 \delta \gamma
170 \end{array}
171 \right\}_{n+1}
172 =
173 \mathbf{R}_n
174 \f]
175 where
176 \f$\varepsilon^p_{n+1}=\varepsilon^p_n + \delta\varepsilon^p_{n+1}\f$,
177 \f$\alpha^p_{n+1}=\alpha_n + \delta\alpha_{n+1}\f$
178 and
179 \f$\Delta\gamma_{n+1} = \Delta\gamma_n+\delta\gamma_{n+1}\f$.
180 Usually it is assumed free energy split, i.e.
181 \f[
182 \psi(\varepsilon^e,\alpha) = \psi^{\varepsilon}(\varepsilon^e) + \psi^\alpha(\alpha)
183 \f]
184 in that case linearized residual takes the form
185 \f[
186 \left[
187 \begin{array}{ccc}
188 -1+\Delta\gamma
189 \frac{\partial\mathbf{n}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
190 &
191 \Delta\gamma
192 \frac{\partial\mathbf{n}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
193 &
194 \mathbf{n}
195 \\
196 \Delta\gamma
197 \frac{\partial\mathbf{h}}{\partial\sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
198 &
199 -1+\Delta\gamma
200 \frac{\partial\mathbf{h}}{\partial\mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
201 &
202 \mathbf{h}\\
203 \frac{\partial f}{\partial \sigma}\frac{\partial^2 \psi}{\partial \varepsilon^e \partial \varepsilon^p}
204 &
205 \frac{\partial f}{\partial \mathbf{q}}\frac{\partial^2 \psi}{\partial \alpha^2}
206 & 0
207 \end{array}
208 \right]_n
209 \left\{
210 \begin{array}{c}
211 \delta \varepsilon^p \\
212 \delta \alpha \\
213 \delta \gamma
214 \end{array}
215 \right\}_{n+1}
216 =
217 -\mathbf{R}_n
218 \f]
219
220 \f[
221 \chi_{n+1} :=
222 \left\{
223 \begin{array}{c}
224 \delta \varepsilon^p \\
225 \delta \alpha \\
226 \delta \gamma
227 \end{array}
228 \right\}_{n+1}
229 =
230 -
231 \mathbf{A}^{-1}
232 \mathbf{R}_n
233 \f]
234
235 \subsubsection small_strain_algebraic_tangent Algorithmic tangent matrix
236
237 \f[
238 \delta \sigma := \mathbf{C}^e(\delta \varepsilon - \delta \varepsilon^p)
239 \f]
240 where \f$\mathbf{C}^e = \frac{\partial^2 \psi^e}{\partial {\varepsilon^e}^2}\f$ and
241 \f$\delta \varepsilon = \textrm{sym}[ \nabla \delta \mathbf{u}]\f$
242
243 \f[
244 \left\{
245 \begin{array}{c}
246 \mathbf{E}^p \\
247 \mathbf{L}^p \\
248 \mathbf{Y}^p
249 \end{array}
250 \right\}
251 \delta \varepsilon
252 =
253 -
254 \mathbf{A}^{-1}
255 \left\{
256 \begin{array}{c}
257 \Delta\gamma \frac{\partial \mathbf{n}}{\partial \sigma}\mathbf{C}^e \\
258 \Delta\gamma \frac{\partial \mathbf{h}}{\partial \sigma}\mathbf{C}^e \\
259 \frac{\partial f}{\partial \sigma}\mathbf{C}^e
260 \end{array}
261 \right\}
262 \delta \varepsilon
263 \f]
264
265 \f[
266 \delta \varepsilon^p = \mathbf{E}^p \delta \varepsilon
267 \f]
268
269 \f[
270 \delta \alpha = \mathbf{L}^p \delta \varepsilon
271 \f]
272
273 \f[
274 \delta \gamma = \mathbf{Y}^p \delta \varepsilon
275 \f]
276
277 \f[
278 \delta \sigma = \mathbf{C}^e(\mathbf{1} - \mathbf{E}^p)\delta \varepsilon
279 \f]
280
281 \f[
282 \delta \sigma = \mathbf{C}^{ep}\delta \varepsilon
283 \f]
284 where \f$\mathbf{C}^{ep} = \mathbf{C}^e-\mathbf{C}^e\mathbf{E}^p\f$
285
286For more details and easy access to all source files related to this see group
287\ref small_strain_plasticity.
288
289 */
290
291/*
292 * This file is part of MoFEM.
293 * MoFEM is free software: you can redistribute it and/or modify it under
294 * the terms of the GNU Lesser General Public License as published by the
295 * Free Software Foundation, either version 3 of the License, or (at your
296 * option) any later version.
297 *
298 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
299 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
300 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
301 * License for more details.
302 *
303 * You should have received a copy of the GNU Lesser General Public
304 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
305
306#ifndef __SMALL_STRAIN_PLASTICITY_HPP
307#define __SMALL_STRAIN_PLASTICITY_HPP
308
309#ifndef WITH_ADOL_C
310 #error "MoFEM need to be compiled with ADOL-C"
311#endif
312
313PetscErrorCode SmallStrainPlasticityfRes(SNES snes,Vec chi,Vec r,void *ctx);
314PetscErrorCode SmallStrainPlasticityfJac(SNES,Vec,Mat,Mat,void *ctx);
315
316/** \brief Small strain finite element implementation
317* \ingroup small_strain_plasticity
318
319*/
321
322
323
324
326
329
330 /// \brief definition of volume element
332
336 addToRule(0) { }
337
338 /** \brief it is used to calculate nb. of Gauss integration points
339 *
340 * for more details pleas look
341 * Reference:
342 *
343 * Albert Nijenhuis, Herbert Wilf,
344 * Combinatorial Algorithms for Computers and Calculators,
345 * Second Edition,
346 * Academic Press, 1978,
347 * ISBN: 0-12-519260-6,
348 * LC: QA164.N54.
349 *
350 * More details about algorithm
351 * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
352 **/
353 int getRule(int order) {
354 // cerr << "order " << order << endl;
355 return 2*(order-1)+addToRule;
356 }
357
358 };
359
360 MyVolumeFE feRhs; ///< calculate right hand side for tetrahedral elements
361 MyVolumeFE feLhs; //< calculate left hand side for tetrahedral elements
362 MyVolumeFE feUpdate; //< update history variables
363
365 mField(m_field),
366 feRhs(m_field),
367 feLhs(m_field),
368 feUpdate(m_field) {
369 }
370
371 PetscErrorCode createInternalVariablesTag();
372
373 /** \brief common data used by volume elements
374 * \ingroup nonlinear_elastic_elem
375 */
376 struct CommonData {
377 map<string,vector<VectorDouble> > dataAtGaussPts;
378 map<string,vector<MatrixDouble> > gradAtGaussPts;
379 vector<VectorDouble> sTress;
380 vector<MatrixDouble> materialTangentOperator;
381 vector<VectorDouble> plasticStrain;
382 vector<VectorDouble> internalVariables;
383 vector<VectorDouble> internalFluxes;
389 bool bBar;
391 bBar(true) {
392 }
393 };
395
397 vector<VectorDouble > &valuesAtGaussPts;
398 vector<MatrixDouble > &gradientAtGaussPts;
400
401 OpGetDataAtGaussPts(const string field_name,
402 vector<VectorDouble > &values_at_gauss_pts,
403 vector<MatrixDouble > &gardient_at_gauss_pts
404 );
405
406 /** \brief operator calculating field data and gradients
407 */
408 PetscErrorCode doWork(
410 );
411 };
412
414 OpGetCommonDataAtGaussPts(const string field_name,CommonData &common_data);
415 };
416
417
419
421
422 OpUpdate(
423 string field_name,
424 CommonData &common_data
425 );
426
427 PetscErrorCode doWork(
429 );
430
431 };
432
433 struct MakeB {
434 PetscErrorCode makeB(const MatrixAdaptor &diffN,MatrixDouble &B);
435 PetscErrorCode addVolumetricB(const MatrixAdaptor &diffN,MatrixDouble &volB,double alpha);
436 };
437
439
443
444
446 string field_name,
447 CommonData &common_data
448 );
449
450 PetscErrorCode doWork(
452 );
453
454 };
455
457
462
463
465 string field_name,
466 CommonData &common_data
467 );
468
469 PetscErrorCode doWork(
470 int row_side,int col_side,
471 EntityType row_type,EntityType col_type,
474 );
475
476 };
477
478 /** \brief Closest point projection algorithm
479 * \ingroup small_strain_plasticity
480
481 */
483
484
491 double tOl;
492 int gG;
493
496 ublas::symmetric_matrix<double,ublas::lower> C,D;
502 ublas::symmetric_matrix<double,ublas::lower> partial2HSigma;
503 ublas::symmetric_matrix<double,ublas::lower> partial2HFlux;
505
506 vector<int> tAgs;
507
509
511 virtual PetscErrorCode setActiveVariablesW();
512
514 virtual PetscErrorCode setActiveVariablesYH();
515
516 double w,y,h;
518 ublas::vector<adouble> a_sTress,a_internalFluxes;
520
521 virtual PetscErrorCode rEcordW();
522
523 virtual PetscErrorCode rEcordY();
524
525 virtual PetscErrorCode rEcordH();
526
529 virtual PetscErrorCode pLayW();
530 virtual PetscErrorCode pLayW_NoHessian();
531
533 virtual PetscErrorCode pLayY();
534 virtual PetscErrorCode pLayY_NoGradient();
535
538 virtual PetscErrorCode pLayH();
539
540 Mat A;
541 ublas::matrix<double,ublas::column_major> dataA;
543 virtual PetscErrorCode createMatAVecR();
544
545 virtual PetscErrorCode destroyMatAVecR();
546
547 virtual PetscErrorCode evaluatePotentials();
548
549 virtual PetscErrorCode cAlculateR(Vec R);
550
551 virtual PetscErrorCode cAlculateA();
552
553 friend PetscErrorCode SmallStrainPlasticityfRes(SNES,Vec,Vec,void *ctx);
554 friend PetscErrorCode SmallStrainPlasticityfJac(SNES,Vec,Mat,Mat,void *ctx);
555
556 SNES sNes;
557 KSP kSp;
558 PC pC;
559 SNESLineSearch lineSearch;
560 PetscErrorCode snesCreate();
561
562 PetscErrorCode snesDestroy();
563
564 PetscErrorCode solveColasetProjection();
565
568
569 PetscErrorCode consistentTangent();
570
571 virtual PetscErrorCode freeHemholtzEnergy() {
572 PetscFunctionBegin;
573 SETERRQ(
574 PETSC_COMM_SELF,
576 "Not implemented (overload class and implement this function)"
577 );
578 PetscFunctionReturn(0);
579 }
580
581 virtual PetscErrorCode yieldFunction() {
582 PetscFunctionBegin;
583 SETERRQ(
584 PETSC_COMM_SELF,
586 "Not implemented (overload class and implement this function)"
587 );
588 PetscFunctionReturn(0);
589 }
590
591 virtual PetscErrorCode flowPotential() {
592 PetscFunctionBegin;
593 SETERRQ(
594 PETSC_COMM_SELF,
596 "Not implemented (overload class and implement this function)"
597 );
598 PetscFunctionReturn(0);
599 }
600
601 };
602
604
608
609 bool initCp;
611
612
613
615 MoFEM::Interface &m_field,
616 string field_name,
617 CommonData &common_data,
619 bool is_linear = true
620 ):
621 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
622 mField(m_field),
623 commonData(common_data),
624 cP(cp),
625 initCp(false),
626 isLinear(is_linear),
627 hessianWCalculated(false) {
628 }
629
632 PetscErrorCode getTags();
633
634 PetscErrorCode setTagsData(
635 const EntityHandle tet,const int nb_gauss_pts,const int nb_internal_variables
636 );
637
638 PetscErrorCode doWork(
640 );
641
642 };
643
644};
645
646#endif //__SMALL_STRAIN_PLASTICITY_HPP
647
648/***************************************************************************//**
649 * \defgroup small_strain_plasticity Small Strain Plasticity
650 * \ingroup user_modules
651 ******************************************************************************/
652
653/***************************************************************************//**
654 * \defgroup user_modules User modules
655 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
PetscErrorCode SmallStrainPlasticityfJac(SNES, Vec, Mat, Mat, void *ctx)
PetscErrorCode SmallStrainPlasticityfRes(SNES snes, Vec chi, Vec r, void *ctx)
EntitiesFieldData::EntData EntData
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
const FTensor::Tensor2< T, Dim, Dim > Vec
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
const double r
rate factor
constexpr auto field_name
Deprecated interface functions.
@ OPROW
operator doWork function is executed on FE rows
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
ublas::matrix< double, ublas::column_major > dataA
ublas::symmetric_matrix< double, ublas::lower > C
ublas::symmetric_matrix< double, ublas::lower > partial2HSigma
friend PetscErrorCode SmallStrainPlasticityfJac(SNES, Vec, Mat, Mat, void *ctx)
ublas::symmetric_matrix< double, ublas::lower > partial2HFlux
ublas::symmetric_matrix< double, ublas::lower > D
friend PetscErrorCode SmallStrainPlasticityfRes(SNES, Vec, Vec, void *ctx)
common data used by volume elements
map< string, vector< MatrixDouble > > gradAtGaussPts
map< string, vector< VectorDouble > > dataAtGaussPts
vector< MatrixDouble > materialTangentOperator
PetscErrorCode makeB(const MatrixAdaptor &diffN, MatrixDouble &B)
PetscErrorCode addVolumetricB(const MatrixAdaptor &diffN, MatrixDouble &volB, double alpha)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssembleLhs(string field_name, CommonData &common_data)
OpAssembleRhs(string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStress(MoFEM::Interface &m_field, string field_name, CommonData &common_data, ClosestPointProjection &cp, bool is_linear=true)
PetscErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
OpGetCommonDataAtGaussPts(const string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const string field_name, vector< VectorDouble > &values_at_gauss_pts, vector< MatrixDouble > &gardient_at_gauss_pts)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating field data and gradients
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpUpdate(string field_name, CommonData &common_data)
Small strain finite element implementation.
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
PetscErrorCode createInternalVariablesTag()
SmallStrainPlasticity(MoFEM::Interface &m_field)