-
Notifications
You must be signed in to change notification settings - Fork 28
/
Copy pathcdpm2vumat.f
3190 lines (2833 loc) · 143 KB
/
cdpm2vumat.f
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
c Concrete damage-plasticity model 2 (CDPM2) _ VUMAT for ABAQUS
c
c The CDPM2 was orgiginally developed by the research group of Dr. Peter Grassl (Univ. of Glasgow, UK)
c and has been implemented in LS-DYNA as MAT CDPM (MAT 273).
c Afterwards, it was revised for use in ABAQUS by Seungwook Seok (a PhD student in Civil Eng. at Purdue Univ., USA)
c updated in 2019-07-21
c
c*****************************************************************
c Support page: http://petergrassl.com/Research/DamagePlasticity/CDPMLSDYNA/index.html
c https://github.com/seungwookseok/ABAQUS-version-CDPM2
c
c Key references:
c 1) P. Grassl, D. Xenos, U. Nyström, R. Rempling, K. Gylltoft.: "CDPM2: A damage-plasticity approach to modelling the failure of concrete". International Journal of Solids and Structures. Volume 50, Issue 24, pp. 3805-3816, 2013.
c 2) P. Grassl, U. Nyström, R. Rempling and K. Gylltoft, "A damage-plasticity model for the dynamic failure of concrete", 8th International Conference on Structural Dynamics, Leuven, Belgium, 2011
c 3) P. Grassl and M. Jirasek: "Damage-plastic model for concrete failure". International Journal of Solids and Structures. Vol. 43, pp. 7166-7196, 2006.
c*****************************************************************
c
c-----------------------------------------------------------------
c # of properties (props) = 13
c # of state variables (state) = 49
c-----------------------------------------------------------------
c User-defined material properties are as follows:
c
c props(1) unitflag: Units flag (US Costumary units = 0, SI units = 1)
c props(2) ym: Young’s modulus
c props(3) pr: Poisson’s ratio
c props(4) fc': Concrete uniaxial compressive strength
c props(5) ft: Concrete uniaxial tensile strength
c props(6) gft: Fracture energy
c props(7) gfc: Crushing energy (not used, to be updated)
c props(8) ah: Parameter of the hardening ductility measure
c props(9) bh: Parameter of the hardening ductility measure
c props(10) ch: Parameter of the hardening ductility measure
c props(11) dh: Parameter of the hardening ductility measure
c props(12) hp: Hardening modulus for qh2
c props(13) as: Parameter of the softening ductility measure
c-----------------------------------------------------------------
c The state variables are stored as:
c
c state(*,1) ---------- kappa
c state(*,2) ---------- equivalent strain
c state(*,3) ---------- plastic strain xx (or 11)
c state(*,4) ---------- plastic strain yy (or 22)
c state(*,5) ---------- plastic strain zz (or 33)
c state(*,6) ---------- plastic strain xy (or 12)
c state(*,7) ---------- plastic strain yz (or 23)
c state(*,8) ---------- plastic strain xz (or 13)
c state(*,9) ---------- kappa tension kdt
c state(*,10) --------- kappa tension 1 kdt1
c state(*,11) --------- kappa tension 2 kdt2
c state(*,12) --------- kappa compression kdc
c state(*,13) --------- kappa compression 1 kdc1
c state(*,14) --------- kappa compression 2 kdc2
c state(*,15) --------- damage variable tension omegaT (wt)
c state(*,16) --------- damage variable tension omegaC (wc)
c state(*,17) --------- strain rate factor used for the dynamic formulation based on quasi-static analysis
c state(*,18) --------- alphac is the compression factor given in paper in IJSS by Grassl et al. in equation 46
c state(*,19) --------- equivalent strain tension eqstrT
c state(*,20) --------- equivalent strain compression eqstrC
c state(*,21) --------- total strain along xx (or 11)
c state(*,22) --------- total strain along yy (or 22)
c state(*,23) --------- total strain along zz (or 33)
c state(*,24) --------- total strain along xy (or 12)
c state(*,25) --------- total strain along yz (or 23)
c state(*,26) --------- total strain along xz (or 13)
c state(*,27) --------- equivalent strain (without rate factor influence)
c state(*,28) --------- element deletion flag: "1" active, "0" inactive (deleted element)
c state(*,29:49) ------ to be updated
c-----------------------------------------------------------------
subroutine vumat(
c Read only (unmodifiable)variables -
& nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
& stepTime, totalTime, dt, cmname, coordMp, charLength,
& props, density, strainInc, relSpinInc,
& tempOld, stretchOld, defgradOld, fieldOld,
& stressOld, stateOld, enerInternOld, enerInelasOld,
& tempNew, stretchNew, defgradNew, fieldNew,
c Write only (modifiable) variables -
& stressNew, stateNew, enerInternNew, enerInelasNew )
include 'vaba_param.inc'
dimension props(nprops), density(nblock), coordMp(nblock,*),
& charLength(nblock), strainInc(nblock,ndir+nshr),
& relSpinInc(nblock,nshr), tempOld(nblock),
& stretchOld(nblock,ndir+nshr),
& defgradOld(nblock,ndir+nshr+nshr),
& fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
& stateOld(nblock,nstatev), enerInternOld(nblock),
& enerInelasOld(nblock), tempNew(nblock),
& stretchNew(nblock,ndir+nshr),
& defgradNew(nblock,ndir+nshr+nshr),
& fieldNew(nblock,nfieldv),
& stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
& enerInternNew(nblock), enerInelasNew(nblock)
character*80 cmname
real*8 eps(6);
integer maxnip,lft,llt
c maxnip ---- variable denoting max number of integration points
integer mx,i,j
real*8 ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
common/cdpmc/ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
c ym --------------- Young's modulus
c pr --------------- Poisson's ratio
c ecc --------------- eccentricity parameter used in the plasticity law. Its calibration is described in Jirasek & Bazant (2002)
c qh0 --------------- value of the derivative of the 1st hardening function at kappa=0 as described in eq. (30) of the IJSS paper by Grassl et al.
c ft --------------- tensile strength
c fc --------------- compressive strength
c hp --------------- value of the derivative of the 2nd hardening function as described in eq. (31) of the IJSS paper by Grassl et al.
c ah,bh,ch,dh -------- hardening parameters used in eq. (33) of the IJSS paper by Grassl et al.
c as,bs,df ----------- softening parameters used in eqs. (56) and (50) of the IJSS paper by Grassl et al.
c fc0 --------------- reference compressive strength used in the definition of the rate factor
c type --------------- softening type used in the formulation of the tensile damage variable
c wf --------------- max crack opening displacement used in the tensile damage law
c wf1 --------------- max crack opening displacement used in the bilinear tensile damage law
c efc --------------- compressive strain used as a parameter in the exponential compressive damage law
c ft1 --------------- tensile stress used in the billinear tensile damage law
c strrateflg --------- variable denoting whether impact phenomena are taken into account in the constitutive law
c failflg ------------ flag denoting when an element should be deleted
c m0 --------------- parameter used in the plasticity law calculated in eq.(20) of the IJSS paper by Grassl et al.
c isoflag------------- flag to denote whether isotropic or anisotropic damage law is used
c Variables used for the evaluation of the plasticity algorithm
real*8 sig(6),sigEff(6),oldStrain(6),convStrain(6),
$ deltaTotStrain(6),
$ tempTotStrain(6),elStrain(6),princStress(3),
$ princDir(3,3),totStrain(6),plastStrain(6),
$ sum,tempTheta,
$ strain(6), sigVTrial,rhoTrial,thetaTrial,apexStress,
$ tempkappaP,yieldval
integer subincCounter,rtype,subincFlag,converged,l,k
c totstrain ----------- total strain vector equal to the sum of plastic and elastic strains
c sigVTrial ----------- trial volumetric stress
c rhoTrial ----------- trial deviatoric stress
c thetaTrial ----------- trial Lode angle
c apexStress ----------- apexstress. Used only in the vertex case of the plasticity algorithm
c yieldval - value of the yield function
c subincCounter - counter of subincrementations performed in the plasticity algorithm
c subincFlag - flag denoting whether or not subincrementation is taking place
c =0 no subincrementation is taking place
c =1 subincrementation is taking place
c rtype - return type of the yield surface
c =0 regular return
c =1 return on the tensile apex of the yield surface
c =2 return on the compressive apex of the yield surface
c converged - integer denoting whether plasticity algorithm has converged
c =0 converged
c =1 not converged
c sigEff(6) ---------- effective stress (no damage included)
c oldStrain(6) ---------- old elastic strain
c convStrain(6) ---------- last elastic strain vector for which the plasticity algorithm has converged
c deltaTotStrain(6) ---------- difference between last and new strain vector
c tempTotStrain(6) ---------- temporary elastic strain
c elStrain(6) ---------- elastic strain
c princStress(3) ---------- array containing principal stresses
c princDir(3,3) ---------- matrix containing principal eigenvectors stored columnwise
c plastStrain(6) ---------- plastic strain vector
c sum ---------- variable used for summation
c tempTheta ---------- temporary Lode angle
c strain(6) ---------- strain array
c l,k ---------- counters used in various iterations
c
c Variables used in the damage algorithm
real*8 alpha,omegaT,omegaC,effStressT(6),effStressC(6),
$ rateFactor,cdpm2u_computeRateFactor,strainrate(6),epsilonT,
$ epsilonC,kappaDT,kappaDT1,kappaDT2,kappaDC,kappaDC1,kappaDC2,
$ length,deltaElStrain(6),sigOld(6), damagedGP
c alpha --------- variable used to identify amount of contribution of compressive stresses and is defined in eq. (46) of the IJSS paper by Grassl et al.
c omegaT --------- tensile damage variable
c omegaC --------- compressive damage variable
c effStressT --------- effective tensile part of the effective stress tensor(no damage included) as described in Section 2.1 of the IJSS paper by Grassl et al.
c effStressC --------- effective compressive part of the effective stress tensor(no damage included) as described in Section 2.1 of the IJSS paper by Grassl et al.
c rateFactor ---------- variable used to incorporate impact effects on the constitutive law
c cdpm2u_computeRateFactor --------- function used to calculate the rate factor
c strainrate(6) ---------- rate of the strain tensor used for the calculation of the rateFacto
c epsilonT ---------- tensile equivalent strain
c epsilonC ---------- compressive equivalent strain
c kappaDT ---------- history parameter kappaDT
c kappaDT1 ---------- history parameter kappaDT1
c kappaDT2 ---------- history parameter kappaDT2
c kappaDC ---------- history parameter kappaDC
c kappaDC1 ---------- history parameter kappaDC1
c kappaDC2 ---------- history parameter kappaDC2
c length ---------- characteristic length used for the formulation of the tensile damage function based on the crack-band approach
c deltaElStrain(6) ---------- elastic strains of the previous step
c damagedGP ---------- number of damaged gausspoints within analysed element
c
c-----------------------------------------------------------------
c
c Input variables
real*8 failflag
c type (Damage type)
c = 0.0: Linear softening
c = 1.0: Bi-Linear softening
c = 2.0: Exponential softening
c = 3.0: NO DAMAGE
c strrateflg (Strain rate parameter)
c = 0.0: No rate effects
c = 1.0: Rate effects included
c failflag
c = 0.0: if not ALL gausspoints of the element are damaged
c = 1.0: if ALL gausspoints of the element are damaged
c unitflag
c = 0.0: US Costumary units
c = 1.0: SI units
c
c Additional variables
real*8 fbfcRatio,fb,epsforecc,maxsubinc
real*8 jacobian(4,4),residuals(4),normResiduals,jacobianPrev(4,4),
$ jacobianOld(4,4),residualsOld(4),normResidualsOld,
$ endflag,jacDet,jacDetOld,
$ evecOld(4,4),evalOld(4),itNoOld,roNoOld,
$ evec(4,4),eval(4),itNo,roNo
c
c-----------------------------------------------------------------
c Default parameters:
type = 1.0
strrateflg = 0.0
failflg = 1 ! 0.0: not active, x>0.0: active and element will erode to 1
! in x percent of the integration points. If x=0.60, 60% of
! all integration points must fail before erosion.
maxnip = 1 ! number of the integration points for C3D8R
isoflag = 0.0
bs = 1.0
df = 0.85
printflag = 0.0 ! means nothing!
maxsubinc = 20
endflag = 0
c Material constants
unitflag = props(1)
ym = props(2)
pr = props(3)
fc = props(4)
ft = props(5)
gft = props(6)
c gfc = props(7)
c gfc = 25 ! (N/mm)
ah = props(8)
bh = props(9)
c bh = 0.0025 ! Bh
c bh = -2.29*pepeak + 0.00046 ! (A.10), Grassl and Jirasek, p.25
ch = props(10)
dh = props(11)
hp = props(12)
as = props(13)
c This is the global tolerance. All other tolerances are made relative to this global tolerance
gTol = 1.e-4
damagedGP = 0.
if (ym<0 ) then
isoflag = 1.0
ym = abs(ym)
end if
c Calculated parameters
if (unitflag == 0.0) then
unitconv = 6.89475908677537 ! Convert ksi to MPa
else if (unitflag == 1.0) then
unitconv = 1
end If
fbfcRatio = 1.5*(fc*unitconv)**(-0.075) ! (10), Papanikolaou and Kappos (2007), p.5
fb = fbfcRatio*fc
c fb = 1.16*fc ! Kupfer et al. (1969)
fc0 = ((fc*unitconv)**(1.855)/60.) / unitconv ! (14), Papanikolaou and Kappos (2007), p.8
qh0 = fc0/fc
epsforecc = ft/fb * (fb**2.-fc**2.)/(fc**2.-ft**2.)
ecc = (1.+epsforecc) / (2.-epsforecc)
m0 = 3. * (fc**2.-ft**2.)/(fc*ft) * ecc/(ecc+1.)
e0 = ft/ym ! e0
wf = Gft/(0.225*ft) ! see p.7 in the paper
wf1 = 0.15*wf ! Jirasek and Zimmermann (1998)
ft1 = 0.3*ft ! Jirasek and Zimmermann (1998)
lft = 1 ! number of first material point?
llt = nblock ! number of last material point?
if (stepTime .eq. 0. .and. lanneal .eq. 0) then
call cdpm2u_printallinputvariables()
end if
do i = lft,llt
stateNew(i,28) = stateOld(i,28)
tempkappaP = stateOld(i,1)
length = charLength(i)
efc = 0.0001
jacobianOld(1,1) = stateOld(i,29)
jacobianOld(1,2) = stateOld(i,30)
jacobianOld(1,3) = stateOld(i,31)
jacobianOld(1,4) = stateOld(i,32)
jacobianOld(2,1) = stateOld(i,33)
jacobianOld(2,2) = stateOld(i,34)
jacobianOld(2,3) = stateOld(i,35)
jacobianOld(2,4) = stateOld(i,36)
jacobianOld(3,1) = stateOld(i,37)
jacobianOld(3,2) = stateOld(i,38)
jacobianOld(3,3) = stateOld(i,39)
jacobianOld(3,4) = stateOld(i,40)
jacobianOld(4,1) = stateOld(i,41)
jacobianOld(4,2) = stateOld(i,42)
jacobianOld(4,3) = stateOld(i,43)
jacobianOld(4,4) = stateOld(i,44)
residualsOld(1) = stateOld(i,45)
residualsOld(2) = stateOld(i,46)
residualsOld(3) = stateOld(i,47)
residualsOld(4) = stateOld(i,48)
normResidualsOld = stateOld(i,49)
c Write strain rate vector
eps(1:6) = strainInc(i,1:6)
do l = 1,6
totStrain(l) = eps(l) + stateOld(i,l+20)
strainrate(l) = eps(l)
plastStrain(l) = stateOld(i,l+2)
oldStrain(l) = stateOld(i,l+20)
convStrain(l) = oldStrain(l)
tempTotStrain(l) = totStrain(l)
deltaTotStrain(l) = eps(l)
enddo
subincCounter=0
subincFlag=0
converged=1
do while ( converged .eq. 1 .or. subincFlag .eq. 1 )
do l=1,6
elStrain(l) = tempTotstrain(l) - plastStrain(l)
enddo
call cdpm2u_computeStressesfromStrains(sigEff,elStrain,
$ ym,pr)
call cdpm2u_computeTrialCoordinates(sigEff,sigVTrial,
$ rhoTrial,tempTheta)
thetaTrial=tempTheta
call cdpm2u_computeYieldValue(yieldval,sigVTrial,rhoTrial,
$ thetaTrial,tempKappaP)
apexStress=0.
if (yieldval .gt. 0.) then
call cdpm2u_checkForVertexCase(apexStress,sigVTrial,
$ tempKappaP,rtype)
if (rtype.eq.1 .or. rtype .eq. 2) then
call cdpm2u_performVertexReturn(sigEff,
$ apexStress,tempKappaP,rtype,converged)
end if
if (rtype.eq.0) then
call cdpm2u_performRegularReturn(sigEff,
$ tempKappaP,converged,ym,pr,gTol,jacobian,
$ residuals,normResiduals,endflag)
end if
else
converged = 0
do l=1,6
plastStrain(l) = stateOld(i,l+2)
enddo
goto 925
end if
if ( converged .eq. 1 ) then
subincCounter=subincCounter+1
if ( subincCounter .gt. maxsubinc ) then
call cdpm2u_computeMatrixDeterminant(4,jacobian,jacDet)
call cdpm2u_computeMatrixDeterminant(4,jacobianOld,jacDetOld)
write(*,*) '*** Perform Plasticity return with'
write(*,*) 'subincrementation methodology'
write(*,*) 'No convergence reached !***'
write(*,*) 'jacobianOld(1,:)', jacobianOld(1,:)
write(*,*) 'jacobianOld(2,:)', jacobianOld(2,:)
write(*,*) 'jacobianOld(3,:)', jacobianOld(3,:)
write(*,*) 'jacobianOld(4,:)', jacobianOld(4,:)
write(*,*) 'residualsOld', residualsOld
write(*,*) 'normResidualsOld', normResidualsOld
write(*,*) 'jacobianPrev(1,:)', jacobianPrev(1,:)
write(*,*) 'jacobianPrev(2,:)', jacobianPrev(2,:)
write(*,*) 'jacobianPrev(3,:)', jacobianPrev(3,:)
write(*,*) 'jacobianPrev(4,:)', jacobianPrev(4,:)
write(*,*) 'jacobian(1,:)', jacobian(1,:)
write(*,*) 'jacobian(2,:)', jacobian(2,:)
write(*,*) 'jacobian(3,:)', jacobian(3,:)
write(*,*) 'jacobian(4,:)', jacobian(4,:)
write(*,*) 'jacDet', jacDet
write(*,*) 'jacDetOld', jacDetOld
write(*,*) 'residuals', residuals
write(*,*) 'normResiduals', normResiduals
call jacobi_eigenvalue (4,jacobianOld,1000,evecOld,
$ evalOld,itNoOld,roNoOld)
call jacobi_eigenvalue (4,jacobian,1000,evec,
$ eval,itNo,roNo)
write(*,*) 'evalOld', evalOld
write(*,*) 'eval', eval
endflag = 1
call cdpm2u_performRegularReturn(sigEff,
$ tempKappaP,converged,ym,pr,gTol,jacobian,
$ residuals,normResiduals,endflag)
stop
else if (subincCounter .gt. maxsubinc-1 .and.
$ tempKappaP .lt. 1.0 ) then
tempKappaP = 1.
end if
subIncFlag = 1
do l = 1,6
deltaTotStrain(l) = deltaTotStrain(l)*0.5
tempTotStrain(l) = convStrain(l)+deltaTotStrain(l)
enddo
jacobianPrev = jacobian
else if ( converged .eq. 0 .and.
$ subIncFlag .eq. 0) then
call cdpm2u_computeStrainsfromStresses(sigEff,elStrain,
$ ym,pr)
do l=1,6
plastStrain(l) = totStrain(l)-elStrain(l)
enddo
else if ( converged .eq. 0 .and.
$ subIncFlag .eq. 1) then
c write(*,*) '*** Subincrementation required',subincCounter
call cdpm2u_computeStrainsfromStresses(sigEff,elStrain,
$ ym,pr)
do l = 1,6
plastStrain(l) = tempTotStrain(l)-elStrain(l)
convStrain(l) = tempTotStrain(l)
deltaTotStrain(l) = totStrain(l)-convStrain(l)
tempTotStrain(l) = totStrain(l)
enddo
subincCounter = 0
subincFlag = 0
converged = 1
end if
end do
925 continue
if (type .eq. 3.0) then
omegaT = 0.0
omegaC = 0.0
epsilonT = 0.0
kappaDT = 0.0
kappaDT1 = 0.0
kappaDT2 = 0.0
kappaDC = 0.0
kappaDC1 = 0.0
kappaDC2 = 0.0
omegaT = 0.0
omegaC = 0.0
rateFactor = 0.0
alpha = 0.0
epsilonT = 0.0
epsilonC = 0.0
do l = 1,6
sig(l) = sigEff(l)
enddo
goto 152
end if
c Initialize parameters used in the damage algorithm
rateFactor = stateOld(i,17)
epsilonT = stateOld(i,19)
epsilonC = stateOld(i,20)
kappaDT = stateOld(i,9)
kappaDT1 = stateOld(i,10)
kappaDT2 = stateOld(i,11)
kappaDC = stateOld(i,12)
kappaDC1 = stateOld(i,13)
kappaDC2 = stateOld(i,14)
omegaT = stateOld(i,15)
omegaC = stateOld(i,16)
alpha = stateOld(i,18)
call cdpm2u_computeAlpha(effStressT,effStressC,sigEff,alpha)
sum = 0.
do l = 1,6
c Compute norm of increment of plastic strains
sum = sum + (plastStrain(l)-stateOld(i,l+2))**2.
deltaElStrain(l) = (stateOld(i,l+20)-stateOld(i,l+2))
enddo
call cdpm2u_computeStressesfromStrains(sigOld,deltaElStrain,
$ ym,pr)
sum = sqrt(sum)
call cdpm2u_computeDamage(omegaC,omegaT,strainrate,
$ rateFactor,alpha,epsilonT,epsilonC,kappaDT,kappaDT1,
$ kappaDT2,kappaDC,kappaDC1,kappaDC2,sigEff,sum,
$ tempKappaP,length,sigOld,
$ stateOld(i,18),stateOld(i,27),stateNew(i,27))
do l = 1,6
if (isoflag .eq. 1.0) then
sig(l) = (1.-omegaT)*sigEff(l)
else
sig(l) = (1.-omegaT)*effStressT(l)
$ +(1.-omegaC)*effStressC(l)
end if
end do
152 continue
c Write the history variable at the end of the routine
stateNew(i,1) = tempkappaP
stateNew(i,2) = epsilonT
stateNew(i,3) = plastStrain(1)
stateNew(i,4) = plastStrain(2)
stateNew(i,5) = plastStrain(3)
stateNew(i,6) = plastStrain(4)
stateNew(i,7) = plastStrain(5)
stateNew(i,8) = plastStrain(6)
stateNew(i,9) = kappaDT
stateNew(i,10) = kappaDT1
stateNew(i,11) = kappaDT2
stateNew(i,12) = kappaDC
stateNew(i,13) = kappaDC1
stateNew(i,14) = kappaDC2
stateNew(i,15) = omegaT
stateNew(i,16) = omegaC
stateNew(i,17) = rateFactor
stateNew(i,18) = alpha
stateNew(i,19) = epsilonT
stateNew(i,20) = epsilonC
stateNew(i,21) = totStrain(1)
stateNew(i,22) = totStrain(2)
stateNew(i,23) = totStrain(3)
stateNew(i,24) = totStrain(4)
stateNew(i,25) = totStrain(5)
stateNew(i,26) = totStrain(6)
stateNew(i,29) = jacobian(1,1)
stateNew(i,30) = jacobian(1,2)
stateNew(i,31) = jacobian(1,3)
stateNew(i,32) = jacobian(1,4)
stateNew(i,33) = jacobian(2,1)
stateNew(i,34) = jacobian(2,2)
stateNew(i,35) = jacobian(2,3)
stateNew(i,36) = jacobian(2,4)
stateNew(i,37) = jacobian(3,1)
stateNew(i,38) = jacobian(3,2)
stateNew(i,39) = jacobian(3,3)
stateNew(i,40) = jacobian(3,4)
stateNew(i,41) = jacobian(4,1)
stateNew(i,42) = jacobian(4,2)
stateNew(i,43) = jacobian(4,3)
stateNew(i,44) = jacobian(4,4)
stateNew(i,45) = residuals(1)
stateNew(i,46) = residuals(2)
stateNew(i,47) = residuals(3)
stateNew(i,48) = residuals(4)
stateNew(i,49) = normResiduals
c Element deletion
c if ( stateNew(i,15).gt. 0.9995 .and. stateNew(i,16).gt. 0.9995 ) then
c damagedGP = damagedGP+1.0
c end if
c if (failflg .gt. 0.0) then
c damagedGP = damagedGP/FLOAT(maxnip) ! e.g., maxnip = 1 for C3D8R
c if (damagedGP .ge. failflg) then
c stateNew(i,28) = 0
c end if
c end if
c Write sig components
stressNew(i,1:6) = sig(1:6)
end do
return
end
subroutine cdpm2u_printallinputvariables()
c Subroutine to check if all model parameters have values. If a parameter does not have any values a default value is provided.
real*8 ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
common/cdpmc/ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
c Write output
c Print all input variables
c
write(*,*) '*****************************************************'
write(*,*) '*****************************************************'
write(*,*) ' ABAQUS is using the usermaterial CDPM2 '
write(*,*) ' ------------------------------------------------- '
write(*,*) ' '
write(*,*) ' Material parameters: '
write(*,2) ' YM (Youngs modulus)............... = ',ym
write(*,2) ' PR (Poissons ratio)............... = ',pr
write(*,2) ' ECC (Eccentricity)................ = ',ecc
write(*,2) ' QH0 (Initial hardening)........... = ',qh0
write(*,2) ' FT (Uniaxial tension strength).... = ',ft
write(*,2) ' FC (Uniaxial compression strength) = ',fc
write(*,2) ' HP (Hardening modulus)............ = ',hp
write(*,2) ' AH (Hardening ductility measure).. = ',ah
write(*,2) ' BH (Hardening ductility measure).. = ',bh
write(*,2) ' CH (Hardening ductility measure).. = ',ch
write(*,2) ' DH (Hardening ductility measure).. = ',dh
write(*,2) ' AS (Damage ductility measure)..... = ',as
write(*,2) ' DF (Dilation constant)............ = ',df
write(*,2) ' FC0 (Initial compression strength) = ',fc0
write(*,2) ' TYPE (Damage type)................ = ',type
write(*,*) ' EQ.0.0: Linear softening '
write(*,*) ' EQ.1.0: Bi-Linear softening '
write(*,*) ' EQ.2.0: Exponential softening '
write(*,*) ' EQ.2.0: No damage '
write(*,2) ' BS (Damage: ductility parameter) = ',bs
write(*,2) ' WF (Damage: disp threshold 0)..... = ',wf
write(*,2) ' WF1 (Damage: disp threshold 1).... = ',wf1
write(*,2) ' FT1 (Damage: stress threshold 1).. = ',ft1
write(*,2) ' EFC (strain threshold in comp).... = ',efc
write(*,2) ' STRPR (Strain rate flag).......... = ',strrateflag
write(*,*) ' EQ.0.0: No rate effects '
write(*,*) ' EQ.1.0: Rate effects included '
write(*,2) ' ISOFLAG (isotropic damage flag).. = ',isoflag
write(*,*) ' EQ.0.0: standard model with two damage params'
write(*,*) ' EQ.1.0: model with one damage param '
write(*,*) '*****************************************************'
write(*,*) ' History variables: '
write(*,*) ' Var #1 = kappa'
write(*,*) ' Var #2 = equivalient strain'
write(*,*) ' Var #3 = plastic strain direction 1'
write(*,*) ' Var #4 = plastic strain dierction 2'
write(*,*) ' Var #5 = plastic strain direction 3'
write(*,*) ' Var #6 = hardening tension kdt'
write(*,*) ' Var #7 = hardening tension kdt1'
write(*,*) ' Var #8 = hardening tension kdt2'
write(*,*) ' Var #9 = hardening compression kdc'
write(*,*) ' Var #10 = hardening compression kdc1'
write(*,*) ' Var #11 = hardening compression kdc2'
write(*,*) ' Var #12 = damage function tension wt'
write(*,*) ' Var #13 = damage function compression wc'
write(*,*) ' Var #14 = (internal book keeping)'
write(*,*) ' Var #15 = compression factor alphac'
write(*,*) ' Var #16 = ratefactor alphar'
write(*,*) ' Var #17 = elastic strain direction 1'
write(*,*) ' Var #18 = elastic strain direction 2'
write(*,*) ' Var #19 = elastic strain direction 3'
write(*,*) ' Var #20 = equivalent strain tension'
write(*,*) ' Var #21 = equivalent strain compression'
write(*,*) ' Var #22-#27 = undamaged stresses'
write(*,*) '*****************************************************'
write(*,*) ' Check that the bulk modulus and shear modulus are '
write(*,*) ' correctly calculated. They MUST be set on the '
write(*,*) ' input card. '
write(*,1) ' With E = ',ym,' and pr = ',pr
write(*,2) ' BLK = E/(3*(1-2pr)) = ',ym/(3.0*(1.-2.*pr))
write(*,2) ' SHR = E/(2*(1+pr)) = ',ym/(2.*(1.+pr))
write(*,*) '*****************************************************'
write(*,*) '*****************************************************'
1 format(1x,A,1pE9.3,A,1pE9.3)
2 format(1x,A,1pE12.5)
return
end
subroutine cdpm2u_computeStrainsfromStresses(stress,strain,ym,pr)
c Subroutine to calculate elastic strains from the stress tensor. Performs operation epsilon = D : sigma
real*8 stress(6),strain(6),ym,pr
c strain(6) -------------- elastic strains
c stress(6) -------------- stress tensor
c pr -------------- Poisson's ratio
c ym -------------- Young's modulus
strain(1)=(stress(1) - pr * stress(2) - pr * stress(3))/ym
strain(2)=(-pr*stress(1) + stress(2) - pr * stress(3))/ym
strain(3)=(-pr*stress(1) - pr * stress(2) + stress(3))/ym
strain(4)=(2. * (1+pr) * stress(4) )/ym
strain(5)=(2. * (1+pr) * stress(5) )/ym
strain(6)=(2. * (1+pr) * stress(6) )/ym
return
end
subroutine cdpm2u_computeStressesfromStrains(stress,strain,ym,pr)
c Subroutine to calculate strains from the elastic strain tensor. Performs operation sigma = C : epsilon
c strain(6) -------------- elastic strains
c stress(6) -------------- stress tensor
c pr -------------- Poisson's ratio
c ym -------------- Young's modulus
real*8 factor, stress(6),strain(6),ym,pr
factor = ym/((1.+pr)*(1.-2.*pr))
stress(1)=factor*((1.-pr)*strain(1) + pr * strain(2) +
$ pr * strain(3))
stress(2)=factor*(pr*strain(1) + (1.-pr) * strain(2) +
$ pr * strain(3))
stress(3)=factor*(pr*strain(1) + pr * strain(2) +
$ (1.-pr)*strain(3))
stress(4)=factor*(((1.-2.*pr)/2.) * strain(4) )
stress(5)=factor*(((1.-2.*pr)/2.) * strain(5) )
stress(6)=factor*(((1.-2.*pr)/2.) * strain(6) )
return
end
c ---------------------------------------------------------VERTEX RETURN FUNCTIONS ---------------------------------------------------
subroutine cdpm2u_checkForVertexCase(apexStress,sigV,tempkappa,
$ rtype)
c Subroutine that check whether the current stress state requires plasticity return to the vertex, at which
c derivative of plastic potential and yield surface are discontinuous.
real*8 ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
common/cdpmc/ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
real*8 apexStress,sigV,tempkappa,qh2,cdpm2u_qh2fun
integer rtype
c sigV ----------- volumetric stress <Input>
c tempKappa ----------- cummulative plastic strain <Input>
c apexStress ----------- sigmaV of the yield surface for the current
c tempkappa and rho=theta=0 <Output>
c rtype ----------- return type of the yield surface <Output>
c =0 regular return
c =1 return on the tensile apex of the yield surface
c =2 return on the compressive apex of the yield surface
c qh2 ----------- variables containing the results of the hardening functions
c cdpm2u_qh2fun ----------- function to calculate the hardening function given in eq. (31) of IJSS paper by P. Grassl et al.
if ( sigV .gt. 0. ) then
rtype = 1
if (tempKappa .lt. 1.) then
apexStress = 0.
else
qh2=cdpm2u_qh2fun(tempKappa,hp)
apexStress=qh2*fc/m0
end if
else if ( sigV .lt. 0. .and. tempKappa .lt. 1.) then
rtype = 2
apexStress = 0.
else
rtype = 0
apexStress=0.
end if
return
end
subroutine cdpm2u_performVertexReturn(stress,apexStress,
$ kappa,rtype,converged)
c Subroutine that performs plasticity return close whenever the stress state needs
c to be returned to the apex. If the stress state is not an actual vertex case
c rtype=0 is returned.
real*8 ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
common/cdpmc/ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
real*8 sigV,rho,apexStress,kappa,stress(6),theta,
$ yieldValue,yieldValueMid,sig2,dSig,sigMid,sigAnswer,
$ ratioPotential,kappa0,tempKappaP, cdpm2u_computeTempKappa,
$ cdpm2u_computeRatioPotential,ratioTrial,yldTol
integer i,j,k,rtype,maxiter,converged
c apexStress ------------------ variable containing the apex to return <Input>
c kappa ------------------ cummulative plastic strain <Input>
c rtype ------------------ return type of the yield surface <Input/Output>
c =0 regular return
c =1 return on the tensile apex of the yield surface
c =2 return on the compressive apex of the yield surface
c converged ------------------ integer denoting whether plasticity algorithm has converged
c =0 converged
c =1 not converged
c kappa0,tempKappa---------------- variables for cummulative plastic strains
c yieldValue,yieldValueMid ------- variables with yield values
c sigV ----------------- volumetric stress
c rho ----------------- deviatoric stress
c sig2,sigMid,sigAnswer,dSig ----- variables containing volumetric stress units
c ratioPotential ----------------- variable containing the ratio of the derivatives of plastic potential with respect to rho and sig multiplied by a parameter to convert strains in stresses
c ratioTrial ----------------- ratio of rho/sigmaV
c yldTol ----------------- tolerance used in bisection method solver
c cdpm2u_computeTempKappa--------- function to calculate tempKappa according to eq.(32) of the IJSS paper by Grassl et al.
c cdpm2u_computeRatioPotential --- function to calculate the ratio of the derivatives of plastic potential with respect to rho and sig multiplied by a parameter to convert strains in stresses
c maxiter ---------------- parameter denoting max number of iterations performed using the bisection method
yldTol=gTol
yieldValue = 0.
yieldValueMid = 0.
sig2 = 0.
kappa0=kappa
tempKappaP=kappa
maxiter=250
call cdpm2u_computeTrialCoordinates(stress,sigV,rho,theta)
sig2 = apexStress
tempKappaP =cdpm2u_computeTempKappa(kappa0, sigV, rho, sigV,ym,pr)
call cdpm2u_computeYieldValue(yieldValue,sigV, 0., 0., tempKappaP)
tempKappaP =
$ cdpm2u_computeTempKappa(kappa0, sigV, rho, sig2,ym,pr)
call cdpm2u_computeYieldValue(yieldValueMid,sig2, 0.,0.,
$ tempKappaP)
if ( yieldValue * yieldValueMid .ge. 0. ) then
converged=1
rtype = 0
goto 501
end if
if ( yieldValue .lt. 0.0 ) then
dSig = sig2 - sigV
sigAnswer = sig2
else
dSig = sigV - sig2
sigAnswer = sig2
end if
do j = 1, maxiter
dSig = 0.5 * dSig
sigMid = sigAnswer + dSig
tempKappaP =cdpm2u_computeTempKappa(kappa0, sigV, rho, sigMid,
$ ym,pr)
call cdpm2u_computeYieldValue(yieldValueMid,sigMid, 0., 0.,
$ tempKappaP)
if ( yieldValueMid .le. 0. ) then
sigAnswer = sigMid
end if
if (abs(yieldValueMid) .lt. yldTol .and.
$ yieldValueMid .le. 0.) then
ratioPotential =
$ cdpm2u_computeRatioPotential(sigAnswer, tempKappaP)
ratioTrial = rho / ( sigV - sigAnswer );
if ( ( ( ( ratioPotential .ge. ratioTrial ) .and.
$ rtype .eq. 1 ) ) .or.
$ ( ( ratioPotential .le. ratioTrial ) .and.
$ rtype .eq. 2 ) ) then
goto 500
else
converged=1
rtype = 0
goto 501
endif
endif
enddo
500 do k = 1, 3
stress(k) = sigAnswer
stress(k+3) = 0.
enddo
kappa=tempKappaP
converged=0
501 continue
return
end
real*8 function cdpm2u_computeTempKappa(kappaInitial,sigV1,rho,
$ sigV2,ym,pr)
c Function to calculate the tempKappa whenever requested from the performVertexReturn function.
c TempKappa is calculated according to eq. (32) of the IJSS paper by P. Grassl et al.
real*8 sigV1,rho,sigV2,kappaInitial,ym,pr,
$ equivalentDeltaPlasticStrain,kM,gM,ducMeas,
$ cdpm2u_computeDucMeas
c sigV1 -------------- volumetric stress in the previous stress state <Input>
c sigV2 -------------- volumetric stress in the current stress state <Input>
c rho -------------- deviatoric stress <Input>
c kappaInitial ------- previous kappaP (cummulative plastic strain) <Input>
c ym -------------- Young's modulus <Input>
c pr -------------- Poisson's ratio <Input>
c kM,gM -------------- bulk and shear moduli
c equivalentDeltaPlasticStrain ----- Increase of the plastic strains
c ducMeas------------- ductility measure in plasticity
c cdpm2u_computeDucMeas------ function to calculate the ductility measure according to eq.(33) of IJSS paper by P. Grassl et al.
kM = ym / ( 3. * ( 1. - 2. * pr ) )
gM = ym / ( 2. * ( 1. + pr ) )
equivalentDeltaPlasticStrain = sqrt( 1. / 9. * (( sigV1 - sigV2 )
$ / kM )** 2. + (rho / ( 2. * gM ))** 2. )
ducMeas = cdpm2u_computeDucMeas(sigV2, 0., 3.141592653589793/3.)
cdpm2u_computeTempKappa=kappaInitial+equivalentDeltaPlasticStrain/
$ ducMeas
return
end
real*8 function cdpm2u_computeRatioPotential(sig ,kappa)
c Function to calculate the ratio of the derivatives of the plastic potential, given in eq.(22) of the IJSS paper by P. Grassl et al. with respect to the deviatoric and volumetric stress respectively.
real*8 ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
common/cdpmc/ym,pr,ecc,qh0,ft,fc,hp,ah,bh,ch,dh,as,df,fc0,type,bs,
1 wf,wf1,efc,ft1,strrateflg,failflg,m0,gTol,isoflag,printflag
real*8 AGParam,BGParam,qh1,qh2,cdpm2u_qh1fun,
$ cdpm2u_qh2fun,R,mQ,kappa,
$ sig,rho,
$ dgdsig,dgdrho,Al,Bl
integer j
c sig ---------------- volumetric stress <Input>
c rho ---------------- deviatoric stress <Input>
c kappa ---------------- cummulative plastic strain kappaP <Input>
c dgdsig,dgdrho --------------- derivatives of the plastic potential
c AGParam,BGParam ------------- components of the plastic potential function
c qh1,qh2 ------------ variables containing the results of the hardening functions
c cdpm2u_qh1fun,cdpm2u_qh2fun - functions to calculate the hardening functions given in eqs. (30) in IJSS paper by Grassl et al.
c Al,Bl ---------------- components of the function given in eq.(23) of the IJSS paper by P. Grassl et al.
c R,mQ ---------------- variables
rho=0.
qh1=cdpm2u_qh1fun(kappa,qh0,hp)
qh2=cdpm2u_qh2fun(kappa,hp)
AGParam = ft * qh2 * 3. / fc + m0 / 2.
BGParam = qh2 / 3. * ( 1. + ft / fc ) /( log(AGParam) +
$ log(df + 1.) - log(2.*df - 1.) - log(3. * qh2 + m0 / 2.) )
R = ( sig - ft / 3. * qh2 ) / fc / BGParam
mQ = AGParam * exp(R)
Bl = sig / fc + rho / ( fc * sqrt(6.) )
Al = ( 1. - qh1 ) * Bl**2. + sqrt(1.5) *rho/fc
dgdsig = 4. * ( 1. - qh1 ) / fc * Al * Bl + qh1**2. * mQ / fc
dgdrho = Al / ( sqrt(6.) * fc ) * ( 4. * (1. - qh1 ) * Bl + 6. ) +
$ m0 * (qh1**2.) / ( sqrt(6.) * fc )
cdpm2u_computeRatioPotential= dgdrho/dgdsig*3.*(1.-2.*pr ) /
$ ( 1. + pr )
return
end
c ------------------------------------ Plasticity Algorithm functions (regular return functions) ---------------------------------
subroutine cdpm2u_performRegularReturn(stress,kappa,converged,ym,
$ pr,yieldTol,jacobian,normalisedResid,normOfResiduals,endflag)
c Subroutine to perform regular plasticity return
real*8 stress(6),tempKappa,ym,pr, resid(4), normalisedResid(4),PI,
$ jacobian(4,4), inverseJac(4,4),increment(4),
$ deltaLambda,normOfResiduals, unknowns(4),
$ trialSig,trialRho,trialTheta,kappa,kappaP,tempKappaP,
$ yieldTol,sig,rho,ddkappadDeltaLambdadInv(2),
$ dgdInv(2),ddgddInv(2,2), dkappadDeltaLambda,
$ dfdkappa, ddgdInvDkappa(2), ddkappadDeltaLambdadKappa,
$ ddkappaddDeltaLambdadInv(2),kM,gM,dfdInv(2),sum,
$ stressPrincipal(3),princDir(3,3),
$ endflag
integer i,j,iterations,totIter,converged,error
c stress(6) ------------------ effective stress components are in order xx,yy,zz,xy,yz,xz <Input/Output>
c kappa ------------------ cummulative plastic strain (kappa) <Input>
c ym ------------------ Young's modulus <Input>
c pr ------------------ Poisson's ratio <Input>
c yieldTol ------------------ tolerance of the N-R solver <Input>
c converged ------------------ integer showing whether solution has converged <Output>
c = 0 solution has converged
c = 1 solution has not converged
c resid(4) ------------------ array containing the residuals of the N-R solver
c normalisedResid(4) -------------- array containing the normalised residuals of the N-R solver
c PI ------------------ parameter defined as the pi=3.14....
c jacobian(4,4) ------------------ matrix containing the jacobian of the problem in order to calculate the solution
c inversJac(4,4) ------------------ matrix containing the inverse matrix of the jacobian of the problem in order to calculate the solution
c increment(4) ------------------ array containing the increments of the unknowns at each N-R iteration
c deltaLambda ------------------ plastic multiplier
c normOfResiduals ----------------- norm the array resid(4)
c unknowns(4) ------------------ array with the unknowns in order sigV,rho,kappa,deltaLambda
c trialSig ------------------ trial volumetric stress (initial guess)
c trialRho ------------------ trial deviatoric stress (initial guess)
c trialTheta ------------------ trial Lode angle (initial guess)
c kappaP ------------------ cummulative plastic strain kappa (initial guess)
c tempKappa ------------------ temporary cummulative plastic strain kappa (each iteration)
c sig ------------------ temporary volumetric stress (each iteration)
c rho ------------------ temporary deviatoric stress (each iteration)
c ddkappadDeltaLambdadInv(2) ------ derivative of the kappa with respect to plastic multiplier and volumetric and deviatoric stress
c dfdInv(2) ------------------ derivative of the yield function with respect to the volumetric and deviatoric stress
c dgdInv(2) ------------------ derivative of the plastic potential with respect to the volumetric and deviatoric stress
c ddgddInv(2,2) ------------------ second derivative of the plastic potential with respect to the volumetric and deviatoric stress
c dkappadDeltaLambda -------------- derivative of kappa with respect to the plastic multiplier
c dfdkappa ------------------ derivative of the yield function with respect to kappa
c ddgdInvDkappa(2) ---------------- derivative of the plastic potential with respect to volumetric and deviatoric stress and to kappa
c ddkappadDeltaLambdadKappa ------- derivative of kappa with respect to the plastic multiplier and kappa
c ddkappaddDeltaLambdadInv(2) ----- derivative of kappa with respect to the plastic multiplier and deviatoric and volumetric stress
c kM ------------------ bulk modulus