-
Notifications
You must be signed in to change notification settings - Fork 5
/
UMAT_SSMCWStrainRates_Zambrano.for
2037 lines (1756 loc) · 76.3 KB
/
UMAT_SSMCWStrainRates_Zambrano.for
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
!*******************************************************************************
!*******************************************************************************
!** UMAT prepared by Luis Zambrano-Cruzatty based on Yerro 2015 **
!** email: [email protected], [email protected] **
!*******************************************************************************
!*******************************************************************************
! Copyright 2021 Luis Zambrano-Cruzatty and Alba Yerro
!Permission is hereby granted, free of charge, to any person obtaining a copy of
!this software and associated documentation files (the "Software"), to deal in the
!Software without restriction, including without limitation the rights to use, copy,
!modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
!and to permit persons to whom the Software is furnished to do so, subject to the
!following conditions:
!The above copyright notice and this permission notice shall be included in all copies
!or substantial portions of the Software.
!THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
!INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
!PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
!FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!DEALINGS IN THE SOFTWARE.
Subroutine UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL, DDSDDT, DRPLDE, &
DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED,&
CMNAME, NDI, NSHR, NTENS, & NSTATEV, PROPS, NPROPS, COORDS, &
PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, &
LAYER, KSPT, KSTEP, KINC)
!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:"UMAT" :: UMAT
INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATEV), &
DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), &
STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), &
PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
! Arguments:
! I/O Type
! PROPS I R() : List with model parameters
! DSTRAN I R() : Strain increment
! DTIME I R() : Time increment
! DDSDDE O R(,) : Material stiffness matrix
! STRESS I/O R() : stresses
! STATEV I/O R() : state variables
!
!--- Local variables
!
Dimension :: DE(6,6), Sig(6), dEpsE(6), dEpsP(6), &
EpsE(6), EpsP(6), dEps(6), Eps(6), ERate(6), Erate0(6)
logical :: switch_smooth, switch_yield
! Local variables:
!
! DE : Linear Elastic constitutive matrix
! dSig : Stress increment vector
! Sig : Stress vector
! dEpsE : Elastic strain increment vector
! dEpsP : Plastic strain increment vector
! dEps : Total strain increment vector
! EpsE : Elastic strain vector
! EpsP : Plastic strain vector
! Eps : Total strain vector
! EpsRate : Total strain rate tensor
!
Rad = 45d0 / datan(1d0)
!
! SSMC w Strain rate model parameters
!
! Contents of PROPS(10) MCSS
! 1 : G shear modulus
! 2 : ENU Poisson's ratio
! 3 : cp peak cohesion
! 4 : cr residual cohesion
! 5 : phip peak friction angle
! 6 : phir residual friction angle
! 7 : psip peak dilation angle
! 8 : psir residual dilation angle
! 9 : eta shape factor
! 10 : alpha_G Shear modulus viscosity factor
! 11 : alpha_K Bulk modulus viscosity factor
! 12 : alpha_chi friction and cohesion viscosity
! 13 : alpha_beta dilation angle viscosity
! 14 : RefRate Reference strain rate
! 15 : Switch_smooth Boolean switch for activating strain rate smoothing
! 16 : N_S Degree of smoothening
G_0 = PROPS(1) ! shear modulus
ENU = PROPS(2) ! Poisson's ratio
cp = PROPS(3) ! peak cohesion
cr = PROPS(4) ! residual cohesion
phip = PROPS(5)/Rad ! peak friction angle (rad)
phir = PROPS(6)/Rad ! residual friction angle (rad)
psip = PROPS(7)/Rad ! peak dilation angle (rad)
psir = PROPS(8)/Rad ! residual dilation angle (rad)
eta = PROPS(9) ! shape factor
alpha_G = PROPS(10) ! Shear modulus viscosity factor
alpha_K = PROPS(11) ! Bulk modulus viscosity factor
if ((alpha_K==0.0d0).and.(alpha_G>0.0d0)) then ! in case user uses 0 it will be computed from nu
alpha_K= alpha_G*2*(1+ENU)/(3*(1-2*ENU))
endif
alpha_chi = PROPS(12) ! friction and cohesion viscosity
if ((alpha_chi==0.0d0).and.(alpha_G>0.0d0)) then ! in case user uses 0 it will be computed from nu
alpha_chi= 0.1*alpha_G
endif
alpha_beta = PROPS(13) ! dilation angle viscosity
if ((alpha_beta==0.0d0).and.(alpha_G>0.0d0)) then ! in case user uses 0 it will be computed from nu
alpha_beta= alpha_G
endif
RefERate= PROPS(14) ! reference strain rate
if (RefERate==0.0d0) then
RefERate=2.7e-5
endif
call dbltobool(PROPS(15), switch_smooth) ! switch for activating strain rate smoothening
N_S=PROPS(16) ! Degree of smoothening
G = STATEV(1) ! Previous shear modulus
bK = STATEV(2) ! Previous bulk modulus
c = STATEV(3) ! cohesion
phi = STATEV(4) ! friction angle
psi = STATEV(5) ! dilatancy angle
call dbltobool(STATEV(6),switch_yield) ! Point is plastifying
Do i = 1,NTENS
Erate0(i) = STATEV(6+i) ! Previous strain rate
end do
Do i = 1,NTENS
EpsP(i) = STATEV(12+i) ! Previous plastic strain
end do
N_i=STATEV(19) ! Current number of strain rate sums
SUM_rate=STATEV(20) ! Current sum of strain rates
!_____Error control state parameters__________________________________________________
Error_Euler_max=0.0d0 !
Error_Yield_last=0.0d0 !
Error_Yield_max=0.0d0 !
! !
!_____________________________________________________________________________________!
X=1/DTIME
if (DTIME==0.0d0) then
ERate= 0.0d0 ! Current strain rate
else
ERate= X*DSTRAN ! Current strain rate
end if
bK_0= 2*G_0*(1+ENU)/(3*(1-2*ENU))
IPL=0 !Variable required to check if the material is already in residual cond.
!***********************************************************************************
!Call the modified Euler algorithm
call SSMC_Strain_Rate(NOEL,F1,F2,G_0, bK_0, G, bK ,cp,cr,phip,phir,psip,psir,eta,c,phi,psi, &
stress,EpsP,DSTRAN,dEpsP,Sig, Erate0, ERate, RefERate, DTIME, alpha_G, &
alpha_K, alpha_chi, alpha_beta, switch_smooth, N_S, N_i, SUM_rate, switch_yield ,IPL& !)
,Error_Euler_max, Error_Yield_last, Error_Yield_max)
!************************************************************************************
!************************************************************************************
!Stress and state variables updating
Do i=1,NTENS
STRESS(i) = Sig(i)
End Do
STATEV(1) = G
STATEV(2) = bK
STATEV(3) = c
STATEV(4) = phi
STATEV(5) = psi
STATEV(6)= logic2dbl(switch_yield)
Do i = 1,NTENS
STATEV(6+i) =Erate(i) ! Current strain rate
end do
Do i = 1,NTENS
STATEV(12+i) = EpsP(i)
end do
STATEV(19)=N_i
STATEV(20)=SUM_rate
!_____Error control state parameters__________________________________________________
! Comment if not wanted !
!_____________________________________________________________________________________!
STATEV(21)=Error_Euler_max
STATEV(22)=Error_Yield_last
STATEV(23)=Error_Yield_max
!************************************************************************************
!************************************************************************************
!Tangent stiffness matrix to be returned done by elastic stiffness
F1 = bK+(4*G/3)
F2 = bK-(2*G/3)
DDSDDE = 0.0
DDSDDE(1:3,1:3) = F2
DDSDDE(1,1) = F1
DDSDDE(2,2) = F1
DDSDDE(3,3) = F1
DDSDDE(4,4) = G
DDSDDE(5,5) = G
DDSDDE(6,6) = G
!*************************************************************************************
!End of UMAT
Return
end subroutine UMAT
Subroutine SSMC_Strain_Rate(IntGlo,D1,D2,G_0, K_0, G, K,cp,cr,phip,phir,psip,psir,eta,c,phi,psi, &
Sig0,EpsP,DEps,dEpsP,SigC, Erate0, ERate, RefRate, Dtime, &
alpha_G, alpha_K, alpha_chi, alpha_beta, switch_smooth, N_S, N_i, SUM_rate, &
switch_yield, IPL,Error_Euler_max, Error_Yield_last, Error_Yield_max)
!**********************************************************************
!
! Elasto-plastic constitutive model with strain softening and strain rate effects,
! based on the MOHR-COULOMB criterion (considering modifications of Abbo & Sloan (1995))
! and the CT. model proposed by Yerro (2015)
! Explicit MODIFIED EULER INTEGRATION SCHEME with automatic error control.
! Final correction of the yield surface drift (END OF STEP CORRECTION).
!
!**********************************************************************
implicit none
!Local variables
integer :: i,j,n,it,m, MAXITER
double precision :: F,F0,F2 !Evaluation of the Yield function
double precision :: alpha, alpha1, alpha0 !Elastic Strain proportion
double precision :: SSTOL !Tolerance Relative Error
double precision :: YTOL !Tolerance Relative Error of the yield function evaluation
double precision :: SPTOL !Tolerance Softening parameters
double precision :: LTOL ! Tolerance for unloading path
logical :: ApplyStrainRateUpdates !Check if the strain rate path crosses the reference line
!double precision :: PTOL ! Tolerance for parallelism between plastic and elastic strain increment
double precision :: Rn !Relative error
double precision :: T,DT,T1,beta,DTmin !Sub-stepping parameters
double precision :: c1,phi1,psi1,c2,phi2,psi2
double precision :: cu, phiu, psiu, Gu, Ku
double precision :: ctol,phitol,psitol !c,phi,psi tolerances
double precision :: Dcr,Dphir,Dpsir !Difference between current and residual values
double precision :: moduleEr,moduleSigDSig
double precision :: EpsPEq,EpsPEq1,EpsPEq2 !Equivalent Plastic Deformation
double precision :: DEpsPEq !Derivative Strain in function of Equivalent Plastic Deformation
double precision :: IErateI, IErate0I, Erate_dev, Erate_dev_0, Erate_vol, Erate_vol_0, dErate_eff ! 2 Norm of the strain rate tensor, time increment
double precision :: dGdErate, dKdErate, dcdErate, dphidErate, dpsidErate, kappa
double precision, dimension(2) :: DummyArg
double precision, dimension(6) :: SigYield, SigYield2
double precision, dimension(6) :: DSigPP,DSigP1,DSigP2
double precision, dimension(6) :: DEpsPP,DEpsPP1,DEpsPP2
double precision, dimension(6) :: DEpsS,DEpsSS,DSigE
double precision, dimension(6) :: EpsP1,EpsP2
double precision, dimension(6) :: dEpsEqdEpsP,dEpsEqdEpsP1
double precision, dimension(6) :: sumSg,Er, dErate, dNratedErate, DErateS, ErateYield,DErateSS
!double precision, dimension(6) :: DIErateIDEpsRate !Derivative of the norm of the strain rate respect to the strain rate
double precision, dimension(3) :: dChisdEpsEq,dChisdEpsEq1 !Variation of softening parameters (c,phi,psi) in function of plastic strain
double precision, dimension(3) :: DChidErateef, DChidErateef1, DchiErate ! variation of state parameters with strain rate norm
double precision, dimension(3,6):: DChidErate
double precision, dimension(6,6):: DE
!double precision, dimension(6,6) :: Q ! rotational matrix to accommodate elastic and plastic strains
logical :: IsElasticUnloading = .false.
logical :: IsFailedStep = .false.
!In variables
integer, intent(in) :: IntGlo, N_S !Global ID of Gauss point or particle
double precision, intent(in) :: G_0, K_0 !Elastic Parameters
double precision, intent(in) :: Dtime !Time step
double precision, intent(in) :: cp,cr,phip,phir,psip,psir,eta !strain Softening parameter
double precision, intent(in) :: alpha_K, alpha_G, alpha_chi, alpha_beta, RefRate
double precision, intent(in), dimension(6) :: Sig0 !Initial Stress
double precision, intent(in), dimension(6) :: DEps !Incremental total strain
double precision, intent(in), dimension(6) :: Erate0 !Strain rate vector
logical, intent(in):: switch_smooth
!Inout variables
double precision, intent(inout):: D1,D2,c,phi,psi, G, K, SUM_rate !cohesion,friction angle and dilatancy angle
double precision, intent(inout), dimension(6) :: EpsP !Accumulated Plastic Strain
double precision, intent(inout), dimension(6) :: SigC !Final Stress
double precision, intent(inout), dimension(6) :: ERate !Incremental Elastic Stress
double precision, intent(inout):: Error_Euler_max, Error_Yield_last, Error_Yield_max
!Out variables
integer, intent(out) :: IPL
double precision, intent(out), dimension(6) :: DEpsP !Incremental plastic strain
logical, intent(inout):: switch_yield
integer, intent(inout):: N_i
!************************************************************************************************************
!Initialization
DEpsPEq = 0.0d0
EpsPEq = 0.0d0
SigYield = 0.0d0
DEpsP = 0.0d0
F = 0.0d0
it = 0
if (G==0.0) then
G=G_0
endif
if (K==0.0) then
K=K_0
endif
if (c==0.0) then
c=cp
endif
if (phi==0.0) then
phi=phip
endif
if (psi==0.0) then
psi=psip
endif
!Tolerances
SSTOL = 0.001d0 !Tolerance Relative Error (10-3 to 10-5)
YTOL = 0.000000001d0 !Tolerance Error on the Yield surface (10-6 to 10-9)
SPTOL = 0.0001d0 !Tolerance Softening Parameters (0.0001d0)
LTOL=0.01d0 !Tolerance for elastic unloading
MAXITER=20
ctol = abs(cp-cr)*SPTOL
phitol = abs(phip-phir)*SPTOL
psitol = abs(psip-psir)*SPTOL
DTmin = 0.000000001d0
call YieldFunctionValue(IntGlo,Sig0,c,phi,F0) ! initial yield function value
!Computes current effective strain rate tensor
call TwoNormTensor(Erate,6,IErateI)
! Compute a smoothed strain rate if switch_smooth is true
if (switch_smooth) then
N_i=N_i+1
if (N_i<N_S) then !Not enough values
SUM_rate=SUM_rate+IErateI !Accumulate the strain rate
IErateI=SUM_rate/N_i !Takes the average
else !Enough values
call TwoNormTensor(Erate0,6,IErate0I) !previous step strain rate
SUM_rate=SUM_rate*(1.0-1.0/N_S) !Approximate sum without one term
SUM_rate=SUM_rate+IErateI
IErateI=SUM_rate/N_S !Averaged strain rate
endif
call TwoNormTensor(Erate,6,IErate0I)! Norm of the uncorrected strain rate
Erate=(IErateI/IErate0I)*Erate !Corrected strain rate tensor
endif
!Computes the increment of strain rate
dErate= Erate-Erate0
!*****************************************************************************************
!Store state variables
Gu=G
Ku=K
phiu=phi
psiu=psi
cu=c
!Compute the norm of dErate
call CalculateEpsPEq(Erate,Erate_dev)
call CalculateEpsPEq(Erate0,Erate_dev_0)
!For the bulk modulus use the volumetric strain rate
Erate_vol=abs(Erate(1)+Erate(2)+Erate(3))
Erate_vol_0=abs(Erate0(1)+Erate0(2)+Erate0(3))
!call TwoNormTensor(Erate0,6,IErate0I) !previous step strain rate
dErate_eff=Erate_dev-Erate_dev_0
call check4crossing(Erate_dev_0, Erate_dev, dErate_eff, RefRate,ApplyStrainRateUpdates)
!Update G and K
if (ApplyStrainRateUpdates) then
call UpdatePardue2StrainRate(alpha_G, Erate_dev_0, Erate_dev, dErate_eff, RefRate, G_0, Gu)
end if
call CalculateEpsPEq(EpsP,EpsPEq) !Compute equivalent plastic strain
!Update state parameters
if (ApplyStrainRateUpdates) then
call UpdateSatePardue2StrainRate(EpsPEq, alpha_chi, alpha_beta, Erate_dev_0, Erate_dev,&
dErate_eff, RefRate,cp, phip, psip, cr, phir, psir,&
eta, cu,phiu, psiu)
end if
dErate_eff=Erate_vol-Erate_vol_0
call check4crossing(Erate_vol_0, Erate_vol, dErate_eff, RefRate,ApplyStrainRateUpdates)
if (ApplyStrainRateUpdates) then
call UpdatePardue2StrainRate(alpha_K, Erate_vol_0, Erate_vol, dErate_eff, RefRate, K_0, Ku)
endif
if (cu < cr.or.phiu < phir.or.psiu < psir) then !Check that the state parameters are never lower than the residual
cu = max(c,cr)
phiu = max(phi,phir)
psiu = max(psi,psir)
end if
!***********************************************************************************
! Fill elastic material matrix
D1 = Ku+(4*Gu/3)
D2 = Ku-(2*Gu/3)
DE = 0.0
DE(1:3,1:3) = D2
DE(1,1) = D1
DE(2,2) = D1
DE(3,3) = D1
DE(4,4) = Gu
DE(5,5) = Gu
DE(6,6) = Gu
!***********************************************************************************
! elastic stress increment
Call MatVec( DE, 6, DEps, 6, DSigE)
! elastic stress
Call AddVec( Sig0, DSigE, 1d0, 1d0, 6, SigC )
!***********************************************************************************
!***********************************************************************************
!Check for elasticity or plasticity loading
!Check the yield function value
call YieldFunctionValue(IntGlo,SigC,cu,phiu,F)
!If F<0 then the behavior is elastic --> Return
if (F <= YTOL) then
G=Gu
K=Ku
c=cu
phi=phiu
psi=psiu
switch_yield=.false.
IPL=0.0d0
return
end if
!If F>0, the behavior is elasto-plastic --> Continue
!***************************************************************************************
Dcr = abs(c - cr)
Dphir = abs(phi - phir)
Dpsir = abs(psi - psir)
!Check if we are in residual conditions or in softening conditions
if (Dcr <= ctol.and.Dphir <= phitol.and.Dpsir <= psitol) then
IPL = 1 !IPL=1 Residual conditions --> no changes of the strength parameters
c = cr
phi = phir
psi = psir
else
IPL = 2 !IPL=2 Softening Conditions --> changes of the strength parameters
end if
!****************************************************************************************
!Determine the proportion (alpha) of the stress increment that lies within the yield function.
SigYield = Sig0
DEpsS=DEps
ErateYield=Erate0
if (F0 < -YTOL) then !In this Time increment there is part of elastic behavior
call NewtonRaphson(SigYield, phip, phir, cp, cr, psip, psir, G_0, K_0, &
alpha_chi, alpha_G, alpha_K, alpha_beta, DEpsS, ErateYield, Erate &
, Erate_dev_0, Erate_dev, Erate_vol, Erate_vol_0,refrate, phi, c, psi, K, G, alpha, &
eta, EpsPEq, F0, MAXITER, YTOL)
else
!Determine Elastic Unloading path
call CheckElasticUnloading(IntGlo, Sig0, DSigE,c,phi,psi, IsElasticUnloading, LTOL)
if (IsElasticUnloading) then
call NewtonRaphson(SigYield, phip, phir, cp, cr, psip, psir, G_0, K_0, &
alpha_chi, alpha_G, alpha_K, alpha_beta, DEpsS, ErateYield, Erate &
, Erate_dev_0, Erate_dev, Erate_vol, Erate_vol_0, refrate, phi, c, psi, K, G, alpha, &
eta, EpsPEq, F0, MAXITER, YTOL)
else
alpha = 0.0d0 !In this time increment all is plastic behavior
end if
end if
!****************************************************************************************
!Determine the elastic portion of the stress increment
!DSigE = alpha * DSigE !Correct Incremental Elastic Stress
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Determine the plastic portion of the stress increment.
!The method used is the MODIFIED EULER INTEGRATION SCHEME with error control
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dErate=Erate-ErateYield
!Sub-stepping parameters
T = 0.0d0
DT = 1.0d0
m = 0 !Counter
!Start the yielding
Do while (T <= 1.0d0)
Rn = 100
!DrealTime=DT*Dtime ! Real time increment of sub-step
call CalculateEpsPEq(EpsP,EpsPEq) !Determine Equivalent plastic Strain (EpsPEq)
!****************************************************************************************
!Calculate first and second estimate of stress increment and compute relative error
!1)Calculation of the portion of the plastic strain increment (DEpsPP)
DEpsSS = DT * DEpsS !Portion of the plastic strain increment
DErateSS= DT * dErate
!Calculate the new increments of strain rate
call CalculateEpsPEq(ErateYield, Erate_dev_0)
Erate_vol_0=abs(ErateYield(1)+ErateYield(2)+ErateYield(3))
!update the strain rate
ErateYield=ErateYield+DErateSS
call CalculateEpsPEq(ErateYield, Erate_dev)
Erate_vol=abs(ErateYield(1)+ErateYield(2)+ErateYield(3))
dErate_eff=Erate_dev-Erate_dev_0
call check4crossing(Erate_dev_0, Erate_dev, dErate_eff, RefRate,ApplyStrainRateUpdates)
!Calculate a first estimate of the associated stress
!hardening/softening parameter changes
call PartialChiToPartialEpsilonEq(Erate_dev_0,alpha_chi, alpha_beta,eta,cp,cr,phip,phir,psip,psir, &
EpsPEq,dChisdEpsEq,RefRate)
call PartialEpsEqToPartialEpsP(EpsP,EpsPEq,dEpsEqdEpsP)
!Actual increment of stress calculation
call DetermineDSigAndDEpsP(IntGlo,G, K,G_0, K_0, D1, D2, alpha_k, alpha_G, alpha_chi, alpha_beta,&
Erate_vol_0,Erate_vol, Erate_dev_0, Erate_dev, RefRate, c,phi,psi, cp, phip, psip, EpsPEq, eta,SigYield,&
dEpsEqdEpsP,dChisdEpsEq, dErate_eff, DEpsSS,DSigP1,DEpsPP1, DchiErate)!First increment of stress
!First plastic strain
EpsP1 = EpsP + DEpsPP1
call CalculateEpsPEq(EpsP1,EpsPEq1) !Determine Equivalent plastic Strain (EpsPEq)
call PartialChiToPartialEpsilonEq(Erate_dev,alpha_chi, alpha_beta,eta,cp,cr,phip,phir,psip,psir, &
EpsPEq1,dChisdEpsEq1, RefRate)
call PartialEpsEqToPartialEpsP(EpsP1,EpsPEq1,dEpsEqdEpsP1)
!Store initial state variables values
c1 = c
phi1 = phi
psi1 = psi
!Update the state variables according to the strain rate
if (ApplyStrainRateUpdates) then
c1=c1+DchiErate(1)
phi1=phi1+DchiErate(2)
psi1=psi1+DchiErate(3)
end if
call UpdateStateParameters(dChisdEpsEq1,dEpsEqdEpsP1,DEpsPP1,c1,phi1,psi1)
!*****************************************************************************************
!2)Calculate a second estimate of the associated stress
!hardening/softening parameter changes
SigYield2 = SigYield + DSigP1
call DetermineDSigAndDEpsP(IntGlo,G, K,G_0, K_0, D1, D2, alpha_k, alpha_G, alpha_chi, alpha_beta,&
Erate_vol_0,Erate_vol,Erate_dev_0, Erate_dev, RefRate, c1,phi1,psi1,cp, phip, psip, EpsPEq1, eta,SigYield2,&
dEpsEqdEpsP1,dChisdEpsEq1, dErate_eff, DEpsSS,DSigP2,DEpsPP2, DchiErate)!First increment of stress
!Second plastic strain
EpsP2 = EpsP + DEpsPP2
call CalculateEpsPEq(EpsP2,EpsPEq2) !Determine Equivalent plastic Strain (EpsPEq)
call PartialChiToPartialEpsilonEq(Erate_dev,alpha_chi, alpha_beta,eta,cp,cr,phip,phir,psip,psir, &
EpsPEq2,dChisdEpsEq1, RefRate)
call PartialEpsEqToPartialEpsP(EpsP2,EpsPEq2,dEpsEqdEpsP1)
!store initial state variables for second estimation
c2 = c
phi2 = phi
psi2 = psi
if (ApplyStrainRateUpdates) then
c2=c2+DchiErate(1)
phi2=phi2+DchiErate(2)
psi2=psi2+DchiErate(3)
end if
call UpdateStateParameters(dChisdEpsEq1,dEpsEqdEpsP1,DEpsPP2,c2,phi2,psi2)
!****************************************************************************************
!Increment os stress is computed as the averages of the two calculations
DSigPP = 0.5d0 * (DSigP1 + DSigP2)
!*****************************************************************************************
!Calculation of the relative error
Er = 0.5d0 * (DSigP1 - DSigP2)
moduleEr = sqrt(Er(1)*Er(1)+Er(2)*Er(2)+ Er(3)*Er(3)+ Er(4)*Er(4)+Er(5)*Er(5)+Er(6)*Er(6))
sumSg = SigYield + DSigPP
moduleSigDSig = sqrt(sumSg(1)*sumSg(1) + sumSg(2)*sumSg(2) + sumSg(3)*sumSg(3)+ &
sumSg(4)*sumSg(4) + sumSg(5)*sumSg(5) + sumSg(6)*sumSg(6))
!Check the relative error (Rn) of the new stresses, with the defined tolerance (SSTOL)
Rn = (moduleEr/moduleSigDSig)
!store error
if (Rn>Error_Euler_max) Error_Euler_max=Rn
!************************************************************************************
!4)If Rn>SSTOL, the loop is not finished and the sub-step is recalculated smaller
if (Rn > SSTOL.and.m < 50) then
beta = max (0.9d0*(sqrt(SSTOL/Rn)), 0.1d0)
DT = max (DT*beta, DTmin)
m = m + 1 !Update counter
IsFailedStep= .true.
ErateYield=ErateYield-DErateSS
else
!Update the accumulated stresses, plastic strain and softening parameters
SigYield = SigYield + DSigPP
DEpsPP = 0.5d0 * (DEpsPP1 + DEpsPP2)
DEpsP = DEpsP + DEpsPP
EpsP = EpsP + DEpsPP
call CalculateEpsPEq(EpsP,EpsPEq) !Determine Equivalent plastic Strain (EpsPEq)
call PartialEpsEqToPartialEpsP(EpsP,EpsPEq,dEpsEqdEpsP)
if (ApplyStrainRateUpdates) then
call UpdatePardue2StrainRate(alpha_G, Erate_dev_0, Erate_dev, dErate_eff, RefRate, G_0, G)
call UpdatePardue2StrainRate(alpha_K,Erate_vol_0, Erate_vol, dErate_eff, RefRate, K_0, K)
c=c+DchiErate(1)
phi=phi+DchiErate(2)
psi=psi+DchiErate(3)
end if
call UpdateStateParameters(dChisdEpsEq,dEpsEqdEpsP,DEpsPP,c,phi,psi)
D1 = K+(4*G/3)
D2 = K-(2*G/3)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!! END OF STEP CORRECTION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!**********************************************************************************
!Check if we are on/in the yield surface, otherwise we are still outside (F>0)
!and a correction is needed.
call YieldFunctionValue(IntGlo,SigYield,c,phi,F)
!store yield function errors
Error_Yield_last=F
if (abs(F)>abs(Error_Yield_max)) Error_Yield_max=F
n=0 !Counter
do while (abs(F) > YTOL.and.n < 10) !The correction is needed
n = n + 1
call CalculateEpsPEq(EpsP,EpsPEq)!Determine Equivalent plastic Strain (EpsPEq)
call PartialChiToPartialEpsilonEq(Erate_dev,alpha_chi, alpha_beta,eta,cp,cr,phip,phir,psip,psir,&
EpsPEq,dChisdEpsEq,RefRate)
call PartialEpsEqToPartialEpsP(EpsP,EpsPEq,dEpsEqdEpsP)
call EndOfStepCorrection(IntGlo,D1,D2,G,IPL,F,SigYield,dChisdEpsEq,dEpsEqdEpsP,EpsP, &
c,phi,psi, DEpsPP,DEpsSS)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*********Psudo time updating*********************************************************
beta = min (0.9d0*(sqrt(SSTOL/Rn)), 1.1d0)
if (IsFailedStep) then !previous step failed
IsFailedStep = .false.
beta = min (beta, 1.0d0)
end if
!The substep is updated
T1 = T + DT
DT = beta * DT
DT = max (DT, DTmin)
DT = min (DT, 1.0d0-T1)
!If T1>1 the calculation is finished
If (T1 >= 1d0) then
SigC = SigYield !Determine Final stresses
switch_yield=.true. ! Because it is yielding, no check on state variables is needed
return
end if
T = T1
!**********************************************************************************
end if
end do !If T=1 the loop is finished
SigC = SigYield
switch_yield=.true. !Because it is yielding, no check on state variables is needed
!end do
end subroutine SSMC_Strain_Rate
Subroutine CheckElasticUnloading(IntGlo, Sig, DSig,c,phi,psi, IsElasticUnloading, LTOL)
!**********************************************************************
!
! It returs true if the stress path correspond to elastic unloading
!
!**********************************************************************
implicit none
!Local Variables
double precision :: p,J,Lode,S3TA, NormDSig, NormDFDSig, DSigdpDFDSig, CTeta
double precision, dimension(6) :: DFDsig, dPdSig
!In Variables
double precision, intent(in), dimension(6) :: Sig, DSig
double precision, intent(in) :: c, phi, psi, LTOL
integer, intent(in) :: IntGlo !Global ID of Gauss point or particle
!Out Variables
Logical, intent(out) :: IsElasticUnloading
!Calculation of the invariants (p',J,Lode)
call CalculateInvariants(IntGlo,Sig,p,J,Lode,S3TA)
!Calulate yield function derivative with respect to stress
call PartialForGtoPartialSigma(Sig,p,J,Lode,S3TA,c,phi,psi,DFDSig,dPdSig)
!calculate the 2-norms of the tensors
call TwoNormTensor(Dsig, 6, NormDSig)
call TwoNormTensor(DFDSig, 6, NormDFDSig)
!Calculate Dsig:DFDSig
call MatAdpMatB(DSig, DFDSig, 6, DSigdpDFDSig)
!Cos of angle between stress path and normal vector to yield surface
CTeta=DSigdpDFDSig/(NormDSig*NormDFDSig)
If (CTeta < -LTOL) then !Elastic Unloading
IsElasticUnloading=.true.
else
IsElasticUnloading=.false.
end if
end subroutine CheckElasticUnloading
Subroutine NewtonRaphson(Sig0, phip, phir, cp, cr, psip, psir, G_0, K_0, &
alpha_chi, alpha_G, alpha_K, alpha_psi, dEps, Erate_0, Erate &
, IErate_0I, IErateI, Erate_vol, Erate_vol_0, refrate, phi, c, psi, K, G, alpha, &
eta, Epspq, F0, MAXITER, FTOL)
!_____________________________________________________________________________________
! This subroutine determine the elastic proportion with consistent state and elastic
! parameters
!_____________________________________________________________________________________
implicit none
!Input variables
double precision, dimension(6), intent(in):: Erate
double precision, intent(in):: phip, phir, cp, cr, psip, psir, G_0, K_0, &
alpha_chi, alpha_G, alpha_K, alpha_psi, refrate, &
IErateI, IErate_0I,Erate_vol, Erate_vol_0, FTOL, Epspq, eta
integer, intent(in):: MAXITER
!output variables
double precision, dimension(6), intent(inout):: Sig0, dEps, Erate_0
double precision, intent(inout):: phi, c, psi, K, G, alpha, F0
!local variables
double precision:: FT, dG, dK, dc, dphi, dpsi, K_u, G_u, phi_u, c_u, psi_u, &
derate_alpha(6), dEps_alpha(6), Erate_alpha(6), IErateI_alpha, dErate_eff, &
dSig_vel(6), DE(6,6), D1, D2, SIg_alpha(6), p, J, Lode, S3TA, DFDSig(6), &
dPPdSig(6), dFdChis(2), F_prime, aux(6,6), dDdK(6,6), dDdG(6,6), dSigdAlpha(6), &
Erate_vol_alpha
integer:: n, I
logical:: ApplyStrainRateUpdates
FT=1000
alpha=1.0d0
n=0
!Compute change due to strain rate
dErate_eff=IErateI-IErate_0I
call check4crossing(IErate_0I, IErateI, dErate_eff, RefRate,ApplyStrainRateUpdates)
if (ApplyStrainRateUpdates) then !Check if the strain rates are enough to cause updating
dG=G_0*alpha_G*log10(IErateI/IErate_0I)
dphi=phip*alpha_chi*log10(IErateI/IErate_0I)*exp(-eta*Epspq)
dc=cp*alpha_chi*log10(IErateI/IErate_0I)*exp(-eta*Epspq)
dpsi=psip*alpha_psi*log10(IErateI/IErate_0I)*exp(-eta*Epspq)
else
dG=0.0d0
dphi=0.0d0
dc=0.0d0
dpsi=0.0d0
endif
dErate_eff=Erate_vol-Erate_vol_0
call check4crossing(Erate_vol_0, Erate_vol, dErate_eff, RefRate,ApplyStrainRateUpdates)
if (ApplyStrainRateUpdates) then
dK=K_0*alpha_K*log10(Erate_vol/Erate_vol_0)
else
dK=0.0d0
endif
do while ((abs(FT)>=FTOL).and.(n<=MAXITER))
!____ store variables_______________________________________________________________________
K_u=K
G_u=G
phi_u=phi
c_u=c
psi_u=psi
F0=FT
n=n+1
!___Compute trial strain and strain rate____________________________________________________
derate_alpha=alpha*(Erate-Erate_0)
dEps_alpha=alpha*dEps
Erate_alpha=Erate_0+derate_alpha
call CalculateEpsPEq(Erate_alpha,IErateI_alpha)
dErate_eff=IErateI_alpha-IErate_0I
!__ Update state and elastic parameters_____________________________________________________
call check4crossing(IErate_0I, IErateI_alpha, dErate_eff, RefRate,ApplyStrainRateUpdates)
if (ApplyStrainRateUpdates) then !Update state parameters due to strain rate
call UpdatePardue2StrainRate(alpha_G, IErate_0I, IErateI_alpha, dErate_eff, RefRate, G_0, G_u)
call UpdateSatePardue2StrainRate(Epspq, alpha_chi, alpha_psi, IErate_0I, IErateI_alpha, dErate_eff, RefRate,&
cp, phip, psip,cr, phir, psir, eta, c_u,phi_u, psi_u)
endif
!compute volumetric strain rates
Erate_vol_alpha=abs(Erate_alpha(1)+Erate_alpha(2)+Erate_alpha(3))
dErate_eff=Erate_vol_alpha-Erate_vol_0
call check4crossing(Erate_vol_0, Erate_vol_alpha, dErate_eff, RefRate,ApplyStrainRateUpdates)
if (ApplyStrainRateUpdates) then
call UpdatePardue2StrainRate(alpha_K,Erate_vol_0, Erate_vol_alpha, dErate_eff, RefRate, K_0, K_u)
endif
!__ Update trial stress________________________________________________________________________
! Fill elastic material matrix
D1 = K_u+(4*G_u/3)
D2 = K_u-(2*G_u/3)
DE = 0.0
DE(1:3,1:3) = D2
DE(1,1) = D1
DE(2,2) = D1
DE(3,3) = D1
DE(4,4) = G_u
DE(5,5) = G_u
DE(6,6) = G_u
! elastic stress increment
Call MatVec( DE, 6, dEps_alpha, 6, dSig_vel)
! elastic stress
Call AddVec( Sig0, dSig_vel, 1d0, 1d0, 6, Sig_alpha )
!__ Evaluate Yield function and derivatives____________________________________________________
call YieldFunctionValue(1,Sig_alpha,c_u,phi_u,FT)
call CalculateInvariants(1,Sig_alpha,p,J,Lode,S3TA)
call PartialForGtoPartialSigma(Sig_alpha,p,J,Lode,S3TA,c_u,phi_u,psi_u,DFDSig,dPPdSig)
call PartialFtoPartialChi(p,J,Lode,S3TA,c_u,phi_u,dFdChis)
!__ Update alpha________________________________________________________________________________
F_prime=dFdChis(1)*dc+dFdChis(2)*dphi
dDdK=0.0
dDdK(1:3,1:3) = 1.0
dDdG=0.0
dDdG(1:3,1:3)=-2./3.
dDdG(1,1)=4./3.
dDdG(2,2)=4./3.
dDdG(3,3)=4./3.
dDdG(4,4)=1.0
dDdG(5,5)=1.0
dDdG(5,5)=1.0
aux=alpha* (dK*dDdK+dG*dDdG)
aux=DE+aux
call MatVec(aux, 6, dEps, 6, dSigdAlpha)
do I=1,6 !dot product
F_prime=F_prime+dFdSig(I)*dSigdAlpha(I)
enddo
alpha=alpha-FT/F_prime
enddo
K=K_u
G=G_u
phi=phi_u
c=c_u
psi=psi_u
Sig0=Sig_alpha
dEps=(1-alpha)*dEps
Erate_0=Erate_0+(1-alpha)*(Erate-Erate_0)
end subroutine NewtonRaphson
Subroutine YieldFunctionValue(IntGlo,Sig,c,phi,F)
!**********************************************************************
!
! In this subroutine the yield function evaluated is a smooth hyperbolic approximation to the
! Mohr-Coulomb yield criterion (Abbo and Sloan, 1995).
!
! The edges of the hexagonal pyramid and the tip have been smoothed.
! There are two parameters aSmooth (smoothes the tip) and ATTRAN(smoothes the edges)
! In this case aSmooth=0.0005*c*cot(phi) and LodeT=25�.
! If aSmooth=0 and LodeT=30� the original Mohr-Coulomb is obtained.
!
!**********************************************************************
implicit none
!Local variables
double precision :: p,J,Lode,S3TA !Invariants
double precision :: COH, SPHI, CPHI, COTPHI, STA, CTA, K, aSmooth, ASPHI2, SGN, A, B
double precision, parameter :: C00001 = 1.0d0 !Parameters
double precision, parameter :: C00003 = 3.0d0
double precision, parameter :: C00P50 = 0.0005d0
double precision, parameter :: C00000 = 0.0d0
double precision, parameter :: C00IR3 = 0.577350269189626d0
double precision, parameter :: C000P1 = 0.00000000001D0
!Constants for rounded K function (for LodeT=25)
!double precision, parameter :: A1 = 1.432052062044227d0
!double precision, parameter :: A2 = 0.406941858374615d0
!double precision, parameter :: B1 = 0.544290524902313d0
!double precision, parameter :: B2 = 0.673903324498392d0
!double precision, parameter :: ATTRAN = 0.436332312998582d0 !Smoothing parameter: LodeT in radians
!Constants for rounded K function (for LodeT=29.5)
double precision, parameter :: A1 = 7.138654723242414d0
double precision, parameter :: A2 = 6.112267270920612d0
double precision, parameter :: B1 = 6.270447753139589d0
double precision, parameter :: B2 = 6.398760841429403d0
double precision, parameter :: ATTRAN = 0.514872129338327d0 !Smoothing parameter: LodeT in radians
!Constants for rounded K function (for LodeT=30)
!double precision, parameter :: A1 = -138300705.446275
!double precision, parameter :: A2 = -138300706.472675
!double precision, parameter :: B1 = -138300706.3123
!double precision, parameter :: B2 = 0.192450089729875
!double precision, parameter :: ATTRAN = 0.523598776 !Smoothing parameter: LodeT in radians
!In variables
double precision, intent(in), dimension(6) :: Sig
double precision, intent(in) :: c,phi
integer, intent(in) :: IntGlo !Global ID of Gauss point or particle
!Out variables
double precision, intent(out) :: F
F = C00000
!Calculation of the invariants (p',J,Lode)
call CalculateInvariants(IntGlo,Sig,p,J,Lode,S3TA)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Evaluation of the yield function with Smoothing!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Material parameters
COH = c !Cohesion
SPHI = sin(phi)
CPHI = cos(phi)
COTPHI = CPHI/SPHI
aSmooth = C00P50*COH*COTPHI !Smoothing parameter
ASPHI2 = aSmooth*aSmooth*SPHI*SPHI
if (abs(phi) == C00000) then
ASPHI2 = C00P50*C00P50*COH*COH*CPHI*CPHI
end if
!Calculate K function
if (abs(Lode) < ATTRAN) then
STA = sin(Lode)
CTA = cos(Lode)
K = CTA - STA*SPHI*C00IR3
else
SGN = SIGN(C00001,Lode)
A = A1 + A2*SGN*SPHI
B = B1*SGN + B2*SPHI
K = A - B*S3TA
end if
!Calculate value of Hyperbolic Yield function
F = p*SPHI + sqrt(J*J*K*K+ASPHI2) - COH*CPHI
end subroutine YieldFunctionValue
Subroutine PartialForGtoPartialSigma(Sig,p,J,Lode,S3TA,c,phi,psi,DFDSig,dPPdSig)
!**********************************************************************
!
! Calculation of the derivatives of the yield function (F) and the plastic potencial punction (P).
! Based on Abbo & Sloan (1995)
!
!**********************************************************************
implicit none
!Local variables
integer :: i
double precision :: COH, SPHI, CPHI, TPHI, COTPHI, STA, CTA, A, B, &
D, aSmooth, ASPHI2, SGN, T3TA, C3TA, J2, psi2
double precision :: K, dKdLode
double precision :: SPSI, CPSI, TPSI, COTPSI, ASPSI2
double precision :: i1, i2, Sx, Sy, Sz
double precision :: DFDp,DFDJ,DFDLode !Derivatives F respect Invariants
double precision :: DPDp,DPDJ,DPDLode !Derivatives P respect Invariants
double precision :: C1, C2, C3
double precision, dimension(6):: DpDSig,DJDSig,DJ3DSig !Derivatives Invariants
double precision, parameter :: C00001 = 1.0D0 !Parameters
double precision, parameter :: C000P5 = 0.5D0
double precision, parameter :: C00P50 = 0.0005D0
double precision, parameter :: C00000 = 0.0D0
double precision, parameter :: C00003 = 3.0D0
double precision, parameter :: C00004 = 4.0D0
double precision, parameter :: C00002 = 2.0D0
double precision, parameter :: CP3333 = 0.333333333333333D0
double precision, parameter :: C00IR3 = 0.577350269189626D0
double precision, parameter :: C0R3I2 = 0.866025403784439D0
double precision, parameter :: C000P1 = 0.000000000000001D0
double precision, parameter :: J0 = 0.001D0
!Constants for rounded K function (for LodeT=25)
!double precision, parameter :: A1 = 1.432052062044227d0
!double precision, parameter :: A2 = 0.406941858374615d0
!double precision, parameter :: B1 = 0.544290524902313d0
!double precision, parameter :: B2 = 0.673903324498392d0
!double precision, parameter :: ATTRAN = 0.436332312998582d0 !Smoothing parameter: LodeT in radians
!Constants for rounded K function (for LodeT=29.5)
double precision, parameter :: A1 = 7.138654723242414d0
double precision, parameter :: A2 = 6.112267270920612d0
double precision, parameter :: B1 = 6.270447753139589d0
double precision, parameter :: B2 = 6.398760841429403d0
double precision, parameter :: ATTRAN = 0.514872129338327d0 !Smoothing parameter: LodeT in radians
!Constants for rounded K function (for LodeT=30)
!double precision, parameter :: A1 = -138300705.446275
!double precision, parameter :: A2 = -138300706.472675
!double precision, parameter :: B1 = -138300706.3123
!double precision, parameter :: B2 = 0.192450089729875
!double precision, parameter :: ATTRAN = 0.523598776 !Smoothing parameter: LodeT in radians
!In variables
double precision, intent(in) :: c,phi,psi !Soft Parameters
double precision, intent(in), dimension(6) :: Sig
!Out variables
double precision, intent(out), dimension(6) :: DFDSig, dPPdSig !Derivatives respect Sigma
!Inout variables
double precision, intent(inout) :: p,J,Lode,S3TA !Invariants
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!