-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathdocumentation.tex
executable file
·3170 lines (2859 loc) · 242 KB
/
documentation.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[10pt,letterpaper,oneside]{report}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{cancel}
\usepackage{color}
\usepackage{hyperref}
\usepackage{tocloft}
\usepackage{textcomp}
\usepackage{parskip}
\usepackage[margin=1in]{geometry}
\usepackage[nodayofweek]{datetime}
\usepackage[usenames,dvipsnames,svgnames,table]{xcolor}
\usepackage[draft]{fixme}
\usepackage{siunitx}
\sisetup{output-exponent-marker=\ensuremath{\mathrm{E}}}
\usepackage{listings}
\lstset{ %
language={[77]Fortran}, % the language of the code
basicstyle=\small\ttfamily, % the size of the fonts that are used for the code
numbers=left, % where to put the line-numbers
numberstyle=\tiny\color{gray}, % the style that is used for the line-numbers
stepnumber=5, % the step between two line-numbers.
numbersep=5pt, % how far the line-numbers are from the code
tabsize=2, % sets default tabsize
aboveskip=16pt, belowskip = 8pt,
breaklines=true, % sets automatic line breaking
breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
showspaces=false, % show spaces adding particular underscores
showstringspaces=false, % underline spaces within strings
showtabs=false, % show tabs within strings adding particular underscores
frame=single, % adds a frame around the code
rulecolor=\color{black}, %
backgroundcolor=\color{White}, % choose the background color. You must add \usepackage{color}
keywordstyle=\color{Blue},
commentstyle=\color{Green},
stringstyle=\color{Red}
}
\usepackage{tikz}
\usetikzlibrary{calc,trees,positioning,arrows,chains,shapes.geometric,decorations.pathreplacing,decorations.pathmorphing,shapes,matrix,shapes.symbols}
\tikzset{
>=stealth',
chain/.style={fill=White, rectangle, rounded corners, draw=black, very thick, node distance=.5cm, text centered, on chain},
line/.style={draw, thick, <-},
lineR/.style={draw, thick, ->},
dottedline/.style={draw, dotted, ->},
double/.style={draw, thick,<->},
every join/.style={->, thick,shorten >=1pt},
input/.style={rectangle, rounded corners, draw=black, very thick, node distance=.5cm, text centered, fill=Dandelion},
local/.style={rectangle, rounded corners, draw=black, very thick, node distance=.5cm, text centered, fill=YellowGreen}
}
\setlength{\parskip}{6 pt}
\pagestyle{headings}
\renewcommand{\leftmark}
% scalars
\newcommand{\sca}[1]{\mbox{\rm{#1}}{}} % roman
\newcommand{\ltr}[1]{\mbox{\sf{#1}}} % sans serif
\newcommand{\bltr}[1]{\mbox{\sffamily{\bfseries{{#1}}}}} % bold sans serif
% matrix
\newcommand{\mat}[1]{\mbox{\bf #1}{}}
% 1st and 2nd order tensors
\newcommand{\ten}[1]{\mbox{\boldmath $#1$}{}}
% 4th order tensors
\newcommand{\tenf}[1]{\mbox{{\sffamily{\bfseries {#1}}}}}
% script sizes
\newcommand{\scas}[1]{\mbox{{\scriptsize{${\rm{#1}}$}}}{}}
\newcommand{\tens}[1]{\mbox{\boldmath{\scriptsize{$#1$}}}{}}
% TITLE PAGE
\newcommand*{\titleGP}{\begingroup
\thispagestyle{empty}
\centering
\vspace*{8\baselineskip} % White space at the top of the page
\rule{\textwidth}{1.6pt}\vspace*{-\baselineskip}\vspace*{2pt}
\rule{\textwidth}{0.4pt}\\[\baselineskip]
{\LARGE The Hitchhiker's Guide to Abaqus}\\[0.2\baselineskip]
\rule{\textwidth}{0.4pt}\vspace*{-\baselineskip}\vspace{3.2pt}
\rule{\textwidth}{1.6pt}\\[\baselineskip] % Thick horizontal line
\scshape
Information, Instructions, Derivations, and Explanations \\
For Use with Abaqus User Subroutines \\
With a Focus on Growth \\ [\baselineskip]
Written August 2012 at Stanford University \par
Last updated \monthname \, \the\year
\vspace*{2\baselineskip} % Whitespace
Written by \\[\baselineskip]
{\Large Dr. Maria Holland\par}
Department of Aerospace \& Mechanical Engineering \\
University of Notre Dame\par
\vspace*{2\baselineskip}
Available under a\par
\href{http://creativecommons.org/licenses/by-nc-sa/4.0/}{
Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License}
\includegraphics{CC_license.png}
\vfill
\endgroup}
\begin{document}
\titleGP
\newgeometry{left=2in,right=2in}
\chapter*{Preface}
\thispagestyle{empty}
This document was born as a compilation of notes and derivations recorded in \LaTeX during my doctoral research on the mechanics of the developing human brain. Over the years and several versions, I have done a lot of additions, consolidations, and corrections. The focus originally was and still is the UMAT subroutine in Abaqus, but as my own research broadened and I received feedback and input from colleagues, collaborators, and other researchers who found this document useful, I have made some attempts to be more comprehensive. Above all, my goal is to present in writing a coherent explanation of some of Abaqus' more esoteric quirks for the benefit of future researchers.
I would like to acknowledge Ellen Kuhl, Manuel Rausch, Adrian Buganza Tepole, Francisco Sahli Costabal, Pim Oomen, and Mohsen Darayi for the input and feedback that has helped me improve this resource over the years.
Additionally, I acknowledge the funding support of the National Science Foundation (DGE-1147470).
I have done my best to accurately attribute my sources; if you find this document useful, please cite the papers in the bibliography where appropriate. If you find any errors or omissions or have suggestions for ways I could improve this resource, please contact me at \href{mailto:[email protected]}{maria-holl{\fontfamily{ptm}\selectfont @}nd.edu}.
\vspace{0.5in}
\begin{center}
\large{\emph{DON'T PANIC.}}
\end{center}
\newpage
\setcounter{tocdepth}{1}
\tableofcontents
\restoregeometry
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Introduction to Abaqus Subroutines}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
User subroutines allow customization for particular purposes. The available subroutines offer a variety of options, from specifying user-defined loading or initial conditions, to the creation of user-defined elements.
% Abaqus subroutines are written in Fortran, and interface with Abaqus via the line \texttt{INCLUDE ABA\_PARAM.INC}, which must be included in every subroutine as the first statement after the argument list.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Variables}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Variables To Be Defined}
The purpose of user subroutines is to perform some calculation as desired by the user instead of Abaqus' built-in functionality. To that end, various subroutines allow users to define variables, including:
\subsubsection{State Variables}
Many user subroutines allow the user to define \texttt{nstatev} state variables, \texttt{statev}. These SDVs, or solution-dependent state variables, are variables that can be defined to evolve with the solution of an analysis as functions of any other variables in the user subroutine. Examples include the growth variable in a growing material, or stretch in a user-defined element. The responsibility for calculating the evolution of the SDVs lies with the user, while ABAQUS simply stores the variables and passes them along from the previous time step and to the next time step.
\subsubsection{Stress}
An important measure in constitutive models, \texttt{stress} is the “true” (Cauchy) stress at the end of the increment. It is a second-order tensor stored as a \texttt{ntens} x 1 Voigt tensor, taking advantage of symmetries.
\subsubsection{Tangent}
The Jacobian matrix of the constitutive model, \texttt{ddsdde}, is also referred to as the tangent or elasticity tensor. Abaqus uses a tangent based on the Jaumann stress tensor. For more information on the Abaqus tangent and its relation to other tangents, see Section \ref{subsec:tangents}. It must be defined accurately if rapid convergence of the overall Newton scheme is to be achieved. In most cases the accuracy of this definition is the most important factor governing the convergence rate. An incorrect definition of the material Jacobian affects only the convergence rate; the results (if obtained) are unaffected. The tangent is a fourth-order tensor stored as a \texttt{ntens} x \texttt{ntens} array (6x6 for fully 3D elements, 4x4 for axisymmetric elements), taking advantage of symmetries.
\subsubsection{Specific Energy}
The specific elastic strain energy (\texttt{sse}), specific plastic dissipation (\texttt{spd}), and specific creep dissipation (\texttt{scd}) should be updated to new values at the end of the increment. They do not affect the solution but are used for energy output.
\subsubsection{Coupled Thermal Analysis}
In the case of full coupled thermal-stress analysis, it is also necessary to define the mechanical volumetric heat generation per unit time (\texttt{rpl}), along with its variation with respect to strain increments (\texttt{drplde}) and temperature (\texttt{drpldt}). Additionally, the variation of the stress increments with respect to the temperature (\texttt{ddsddt}) is needed for the tangent.
\subsection{Variables Passed to the UMAT}
To aid in the calculations performed by user subroutines, various variables are made available to the user, depending on the subroutine being used. To use them, make sure they are contained in all function calls and declarations and are properly declared at the beginning of the subroutine. Please note that this is not comprehensive.
\subsubsection{Names}
User-specified names, such as \texttt{cmname} in UMAT and \texttt{orname} in ORIENT, are useful for distinguishing between two different versions of the same type of subroutine, as discussed in Subsection \ref{subsec:multiple}.
\subsubsection{Properties}
Mechanical and other properties are defined in a \texttt{props} array of length \texttt{nprops} x 1. These properties are not able to be changed during analysis.
\subsubsection{Time and Temperature}
The progress of a simulation is captured in several variables. The step and increment number are stored in \texttt{kstep} and \texttt{kinc}, respectively, while the actual time (and temperature, for thermal coupled simulations), are stored in \texttt{time} and \texttt{temp}, with their increments \texttt{dtime} and \texttt{dtemp}. Note that \texttt{time} is actually a 2 x 1 array, where \texttt{time(1)} is the \emph{step} time at the beginning of the increment, and \texttt{time(2)} is the \emph{total} time at the beginning of the increment. In the case of a simulation with only one step, these two values will be equal.
\subsubsection{Element and Node Numbers}
The information of the integration point where a given calculation is taking place can be found via \texttt{noel} (element number) and \texttt{npt} (integration point number). Its surrounding context can be found via \texttt{nnodes} (number of adjacent nodes) and \texttt{jnnum} (array containing the numbers of adjacent nodes).
\subsubsection{Dimensions}
To write code that is robust enough to run for 3D elements as well as axisymmetric and plane strain elements, the number of direct stress (\texttt{ndi}) and shear (\texttt{nshr}) components, along with the size of the stress/strain component array (\texttt{ntens} = \texttt{ndi} + \texttt{nshr}) should be used to dimension \texttt{stress} and \texttt{ddsdde}. Similarly, the number of state variables (\texttt{nstatv}) and properties (\texttt{nprops}) are user-defined quantities passed in for dimensioning purposes.
\subsubsection{Coordinates}
The coordinates of the integration point can be found in the array \texttt{coords}, while the coordinates of the adjacent nodes found in \texttt{jnnum} are contained in \texttt{cnodes}. In the case geometric nonlinearity, \texttt{coords} are given in the deformed configuration, while \texttt{cnodes} are always the original coordinates.
\subsubsection{Deformation and Strains}
The total strains at the beginning of the increment and the strain increments are passed in as \texttt{stran} and \texttt{dstran}. In the case of thermal expansion, only the mechanical strain is passed in. Alternately, the deformation gradient at both the beginning of the increment (\texttt{dfgrd0}) and at the end of the increment (\texttt{dfgrd1}) are available. If a local orientation is defined at the material point, the deformation gradient is expressed in that coordinate system (see Section \ref{subsec:orientation} for more information).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Select Constitutive User Subroutines}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The following user subroutines, among others, can be used to modify the constitutive behavior of a material. Special attention is given to those which are able to model growing materials.
\subsection{UHYPER, UANISOHYPER\_STRAIN, and UANISOHYPER\_INV}
The UHYPER subroutines defines the strain energy of hyperelastic materials. UHYPER is used for isotropic materials, while UANISOHYPER\_STRAIN and UANISOHYPER\_INV are used for anisotropic materials in terms of components of the Green strain tensor and scalar invariants, respectively. In addition to the strain energy density function, the user must specify the deviatoric part and several derivatives of the function with respect to the strain invariants.
\emph{Has access to:} temperature, SDVs, strain (components or invariants)
\emph{Defines:} strain energy density function and its derivatives
% \subsection{UEXPAND}
% The UEXPAND subroutine is used to define thermal strains that are dependent on temperature, predefined field variables, or state variables.
% \emph{Has access to:} time, temperature, SDVs
% \emph{Defines:} thermal strains and their variation
\subsection{UMAT}
The UMAT subroutine is used to define the mechanical behavior of a material. It can be used for nearly any constitutive model, including mechanical-thermal coupling. It updates the solution-dependent state variables and provides the stress and material Jacobian matrix according to the constitutive model.
\emph{Has access to:} time, temperature, SDVs, deformation gradient
\emph{Defines:} stress and tangent
% \subsection{UEL}
% The UEL subroutine is the most high-level subroutine available within Abaqus. By defining a user element, it is possible to define nearly everything while still using Abaqus' numerical solver. One drawback of the UEL is that it is not possible to display results within Abaqus CAE.
% \emph{Has access to:} ?
% \emph{Defines:} ?
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Other User Subroutines}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{SDVINI}
State variables are initialized to zero by default. Initialization to another value can be done in the input file for simple cases using keyword \texttt{*INITIAL CONDITIONS}. For more complicated scenarios, initialization can be done in the \texttt{SDVINI} subroutine.
\subsection{ORIENT}
Global orientations are used by default, but local material coordinates can be defined if desired. This can be done in the input file for simple cases using keyword \texttt{*ORIENTATION}, but for more complicated scenarios a local coordinate system can be created at each integration point using an \texttt{ORIENT} subroutine to define a 3 x 3 matrix
\begin{equation}
\ten{T} = \left[ \begin{array}{ccc}
\{ \ten{x}' \} & \{ \ten{y}' \} & \{ \ten{z}' \}
\end{array} \right]
= \left[ \begin{array}{ccc}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2 \\
x_3 & y_3 & z_3
\end{array} \right] \, ,
\end{equation}
where $\ten{x}$, $\ten{y}$, and $\ten{z}$ are the local material directions. Only the first two directions must be defined; Abaqus will orthogonalize the second direction with respect to the first and will determine the third direction from the cross product of the first and second.
If local coordinate systems exist, values will be rotated to the current local coordinate system; see Section \ref{subsec:orientation} for more information.
One drawback to this method is that these directions are not accessible to other post-processing programs like Tecplot or Paraview. An inelegant workaround to this is to define the material coordinate system as a material property or a solution-dependent variable at each integration point.
% \subsection{HETVAL}
% The \texttt{HETVAL} subroutine can be used to define a heat flux due to internal heat generation in a material, via dependence on solution-dependent state variables. The variables to be defined are \texttt{FLUX(1)}, the heat flux, with units of energy per time per volume; and \texttt{FLUX(2)}, the rate of change of heat flux per temperature.
% The solution-dependent state variables can be updated, but in a fully coupled temperature-displacement analysis, \texttt{UMAT} is called before \texttt{HETVAL}.
% \subsection{SIGINI}
% Definite initial stress field
% \subsection{USDFLD}
% Redefine field variables
% \subsection{UFIELD}
% Specify predefined field variables
% \subsection{Matrix Operation Subroutines}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Using User Subroutines}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
User subroutines can be called from the input file or from Abaqus/CAE. As the actions in Abaqus/CAE are later translated into an input file, this is often a useful way of producing an input file that can later be manipulated as desired. Some options, including user initial conditions, are not supported in Abaqus/CAE. In this case, the relevant lines must be added manually through \textbf{Model} \textrightarrow \textbf{Edit Keywords} at the end.
When submitting a job in the command line, append \texttt{user=filename} to the job. The subroutine file must be in the same directory as the input file. The \texttt{.f} extension is optional.
When creating a job in Abaqus/CAE, include user subroutines by browsing for the file in \textbf{General} \textrightarrow \textbf{User subroutine file}.
\subsection{Using Multiple Subroutines}
\label{subsec:multiple}
To use multiple subroutines of different kinds in one analysis, they should all be contained in the same subroutine file, in any order. To use multiple subroutines of the same kind - for example, two \texttt{UMAT} subroutines for two different materials - an additional step must be added. The main \texttt{UMAT} subroutine should act as a directory, testing some variable (\texttt{CMNAME}, for example) for different values to distinguish between the cases. A different subroutine (for instance, \texttt{UMAT\_MAT1} or \texttt{UMAT\_MAT2}) could then be called:
\newpage
\begin{quote} \begin{lstlisting}
c... ---------------------------------------------------------------
subroutine UMAT(cmname,...)
c... ---------------------------------------------------------------
include 'aba_param.inc'
character*80 cmname
if (cmname(1:4) .eq. 'mat1') then
call UMAT_MAT1(...)
else if(cmname(1:4) .eq. 'mat2') then
call UMAT_MAT2(...)
end if
return
end
c... ---------------------------------------------------------------
subroutine UMAT_MAT1(...)
c... ---------------------------------------------------------------
include 'aba_param.inc'
...
\end{lstlisting} \end{quote}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Abaqus User Subroutine Conventions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Voigt Notation}
The Abaqus/Standard convention for ordering stress and strain components is:
\[ \ten{v} = \left[ \begin{array}{c}
\sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\
\tau_{12} \\ \tau_{13} \\ \tau_{23}
\end{array} \right] \]
so that the full matrix equivalent is:
\[ \ten{\sigma} =
\left[ \begin{array}{ccc}
v_1 & v_4 & v_5 \\
v_4 & v_2 & v_6 \\
v_5 & v_6 & v_3
\end{array} \right] \]
Note that this is a different convention from FEAP and Abaqus/Explicit, and it may differ from other programs as well.
Fourth-order tensors with minor symmetries can be written in Voigt notation and are stored as 6x6 arrays, numbered:
\[\left[ \begin{array}{cccccc}
D_{1111} & D_{1122} & D_{1133} & D_{1112} & D_{1113} & D_{1123} \\
D_{2211} & D_{2222} & D_{2233} & D_{2212} & D_{2213} & D_{2223} \\
D_{3311} & D_{3322} & D_{3333} & D_{3312} & D_{3313} & D_{3323} \\
D_{1211} & D_{1222} & D_{1233} & D_{1212} & D_{1213} & D_{1223} \\
D_{1311} & D_{1322} & D_{1333} & D_{1312} & D_{1313} & D_{1323} \\
D_{2311} & D_{2322} & D_{2333} & D_{2312} & D_{2313} & D_{2323}
\end{array} \right] \]
Fourth-order tensors with major symmetries appear symmetric in Voigt notation.
\subsection{Stress and Strain Components}
Because different types of elements have different number of stress and strain components, user subroutines should be written to accomodate all the element types with which they will be used. For instance, the stress and tangent should be defined in a loop from 1 to \texttt{ntens}. If it is necessary to hard-code the definition of the stress or tangent, place elements beyond (4,4) in an extra \texttt{if} loop.
\subsection{Orientation}
\label{subsec:orientation}
If a local orientation is used at the same point as most user subroutines, the stress and strain components will be in the local rotated system. This coordinate system, with basis vectors $\ten{e}'_i$ is related to the original coordinate system, with basis vectors $\ten{e}_i$, by the rotation tensor $\ten{R}$, which is one of the components of the polar decomposition of the deformation gradient, $\ten{F} = \ten{v} \cdot \ten{R}$.
\subsubsection{Basis Transformations}
The relationship between the basis vectors of two coordinate systems, $\ten{e}'_i$ and $\ten{e}_i$, is expressed by a rotation tensor $\ten{R}$,
\begin{align}
\ten{e}_i = \ten{R}^{\scas{T}} \cdot \ten{e}'_i = R_{ij} \ten{e}'_j
\qquad \text{and} \qquad
\ten{e}'_i = \ten{R} \cdot \ten{e}_i = R_{ji} \ten{e}_j
\end{align}
where $R_{ij}$ are the $i,j$ components of tensor $\ten{R}$ in a given basis. (See Holzapfel, 1.5 for more information, where $\ten{R} = \ten{Q}$ in their notation).
% % proof of previous equation
% \begin{align}
% \ten{R} \cdot \ten{e}_i &= \left[ R_{jk} \ten{e}_j \otimes \ten{e}_k \right] \cdot \ten{e}'_i \\
% \ten{R} \cdot \ten{e}_i &= R_{jk} \left[ \ten{e}_j \otimes \ten{e}_k \right] \cdot \ten{e}_i \\
% \ten{R} \cdot \ten{e}_i &= R_{jk} \ten{e}_j \left[ \ten{e}_k \cdot \ten{e}_i \right] \\
% \ten{R} \cdot \ten{e}_i &= R_{jk} \ten{e}_j \delta_{ki} \\
% \ten{R} \cdot \ten{e}_i &= R_{ji} \ten{e}_j
% \end{align}
Thus for a vector $\ten{u}$, its components, $u_i$ and $u'_i$, in the two bases, $\lbrace \ten{e}_i \rbrace$ and $\lbrace \ten{e}'_i \rbrace$, respectively, are related by
\begin{align}
u_i &= \ten{u} \cdot \ten{e}_i = \ten{u} \cdot \ten{R}^{\scas{T}} \cdot \ten{e}'_i = \ten{u} \cdot \left[ R_{ij} \ten{e}'_j \right] = R_{ij} \ten{u} \cdot \ten{e}'_j = R_{ij} u'_j
\\ \mbox{and} \qquad
u'_i &= \ten{u} \cdot \ten{e}'_i = \ten{u} \cdot \ten{R} \cdot \ten{e}_i = \ten{u} \cdot \left[ R_{ji} \ten{e}_j \right] = R_{ji} \ten{u} \cdot \ten{e}_j = R_{ji} u_j
\end{align}
or
\begin{align}
[ \ten{u} ] = [ \ten{R} ] \cdot [ \ten{u}' ] \qquad \text{and} \qquad [ \ten{u}' ] = [ \ten{R} ]^{\scas{T}} \cdot [ \ten{u} ]
\end{align}
but note that these equations describe the components of the \emph{same vector in two different bases}, which is different from the equations $\ten{u} = \ten{R} \cdot \ten{u}'$ and $\ten{u}' = \ten{R}^{\scas{T}} \cdot \ten{u}$, which describe \emph{two different vectors}, $\ten{u}$ and $\ten{u}'$.
Thus, a given vector can be represented in several different ways:
\begin{align}
\ten{u} = u_i \ten{e}_i = u'_i \ten{e}'_i
\end{align}
Similarly, tensors can also be represented in different coordinate systems:
\begin{align}
\ten{A} = A_{ij} \ten{e}_i \otimes \ten{e}_j = A'_{ij} \ten{e}'_i \otimes \ten{e}'_j
\end{align}
The tensorial transformation law is:
\begin{align}
A_{ij} &
= \ten{e}_i \cdot \ten{A} \cdot \ten{e}_j
= \ten{R}^{\scas{T}} \cdot \ten{e}'_i \cdot \ten{A} \cdot \ten{R}^{\scas{T}} \cdot \ten{e}'_j
= R_{ik} \ten{e}'_k \cdot \ten{A} \cdot R_{jl} \ten{e}'_l
= R_{ik} R_{jl} \ten{e}'_k \cdot \ten{A} \cdot \ten{e}'_l
= R_{ik} R_{jl} A'_{kl}
\\ \mbox{and} \qquad
A'_{ij} &
= \ten{e}'_i \cdot \ten{A} \cdot \ten{e}'_j
= \ten{R} \cdot \ten{e}_i \cdot \ten{A} \cdot \ten{R} \cdot \ten{e}_j
= R_{ki} \ten{e}_k \cdot \ten{A} \cdot R_{lj} \ten{e}_l
= R_{ki} R_{lj} \ten{e}_k \cdot \ten{A} \cdot \ten{e}_l
= R_{ki} R_{lj} A_{kl}
\end{align}
or
\begin{align}
[ \ten{A} ] = [ \ten{R} ] [ \ten{A}' ] [ \ten{R} ]^{\scas{T}}
\qquad \text{and} \qquad
[ \ten{A}' ] = [ \ten{R} ]^{\scas{T}} [ \ten{A} ] [ \ten{R} ]
\end{align}
\subsubsection{Rotated Coordinate Systems in Abaqus}
When a local orientation is used with a user subroutine UMAT, the stress and strain components will be in the local rotated system, $\ten{e}'_i$. Vectors passed into the UMAT, however, will generally be defined in reference to the original coordinate system. Therefore, they will have the components of the desired vector ($u_i$), but because they are defined in the rotated coordinate system ($\ten{e}'_i$), they will represent a different vector entirely, $\ten{w} = u_i \ten{e}'_i \neq \ten{u}$. The goal, then, is to calculate the components of the desired vector in the rotated coordinate system ($u'_i$, such that $\ten{u} = u'_i \ten{e}'_i$).
Additionally, vectors have both material (reference, or undeformed) and spatial (current, or deformed) configurations in continuum mechanics. These are generally represented in the same coordinate systems, related by the deformation gradient, $\ten{u} = \ten{F} \cdot \ten{U}$. Thus the complete set of relevant vectors and tensors in Abaqus user subroutines with large deformation and local orientation is:
\begin{itemize}
\item $\ten{U} = U_i \ten{e}_i$, undeformed vector in the original coordinate system
\item $\ten{U} = U'_i \ten{e}'_i$, undeformed vector in the rotated coordinate system \\
\textbf{This is the reference vector to be used in the UMAT}
\item $\ten{W} = U_i \ten{e}'_i$, an entirely different vector with components of undeformed vector in the original coordinate system, in the rotated coordinate system. \textbf{This is the vector that gets passed into the UMAT.}
\item $\ten{u} = u_i \ten{e}_i$, deformed vector in the original coordinate system.
\item $\ten{u} = u'_i \ten{e}'_i$, deformed vector in the rotated coordinate system. \\
\textbf{This is the deformed vector to be used in the UMAT.}
\item $\ten{F} = F_{ij} \ten{e}_i \otimes \ten{e}_j $, deformation gradient in the original coordinate system
\item $\ten{F}' = F'_{ij} \ten{e}'_i \otimes \ten{e}'_j $, deformation gradient in the rotated coordinate system. \\
\textbf{This is the deformation tensor available within the UMAT.}
\item $\ten{R}$, rotation tensor between the original and rotated coordinate systems
\item $\ten{v} = v_{ij} \ten{e}_i \otimes \ten{e}_j $, stretch tensor in the original coordinate system
\item $\ten{v'} = v'_{ij} \ten{e}'_i \otimes \ten{e}'_j $, stretch tensor in the rotated coordinate system
\end{itemize}
where capital letters refer to vectors in the reference configuration, lower-case letters refer to vectors in the current configuration, prime indicates the rotated coordinate system, and $\ten{W}$ is an entirely different vector.
We start with the relation $\ten{u} = \ten{F} \cdot \ten{U}$, and prove that an equation that holds in one coordinate system holds in all other coordinate systems,
\begin{align}
\ten{u} ={}& \ten{F} \cdot \ten{U}
\notag \\
u_i \ten{e}_i ={}& F_{kl} \left[ \ten{e}_k \otimes \ten{e}_l \right] \cdot U_j \ten{e}_j
\notag \\
u_i \left[ R_{ia} \ten{e}'_a \right] ={}& \left[ F_{kl} \left[ R_{kc} \ten{e}'_c \right] \otimes \left[ R_{ld} \ten{e}'_d \right] \right] \cdot U_j \left[ R_{jb} \ten{e}'_b \right]
\notag \\
\left[ R_{im} u'_m \right] \left[ R_{ia} \ten{e}'_a \right] ={}& \left[ \left[ R_{kp} R_{lq} F'_{pq} \right] \left[ R_{kc} \ten{e}'_c \right] \otimes \left[ R_{ld} \ten{e}'_d \right] \right] \cdot \left[ R_{jn} U'_n \right] \left[ R_{jb} \ten{e}'_b \right]
\notag \\
R_{im} R_{ia} u'_m \ten{e}'_a ={}& R_{kp} R_{kc} R_{lq} R_{ld} R_{jn} R_{jb} F'_{pq} U'_n \left[ \ten{e}'_c \otimes \ten{e}'_d \right] \cdot \ten{e}'_b
\\
\delta_{am} u'_m \ten{e}'_a ={}& \delta_{pc} \delta_{qd} \delta_{nb} F'_{pq} U'_n \left[ \ten{e}'_c \otimes \ten{e}'_d \right] \cdot \ten{e}'_b
\notag \\
u'_a \ten{e}'_a ={}& F'_{cd} \left[ \ten{e}'_c \otimes \ten{e}'_d \right] \cdot U'_b \ten{e}'_b
\notag \\
\ten{u} ={}& \ten{F} \cdot \ten{U} \, .
\notag
\end{align}
Thus, we can represent this as either
\begin{align}
u_i \ten{e}_i = F_{kl} \left[ \ten{e}_k \otimes \ten{e}_l \right] \cdot U_j \ten{e}_j
\qquad \text{or} \qquad
u'_a \ten{e}'_a = F'_{cd} \left[ \ten{e}'_c \otimes \ten{e}'_d \right] \cdot U'_b \ten{e}'_b \, .
\end{align}
Starting in the rotated coordinate system,
\begin{align}
\ten{u} ={}& \ten{F} \cdot \ten{U}
\\
u'_i \ten{e}'_i ={}& F'_{kl} \left[ \ten{e}'_k \otimes \ten{e}'_l \right] \cdot U'_j \ten{e}'_j
\notag \\
u'_i \ten{e}'_i ={}& \left[ R_{mk} R_{nl} F_{mn} \right] \left[ \ten{e}'_k \otimes \ten{e}'_l \right] \cdot \left[ R_{aj} U_a \right] \ten{e}'_j
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} R_{nl} R_{aj} F_{mn} U_a \left[ \ten{e}'_k \otimes \ten{e}'_l \right] \cdot \ten{e}'_j
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} R_{nl} R_{aj} v_{mp} R_{pn} U_a \, \ten{e}'_k \left[ \ten{e}'_l \cdot \ten{e}'_j \right]
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} R_{nl} R_{aj} v_{mp} R_{pn} U_a \, \ten{e}'_k \delta_{lj}
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} R_{nj} R_{aj} v_{mp} R_{pn} U_a \, \ten{e}'_k
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} \delta_{na} v_{mp} R_{pn} U_a \, \ten{e}'_k
\notag \\
u'_i \ten{e}'_i ={}& R_{mk} v_{mp} R_{pn} U_n \, \ten{e}'_k
\notag \\
u'_i \ten{e}'_i ={}& v'_{kn} U_n \, \ten{e}'_k
\end{align}
Or, ignoring the bases (as they are not changing), and dealing only with the components,
\begin{align}
[ \ten{u}' ] = [ \ten{F}' ] [ \ten{U}' ] = [ \ten{F}' ] [ \ten{R} ]^{\scas{T}} [ \ten{U} ] = [ \ten{v}' ] [ \ten{R} ] [ \ten{R} ]^{\scas{T}} [ \ten{U} ] = [ \ten{v}' ] [ \ten{U} ]
\end{align}
Essentially, this means that the vector passed in to the UMAT, purportedly the reference vector, has already been rotated to the spatial configuration due to Abaqus' use of the rotated coordinate system. To obtain the deformed vector in the current configuration, it merely needs to be stretched by $\ten{v}'$ \cite{SahliCostabal2017}.
To obtain the reference vector in the rotated coordinate system,
\begin{align}
\ten{U} = U'_i \ten{e}'_i = R_{ji} U_j \ten{e}'_i = U_j R_{ji} \ten{e}'_i
\end{align}
or
\begin{align}
[ \ten{U}' ] = [ \ten{R} ]^{\scas{T}} \cdot [ \ten{U} ] = [ \ten{F}' ]^{-1} \cdot [ \ten{v}' ] \cdot [ \ten{U} ] = [ \ten{F}' ]^{-1} \cdot [ \ten{u}' ]
\end{align}
\subsubsection{Rotating Vectors in UMAT}
The polar decomposition of a tensor can be calculated via the principal value decomposition. Fortran has a built-in subroutine, \texttt{SPRIND}, that returns the principal values and directions of a given tensor :
\begin{center}
\texttt{call SPRIND(tensor, values(3), directions(3,3), S, NDI, NSHR)}
\end{center}
where
\begin{itemize}
\item \texttt{tensor(3,3)} is the tensor to be decomposed;
\item \texttt{values(3)} is a vector for the output of principal values;
\item \texttt{directions(3,3)} is a tensor containing the output of principal directions, where directions(i,j) contains the direction cosines of the principal direction corresponding to the i-th principal value;
\item \texttt{S} is an identifier equal to 1 if the tensor contains stresses and 2 if it contains strains;
\item \texttt{NDI} is the number of direct components (passed into the UMAT); and
\item \texttt{NSHR} is the number of shear components (passed into the UMAT)
\end{itemize}
The left, or spatial, stretch tensor, $\ten{v}$, can then be calculated from the principal values and directions as
\begin{align}
\ten{v} = \sum_{i=1}^3 \lambda_i \ten{n}_i \otimes \ten{n}_i
\end{align}
where $\lambda_i^2 = \mathtt{values(i)}$ and $\ten{n}_i = \mathtt{directions(i,:)}$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Abaqus Simulations}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Abaqus Input Files}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Input files contain the necessary geometric, material, and other information needed to run an analysis in Abaqus.
\subsection{Abaqus Input Files Language}
Abaqus input files have their own syntax, and tend to be very sensitive to spaces, lines, and punctuation. Some important points to be aware of:
\begin{itemize}
\item keywords and options begin with an asterisk * and are often followed by parameters. Some parameters are required, while others are optional. Keywords can be looked up in the Abaqus Keywords Reference Manual in the Abaqus documentation.
\item some keywords are followed by option blocks that contain lists of data (nodal coordinates or connectivity, etc.). Each value must be separated by a comma. If there is only one value in an option block, it should end with a comma as well.
\item ** at the beginning of line indicates a comment
\item lines are limited to 256 characters, or 8 parameters
\item Abaqus does not deal with units, so make sure you're using a \href{https://caeai.com/blog/tips-tricks-consistent-fea-units}{consistent set of units}. Examples include:
\begin{itemize}
\item m, N, kg, s, Pa
\item mm, N, kg, s, MPa
\item ft, lbf, slug, s, lbf/ft$^2$
\end{itemize}
\end{itemize}
\subsection{Writing an Input File}
It is possible to write a simple input file from scratch. It may be easier to start with an existing input file, however, and make necessary changes - such as nodes, elements, materials, and loading.
It is also possible to create an input file by creating a model in Abaqus/CAE to the desired level of completion, creating a job associated with the model, and, in the Job Manager, selecting ``Write Input''. This will create an input file that, when run, recreates the contents of your model. This is convenient because Abaqus/CAE makes building complex geometries and meshing easier, while changing options and defining materials and loading is often more efficient in the input file than navigating through windows and menus.
Another way to create input files efficiently is to use macros and Python scripting in conjunction with Abaqus/CAE. Macros allow actions in Abaqus/CAE to be recorded and written to a file in Python; this file can then be edited and run using the Python interface for Abaqus. To make the resulting file more readable, it is advisable to enable Abaqus' coordinate convention for naming features in the model, by entering this code into the terminal in Abaqus/CAE:
\begin{quote} \begin{lstlisting}
session.journalOptions.setValues(replayGeometry=COORDINATE, recoverGeometry=COORDINATE)
\end{lstlisting} \end{quote}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Debugging in Abaqus}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{center}
\emph{“If debugging is the process of removing software bugs, \\ then programming must be the process of putting them in.”}
\end{center}
\subsection{Abaqus Job Files}
When Abaqus runs a job, it creates several files that contain valuable information about the job. When a job exits with an error, the associated files can likely point you in the right direction as you try to solve the problem.
\begin{itemize}
\item \textbf{.com}: command file created by the Abaqus execution procedure. It lists all the options for the job, so you can verify that you're running what you want to run.
\item \textbf{.dat}: printed output file written by the \texttt{analysis}, \texttt{datacheck}, \texttt{parametercheck}, and \texttt{continue} options. It contains the options processed from the input file; the number of elements, nodes, and degrees of freedom; the computing time and estimations of memory used; and requested variable outputs.
\item \textbf{.log}: log file, which contains start and end times for modules run by the current ABAQUS execution procedure, and will display error messages if the procedure did not run successfully.
\item \textbf{.msg}: message file. This file often contains error message, so check here first usually. It also contains information on time steps, convergence parameters, convergence information for every iteration of every increment of every step, and the total number of increments and computing time.
\item \textbf{.odb}: output database written during analysis and is read by the Visualization module in Abaqus/CAE. It is the only file necessary for viewing results.
\item \textbf{.sta}: status file. Contains summaries of every increment, including the step time and number of attempts and iterations.
\item \textbf{.lck}: lock file, which indicates the analysis is running. When this file disappears, the analysis has ended. It may remain after an interrupted analysis, and you will have to delete it before being able to reattempt the analysis.
% \item \textbf{.par}: Modified version of original parametrized input file showing input parameters and their values.
% \item \textbf{.pes}: Modified version of original parametrized input file showing input free of parameter information (after input parameter evaluation and substitution has been performed
% \item \textbf{.pmg}: Parameter evaluation and substitution message file. It is written when the input file is parametrized
\item \textbf{.prt}: part file. This file is used to store part and assembly information and is created even if the input file does not contain an assembly definition.
\end{itemize}
\subsection{Some Tips}
If the message file says that convergence is unlikely, there is most likely an error with the tangent.
If the job completes, but the deformed configuration doesn't display, there was likely an error somewhere along the way that resulted in values of \texttt{NaN}, or Not a Number.
If there is a floating point error, there is probably a division by zero somewhere. Print any values that are in denominators to make sure they are not zero!
User subroutines can write debut output to the printed output (\texttt{.dat}) file or the message (\texttt{.msg}) file.
% \section{Verifying Subroutines}
% Before using a user-written subroutine, it should be verified on a one-element input file. First run tests - uniaxial stretch, uniaxial with finite rotation, and finite shear - with all displacements prescribed to verify the integration algorithm for stresses and state variables. Then run similar tests with prescribed loads to verify the accuracy of the Jacobian.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Convergence Controls}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Accuracy}
Solution control parameters can be used to control the accuracy of solutions in nonlinear analysis. The default parameters are chosen to optimize accuracy and efficiency for most nonlinear problems, but they can be changed using the \texttt{Controls} keyword and the parameters \texttt{GLOBAL}, \texttt{DISPLACEMENT}, or \texttt{ROTATION}. For instance:
\begin{quote}
\begin{lstlisting}
*CONTROLS, PARAMETERS=FIELD, FIELD=GLOBAL
\end{lstlisting}
$R_n^\alpha, C_n^\alpha, q_0^\alpha, q_u^\alpha, R_P^\alpha, \epsilon^\alpha, C_\epsilon^\alpha, R^\alpha_I, C_f, \epsilon^\alpha_I, \epsilon^\alpha_d$
\end{quote}
where two important convergence criterion are
\begin{itemize}
\item $R_n^\alpha = \num{5E-3}$, the ratio of the largest residual to the corresponding average flux norm
\item $R_P^\alpha = \num{10E-2}$, the alternate residual force criterion for a nonlinear problem
\end{itemize}
\subsection{Time Stepping}
The time step can be specified using direct time step control, or the step size can be allowed to vary within given parameters using automatic time stepping. When automatic time stepping is used, Abaqus increases or decreases the time step based on ease of convergence. This behavior is controlled by a large number of parameters, the default values of which can be changed via the General Solution Controls of the Step module in Abaqus/CAE, or with the \texttt{Controls} keyword in the input file:
\begin{quote}
\begin{lstlisting}
*CONTROLS, PARAMETERS=TIME INCREMENTATION
\end{lstlisting}
$I_0, I_R, I_P, I_C, I_L, I_G, I_S, I_A, I_J, I_T, I^C_S, I^C_J, I^C_A \\
D_f, D_C, D_B, D_A, D_S, D_H, D_D, W_G \\
D_G, D_M, D_M^{\scas{dyn}}, D_M^{\scas{diff}}, D_L, D_E, D_R, D_F$
\end{quote}
where the parameters are defined as
\begin{itemize}
\item $I_0 = 4$. After $I_0$ iterations, the behavior of the largest residuals will be checked; if they fail to decrease over two consecutive iterations the increment will be abandoned. This may have to be increased if initial convergence is nonmonotonic. This must be $> 3$.
\item $I_R = 8$. If convergence has not been achieved after $I_R$ iterations and the logarithmic rate of convergence suggests that too many iterations will be required, the time step is reduced. This should be increased if convergence is nonquadratic.
\item $I_P = 9$. After $I_R$ iterations, the alternate residual criterion is used
\item $I_C = 16$. Without quadratic convergence, $I_C$ is the maximum number of expected iterations allowed for logarthmic convergence.
\item $I_L = 10$. If an increment converges in more than $I_L$ iterations, the next time increment will be reduced.
\item $I_G = 4$. If no more than $I_G$ iterations are required in two consecutive increments, the time increment may be increased.
\item $I_S = 12$. Maximum number of iterations caused by severe discontinuities allowed in an increment.
\item $I_A = 5$. After an increment is reduced $I_A$ times, the simulation will be abandoned.
\item $I_J =6 $. Maximum number of severe discontinuity iterations allowed in two consecutive increments for the time step to be increased.
\item $I_T = 3$. the number of prior increments that must have completed without cutbacks in order to increase the time step
\item $I_S^C = 50$. Maximum number of severe discontinuity iterations in an increment.
\item $I_J^C = 50$. Maximum number of severe discontinuity iterations in two consecutive increments for time increment increase.
\item $I_A^C = 50$. Maximum number of contact augmentations.
\item $D_f = 0.25$. Cutback factor used when an increment appears to be diverging (due to non-decreasing maximum residuals).
\item $D_C = 0.5$. Cutback factor used when the increment is abandoned due to too many equilibrium iterations ($ > I_C$).
\item $D_B = 0.75$. Cutback factor used for the next increment when an increment takes $> I_L$ iterations to converge
\item $D_A = 0.85$. Cutback factor used when the time integration accuracy tolerance is exceeded.
\item $D_S = 0.25 $ Cutback factor used when an increment is abandoned after too many severe discontinuity iterations.
\item $D_H = 0.25$. Cutback factor used in case of problems such as excessive distortion in large-displacement problems.
\item $D_D = 1.5$. Increase factor used for the next increment when two consecutive increments converge in less than $I_G$ iterations.
\item $W_G = 0.75$, ratio of average time integration accuracy measure over $I_T$ increments to the corresponding tolerance for the next time increment to be increased.
\item $D_G = 0.8$. Increase factor for the next time increment when the time integration accuracy measure is less than $W_G$ of the tolerance during $I_T$ consecutive increments.
\item $D_M = 1.5$, maximum allowable increase in time step
\item $D_M^{dyn} = 1.25$, maximum allowable increase in time step for dynamic stress analysis.
\item $D_M^{diff} = 2.0$, maximum allowable increase in time step for diffusion-dominated analysis.
\item $D_L = 0.95$. Minimum ratio of $\Delta t_{proposed}$ to $D_M \Delta t_{current}$ for a linear transient problem.
\item $D_E = 0.1$. Minimum time increment ratio ($\delta t_i / \delta t_{i-1}$ for extrapolation to occur
\item $D_R = 1.0$. Maximum ratio of time increment to stability limit for conditionally stable time integration procedures.
\item $D_F = 0.95$. Fraction of stability limit used as time increment when the proposed time limit exceeds the stability limit.
\end{itemize}
% \texttt{inc=100000} option overrides the default number of iterations allowed by Abaqus; complex jobs will require this.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Constitutive Laws}
\label{chap:const}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Strain and Stress}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Strain}
Strain is a measure of relative motion in a body. There are many such measures of deformation, beginning from the simple concept of stretch, quantified as
\begin{equation}
\lambda = \ten{N} \cdot \ten{C} \cdot \ten{N} \, .
\end{equation}
\subsubsection{Strain in one dimension}
Other strain measures are a \emph{function} of stretch,
\begin{equation}
\epsilon = f \left( \lambda \right) \, .
\end{equation}
By expanding a Taylor series around the unstrained state,
\begin{equation}
\epsilon = f \left( 1 \right) + \left[ \lambda - 1 \right] \frac{\partial f}{\partial \lambda} + \frac{1}{2} \left[ \lambda - 1 \right]^2 \frac{\partial^2 f}{\partial \lambda^2} + \cdots
\end{equation}
and requiring that $f \left( 1 \right) = 0$ (strain is zero in the unstretched state), $\partial f / \partial \lambda = 1$ when $\lambda = 1$ (so that small strain is defined as the change in length per unit length), and $\partial f / \partial \lambda > 0$ for all $\lambda$, several common strain measures can be derived:
\begin{itemize}
\item Nominal (Biot's) strain: $ f \left( \lambda \right) = \lambda - 1$
\item Logarithmic strain: $ f \left( \lambda \right) = \ln \lambda$
\item Green's strain: $ f \left( \lambda \right) = 1 / 2 \left[ \lambda^2 - 1 \right]$
\end{itemize}
\subsubsection{Strain in three dimensions}
The basic concept of strain can be generalized to three dimensions by writing the strain matrix using the principal strains,
\begin{equation}
\ten{ \epsilon } = \epsilon_i \ten{N}_i \, ,
\end{equation}
where $\epsilon$ is again a function of stretch. Certain choices of the strain measure, including Green's strain, make calculation of strain possible directly from the deformation gradient, $\ten{F}$.
\subsection{Deformation Rates}
Deformation rates show how tensor fields change with space. The spatial velocity gradient, for instance, tells how the velocity varies over the spatial domain,
\begin{equation}
\ten{l} (\ten{x}, t)
= \frac{\partial \ten{v} (\ten{x}, t) }{\partial\ten{x} }
= \dot{\ten{F}} \ten{F}^{-1} \, .
\end{equation}
It can be split into symmetric and non-symmetric parts,
\begin{equation}
\ten{l} = \ten{d} + \ten{w} \, .
\end{equation}
The symmetric part, $\ten{d}$, is known as the rate of deformation tensor, while the non-symmetric part, $\ten{w}$ is known as the spin tensor and characterizes the rate of rotation,
\begin{equation}
\ten{d} = \frac{1}{2} \left[ \ten{l} + \ten{l}^{\scas{T}} \right]
\qquad \mbox{and} \qquad
\ten{w} = \frac{1}{2} \left[ \ten{l} - \ten{l}^{\scas{T}} \right] \, .
\end{equation}
The rate of deformation $\ten{d}$ can be shown to be the push forward of the time rate of change of the Green-Lagrange strain tensor $\ten{E}$
\begin{equation}
\dot{\ten{E}} = \ten{F}^{\scas{T}} \cdot \ten{d} \cdot \ten{F} \, .
\end{equation}
\subsection{Stress Rates}
\subsubsection{Truesdell Stress Rates}
Tensorial quantities in solid mechanics, including stress and stress rates, must satisfy the criterion of objectivity, which means the quantities must remain unchanged regardless of the position of the observer. The simplest objective stress rate is the Truesdell stress rate, $\ten{\sigma}^{\circ}$, which is the Piola transformation (push forward) of the time derivative of the second Piola-Kirchhoff stress. Using the expressions
\begin{equation}
\frac{\partial }{\partial t} \ten{F} = \ten{l} \cdot \ten{F} \, ,
\qquad
\frac{\partial }{\partial t} \ten{F}^{-1} = - \ten{F}^{-1} \cdot \ten{l} \, ,
\qquad \mbox{and} \qquad
\frac{\partial }{\partial t} \ten{F}^{-\scas{T}} = - \ten{l}^{\scas{T}} \cdot \ten{F}^{-\scas{T}} \, ,
\end{equation}
along with
\begin{equation}
\dot{J} = J \, \mbox{div} (\ten{v})
\qquad \mbox{and} \qquad
\mbox{div} (\ten{v}) = \mbox{tr} (\nabla \ten{v} ) = \mbox{tr} (\ten{l}) \, ,
\end{equation}
the Truesdell stress rate of the Cauchy stress simplifies to
\begin{align}
\ten{\sigma}^{\circ } ={}& J^{-1} \ten{F} \cdot \left[ \frac{\partial }{\partial t} \ten{S} \right] \cdot \ten{F}^{\scas{T}}
\\
% % begin derivations
% ={}& J^{-1} \ten{F} \cdot \left[ \frac{\partial }{\partial t} \left[ J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} \right] \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& J^{-1} \ten{F} \cdot \left[ \dot{J} \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \frac{\partial }{\partial t} (\ten{F}^{-1}) \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \frac{\partial }{\partial t} (\ten{\sigma}) \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \frac{\partial }{\partial t} (\ten{F}^{-\scas{T}}) \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& J^{-1} \ten{F} \cdot \left[ \dot{J} \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} - J \ten{F}^{-1} \cdot \ten{l} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \dot{\ten{\sigma}} \cdot \ten{F}^{-\scas{T}} - J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{l}^{\scas{T}} \cdot \ten{F}^{-\scas{T}} \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& J^{-1} \dot{J} \ten{\sigma} - \ten{l} \cdot \ten{\sigma} + \dot{\ten{\sigma}} - \ten{\sigma} \cdot \ten{l} ^{\scas{T}}
% \notag \\
% % end derivations
\cdots{}& \notag \\
\ten{\sigma}^{\circ }
={}& \dot{\ten{\sigma}} - \ten{l} \cdot \ten{\sigma} - \ten{\sigma} \cdot \ten{l}^{\scas{T}} + \mbox{tr} (\ten{l}) \ten{\sigma} \, .
\end{align}
The Truesdell rate of the Cauchy stress can also be interpreted as the Lie derivative of the Kirchhoff stress; that is,
\begin{align}
J \ten{\sigma}^{\circ} = \mathcal{L}_{\phi} \left( \ten{\tau} \right) = \ten{ \tau}^{\circ}
% % begin derivations
% ={}& \ten{F} \left[ \frac{\partial }{\partial t} \ten{S} \right] \ten{F}^{\scas{T}}
% \notag \\
% ={}& \ten{F} \cdot \left[ \frac{\partial }{\partial t} \left[ J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} \right] \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& \ten{F} \cdot \left[ \dot{J} \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \frac{\partial }{\partial t} (\ten{F}^{-1}) \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \frac{\partial }{\partial t} (\ten{\sigma}) \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \frac{\partial }{\partial t} (\ten{F}^{-\scas{T}}) \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& \ten{F} \cdot \left[ \dot{J} \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} - J \ten{F}^{-1} \cdot \ten{l} \cdot \ten{\sigma} \cdot \ten{F}^{-\scas{T}} + J \ten{F}^{-1} \cdot \dot{ \ten{\sigma} } \cdot \ten{F}^{-\scas{T}} - J \ten{F}^{-1} \cdot \ten{\sigma} \cdot \ten{l}^{\scas{T}} \cdot \ten{F}^{-\scas{T}} \right] \cdot \ten{F}^{\scas{T}}
% \notag \\
% ={}& \dot{J} \ten{\sigma} - J \ten{l} \cdot \ten{\sigma} + J \dot{ \ten{\sigma} } - J \ten{\sigma} \cdot \ten{l}^{\scas{T}}
% \notag \\
% ={}& \underbrace{ \dot{J} \ten{\sigma} + J \dot{ \ten{\sigma} } }_{\dot{ \ten{\tau} } } - J \ten{l} \cdot \ten{\sigma} - J \ten{\sigma} \cdot \ten{l}^{\scas{T}}
% \notag \\
% % end derivations
= \cdots{}
= \dot{ \ten{\tau} } - \ten{l} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{l}^{\scas{T}} \, .
\end{align}
(Note: Holzapfel defines the Lie derivative of a spatial stress field to be the Oldroyd stress rate, and gives the same definition to the result.)
\subsubsection{Jaumann Stress Rates}
Another important objective stress rate is the \hypertarget{Jaumann}{Jaumann-Zaremba rate}.
\begin{equation}
\ten{A}^{\nabla} = \dot{\ten{A}} - \ten{w} \cdot \ten{A} + \ten{A} \cdot \ten{w}
\end{equation}
It is essentially a simplification of the Truesdell stress, assuming only rigid-body rotation ($\ten{F} = \ten{R}$); that is, the Truesdell and Jaumann rates of a stress are equivalent if $\ten{d} = \ten{0}$.
\begin{align}
\ten{\sigma}^{\nabla} ={}& \dot{\ten{\sigma}} - \ten{w} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{w}
\\
\ten{ \tau}^{\nabla} ={}& \dot{\ten{\tau}} - \ten{w} \cdot \ten{\tau} + \ten{\tau} \cdot \ten{w}
\end{align}
The Jaumann stress is used because it is easy to deal with, and leads to a symmetric tangent.
\subsubsection{Relations Between Stress Rates}
The Jaumann stress is related to the Truesdell stress by:
\begin{align}
\ten{ \tau}^{\circ} ={}& \dot{ \ten{\tau} } - \ten{l} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{l}^{\scas{T}}
\\
={}& \dot{ \ten{\tau} } - \left[ \ten{d} + \ten{w} \right] \cdot \ten{\tau} - \ten{\tau} \cdot \left[ \ten{d} + \ten{w} \right]^{\scas{T}}
\notag \\
={}& \dot{ \ten{\tau} } - \ten{d} \cdot \ten{\tau} - \ten{w} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{d}^{\scas{T}} - \ten{\tau} \cdot \ten{w}^{\scas{T}}
\notag \\
={}& \dot{ \ten{\tau} } - \ten{d} \cdot \ten{\tau} - \ten{w} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{d} + \ten{\tau} \cdot \ten{w}
\notag \\
={}& \underbrace{ \dot{ \ten{\tau} } - \ten{w} \cdot \ten{\tau} + \ten{\tau} \cdot \ten{w} }_{ \ten{ \tau}^{\nabla} } - \ten{d} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{d}
\notag \\
\ten{ \tau}^{\circ} ={}& \ten{ \tau}^{\nabla} - \ten{d} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{d}
\end{align}
The four stress rates introduced here are related by:
\pgfdeclarelayer{background}
\pgfsetlayers{background,main}
\begin{center} \begin{tikzpicture}[start chain=going below]
\node [input, fill=white] (sigT) {$\ten{\sigma}^{\circ} = \mathcal{L}_{\phi} \left( \ten{\sigma} \right) = \dot{\ten{\sigma}} - \ten{l} \cdot \ten{\sigma} - \ten{\sigma} \cdot \ten{l} ^{\scas{T}} + \mbox{tr} (\ten{l}) \ten{\sigma} $ };
\node [input, fill=white, right=3cm of sigT] (sigJ) {$\ten{\sigma}^{\nabla} = \dot{\ten{\sigma}} - \ten{w} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{w}$};
\node [input, fill=white, below=3cm of sigJ] (tauJ) {$\ten{\tau}^{\nabla} = \dot{\ten{\tau}} - \ten{w} \cdot \ten{\tau} + \ten{\tau} \cdot \ten{w} $};
\node [input, fill=white, below=3cm of sigT] (tauT) {$\ten{\tau}^{\circ} = \mathcal{L}_{\phi} \left( \ten{\tau} \right) = \dot{ \ten{\tau} } - \ten{l} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{l}^{\scas{T}} $};
\draw [lineR] (sigT) to node [yshift=0.5cm] {$\ten{\sigma}^{\nabla} = \left. \ten{\sigma}^{\circ} \right|_{\tens{d} = 0} $ } (sigJ);
\draw [lineR] (tauT) to node [yshift=-0.5cm] {$\ten{\tau}^{\nabla} = \left. \ten{\tau}^{\circ} \right|_{\tens{d} = 0} $ } (tauJ);
\draw [lineR] (sigT) to node [xshift=-1.0cm] {$\ten{\tau}^{\circ} = J \ten{\sigma}^{\circ} $ } (tauT);
\end{tikzpicture} \end{center}
For more information on stress rates, see Bonet \& Wood, p152, or Holzapfel, p192.
\subsection{Tangents}
\label{subsec:tangents}
Tangents relate the change in stress (or stress rate) to the change in strain (or strain rate).
\subsubsection{Lagrangian Elasticity Tensor}
The Lagrangian elasticity tensor or tangent arises from the linearization of the stress,
\begin{align}
D \ten{S}_{IJ} ={}& \left. \frac{\partial }{\partial \varepsilon} \left( \ten{S}_{IJ} \Big( \ten{E} (\ten{u} + \varepsilon \ten{\delta u} ) \Big) \right) \right|_{\tens{\varepsilon} = 0}
\\
={}& \frac{ \partial \ten{S}_{IJ} }{ \partial \ten{E}_{KL} } \left. \frac{\partial }{\partial \varepsilon} \ten{E}_{KL} (\ten{u} + \varepsilon \ten{\delta u} ) \right|_{\tens{\varepsilon} = 0}
\notag \\
={}& \frac{ \partial \ten{S}_{IJ} }{ \partial \ten{E}_{KL} } D \ten{E}_{KL} (\ten{u} )
\notag \\
D \ten{S}_{IJ} ={}& \tenf{C}_{IJKL} D \ten{E}_{KL} (\ten{u} ) \, .
\end{align}
Thus we have the Lagrangian elasticity tensor,
\begin{equation}
\tenf{C}_{IJKL} = \frac{ \partial \ten{S}_{IJ} }{ \partial \ten{E}_{KL} } = 2 \frac{ \partial \ten{S}_{IJ} }{ \partial \ten{C}_{KL} } \, ,
\end{equation}
which has both minor symmetries but not necessarily major symmetry.
\subsubsection{Eulerian Elasticity Tensor}
The spatial equivalent is ``intractable'' so instead the rate form is used. (Essentially, the terms $\ten{S}$ and $\ten{E}$ can be linearized in the direction of $\ten{v}$ instead of $\ten{u}$.) This is done using the fact that the time derivative of the second Piola-Kirchhoff stress is the pull back of the Truesdell stress rate and the time derivative of the Green-Lagrange strain tensor is the pull back of the rate of deformation,
\begin{equation}
\dot{\ten{S}}
= \ten{F}^{-1} \cdot \ten{ \tau }^{\circ} \cdot \ten{F}^{-\scas{T}}
= J \ten{F}^{-1} \cdot \ten{ \sigma }^{\circ} \cdot \ten{F}^{-\scas{T}}
\qquad \mbox{and} \qquad
\dot{\ten{E}} = \ten{F}^{\scas{T}} \cdot \ten{d} \cdot \ten{F} \, .
\end{equation}
Thus,
\begin{align}
\dot{\ten{S}}_{IJ} ={}& \tenf{C}_{IJKL} \dot{\ten{E}}_{KL}
\\
\ten{F}^{-1}_{Ii} \ten{ \tau }^{\circ}_{ij} \ten{F}^{-1}_{Jj} ={}& \tenf{C}_{IJKL} \ten{F}_{kK} \ten{d}_{kl} \ten{F}_{lL}
\notag \\
\ten{ \tau }^{\circ}_{ij} ={}& \underbrace{ \tenf{C}_{IJKL} \ten{F}_{iI} \ten{F}_{jJ} \ten{F}_{kK} \ten{F}_{lL} }_{J \tenf{c}_{ijkl}} \ten{d}_{kl}
\notag \\
\ten{ \tau }^{\circ}_{ij} ={}& J \tenf{c}_{ijkl} \ten{d}_{kl}
\\
\ten{ \sigma }^{\circ}_{ij} ={}& \tenf{c}_{ijkl} \ten{d}_{kl} \, ,
\end{align}
and the Eulerian elasticity tensor is defined as
\begin{equation}
\tenf{c}_{ijkl} = \frac{1}{J} \tenf{C}_{IJKL} \ten{F}_{iI} \ten{F}_{jJ} \ten{F}_{kK} \ten{F}_{lL}
\qquad \mbox{or} \qquad
\tenf{c} = \frac{1}{J} \big[ \ten{F} \overline{\otimes} \ten{F} \big] : \tenf{C} : \big[ \ten{F}^{\scas{T}} \overline{\otimes} \ten{F}^{\scas{T}} \big] \, ,
\end{equation}
which similarly has both minor symmetries but not necessarily major symmetry.
\subsubsection{Jaumann Elasticity Tensor}
The objective stress rate used by Abaqus is actually the \hyperlink{jaumann}{Jaumann} stress rate instead of the Truesdell stress rate, given by:
\begin{equation}
{C}^{\scas{Abaqus}} = \frac{1}{J} \frac{\partial \Delta \left( J \ten{\sigma} \right)}{\partial \Delta \ten{\varepsilon}}
\end{equation}
The format of the elasticity tensor can be found using the relationship between the two stress rates,
\begin{align}
\ten{ \tau }^{\circ} ={}& J \tenf{c} : \ten{d}
\\
\ten{\tau}^{\nabla } - \ten{d} \cdot \ten{\tau} - \ten{\tau} \cdot \ten{d} ={}& J \tenf{c} : \ten{d}
\notag \\
\ten{\tau}^{\nabla } ={}& J \tenf{c} : \ten{d} + \ten{d} \cdot \ten{\tau} + \ten{\tau} \cdot \ten{d}
\notag \\
\ten{\tau}^{\nabla } ={}& J \left[ \tenf{c} : \ten{d} + \ten{d} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{d} \right] \, .
\end{align}
We are seeking a tangent that will satisfy the equation
\begin{equation}
\ten{\tau}^{\nabla } = J \tenf{c}^\text{Abaqus} : \ten{d} \, .
\end{equation}
That is,
\begin{align}
J \tenf{c}^\text{Abaqus} : \ten{d} ={}& J \left[ \tenf{c} : \ten{d} + \ten{d} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{d} \right]
\notag \\
\tenf{c}^\text{Abaqus} : \ten{d} ={}& \tenf{c} : \ten{d} + \ten{d} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{d}
\notag \\
\tenf{c}^\text{Abaqus} : \ten{d} ={}& \left[ \tenf{c} + \tenf{c}' \right] : \ten{d} \, ,
\end{align}
where $\tenf{c}'$ is some fourth-order tangent that satisfies
\begin{equation}
\tenf{c}' : \ten{d} = \ten{d} \cdot \ten{\sigma} + \ten{\sigma} \cdot \ten{d} \, .
\end{equation}
Because the double contraction with a symmetric tensor is not unique (see Section \ref{sec:tensorderivatives} for an explanation), there are an infinite number of tensors that could satisfy this. We restrict this further by requiring major and minor symmetries and find the form to be
\begin{align}
\tenf{c}'_{ijkl} ={}& \frac{1}{2} \left[ \delta_{ik}\sigma_{jl} + \delta_{il}\sigma_{jk} + \delta_{jk}\sigma_{il} + \delta_{jl}\sigma_{ik} \right]
\\
\tenf{c}' ={}& \frac{1}{2} \left[ \ten{I} \overline{\otimes} \ten{ \sigma} + \ten{I} \underline{\otimes} \ten{\sigma} + \ten{\sigma} \underline{\otimes} \ten{I} + \ten{\sigma} \overline{\otimes} \ten{I} \right]
\\
\tenf{c}' ={}& \left[ \begin{array}{cccccc}
2\sigma_1 & 0 & 0 & \sigma_4 & \sigma_5 & 0 \\
& 2\sigma_2 & 0 & \sigma_4 & 0 & \sigma_6 \\
& & 2\sigma_3 & 0 & \sigma_5 & \sigma_6 \\
& & & \frac{1}{2} (\sigma_1 + \sigma_2) & \frac{1}{2} \sigma_6 & \frac{1}{2} \sigma_5 \\
& sym. & & & \frac{1}{2} (\sigma_1 + \sigma_3) & \frac{1}{2} \sigma_4 \\
& & & & & \frac{1}{2} (\sigma_2 + \sigma_3)
\end{array} \right], .
\end{align}
This extra term, the geometric tangent, can be reduced to purely kinematic terms. Thus we have the Abaqus tangent \cite{Prot2009,Bower},
\begin{equation}
\tenf{c}^\text{Abaqus} = \tenf{c} + \tenf{c}' \, ,
\end{equation}
which relates the volume-normalized Jaumann rate of the Kirchhoff stress to the symmetric rate of deformation tensor,
\begin{equation}
\frac{1}{J} \ten{\tau}^{\nabla } = \tenf{c}^\text{Abaqus} : \ten{d} \, .
\end{equation}
% This can be confirmed:
% \begin{align}
% \tenf{c}^\text{Abaqus}_{ijkl} ={}& \frac{1}{J} \frac{ \partial \ten{\tau}^{\nabla }_{ij} } { \partial \ten{d}_{kl} }
% \\
% ={}& \frac{1}{J} \frac{ \partial \left( \ten{\tau}^{\circ}_{ij} + \ten{d}_{im} \ten{\tau}_{mj} + \ten{\tau}_{im} \ten{d}_{mj} \right) } { \partial \ten{d}_{kl} }
% \notag \\
% ={}& \frac{1}{J} \frac{ \partial \ten{\tau}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{1}{J} \frac{ \partial \left( \ten{d}_{im} \ten{\tau}_{mj} \right) } { \partial \ten{d}_{kl} }
% + \frac{1}{J} \frac{ \partial \left( \ten{\tau}_{im} \ten{d}_{mj} \right) } { \partial \ten{d}_{kl} }
% \notag \\
% ={}& \frac{1}{J} \frac{ \partial \ten{\tau}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{1}{J} \frac{ \partial \ten{d}_{im} } { \partial \ten{d}_{kl} } \ten{\tau}_{mj}
% + \ten{\tau}_{im} \frac{1}{J} \frac{ \partial \ten{d}_{mj} } { \partial \ten{d}_{kl} }
% \notag \\
% ={}& \frac{ \partial \ten{\sigma}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{ \partial \ten{d}_{im} } { \partial \ten{d}_{kl} } \ten{\sigma}_{mj}
% + \ten{\sigma}_{im} \frac{ \partial \ten{d}_{mj} } { \partial \ten{d}_{kl} }
% \notag \\
% ={}& \frac{ \partial \ten{\sigma}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{1}{2} \left( \delta_{ik} \delta_{ml} + \delta_{il} \delta_{mk} \right) \ten{\sigma}_{mj}
% + \frac{1}{2} \ten{\sigma}_{im} \left[ \delta_{mk} \delta_{jl} + \delta_{ml} \delta_{jk} \right]
% \notag \\
% ={}& \frac{ \partial \ten{\sigma}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{1}{2} \left[ \delta_{ik} \ten{\sigma}_{lj} + \delta_{il} \ten{\sigma}_{kj} \right]
% + \frac{1}{2} \left[ \delta_{jl} \ten{\sigma}_{ik} + \delta_{jk} \ten{\sigma}_{il} \right]
% \notag \\
% ={}& \frac{ \partial \ten{\sigma}^{\circ}_{ij} } { \partial \ten{d}_{kl} }
% + \frac{1}{2} \left[ \delta_{ik} \ten{\sigma}_{lj} + \delta_{il} \ten{\sigma}_{kj} + \delta_{jl} \ten{\sigma}_{ik} + \delta_{jk} \ten{\sigma}_{il} \right]
% \notag \\
% \tenf{c}^\text{Abaqus}_{ijkl} ={}& \tenf{c}_{ijkl} + \frac{1}{2} \left[ \delta_{ik} \ten{\sigma}_{lj} + \delta_{il} \ten{\sigma}_{kj} + \delta_{jl} \ten{\sigma}_{ik} + \delta_{jk} \ten{\sigma}_{il} \right]
% \\
% \tenf{c}^\text{Abaqus} ={}& \tenf{c} + \frac{1}{2} \left[ \ten{I} \overline{\otimes} \ten{\sigma} + \ten{I} \underline{\otimes} \ten{\sigma} + \ten{\sigma} \underline{\otimes} \ten{I} + \ten{\sigma} \overline{\otimes} \ten{I} \right]
% \end{align}
\subsubsection{Mixed Formulation}
Some calculations are performed in the spatial domain (as in Abaqus), some in the reference domain, and some in a `mixed domain' (as in the Kuhl lab Matlab code). The latter, the mixed domain formulation, uses the First Piola-Kirchhoff stress
and the mixed-domain tangent, defined as
\begin{equation}
\ten{P} = \frac{ \partial W_0 }{ \partial \ten{F} }
\qquad \mbox{and} \qquad
\tenf{A} = \frac{ \partial \ten{P} }{ \partial \ten{F} }
\end{equation}
These two quantities are related to their reference counterparts by
\begin{equation}
\ten{P} = \ten{F} \cdot \ten{S}
\qquad \mbox{or} \qquad
\ten{P} = \ten{F}_{iI} \ten{S}_{IJ}
\end{equation}
and
\begin{equation}
\tenf{A} = \left( \ten{F} \overline{\otimes} \ten{I} \right) : \tenf{C}_{IJKL} : \left( \ten{F}^{\scas{T}} \overline{\otimes} \ten{I} \right) + \ten{I} \overline{\otimes} \ten{S}
\qquad \mbox{or} \qquad
\tenf{A} = \ten{F}_{iI} \tenf{C}_{IJKL} \ten{F}_{kK} + \delta_{ik} \ten{S}_{JL}
\end{equation}