forked from happyhappysyr/Phase_field_ductile
-
Notifications
You must be signed in to change notification settings - Fork 0
/
2D_Powerlaw.for
314 lines (313 loc) · 10.1 KB
/
2D_Powerlaw.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
c Author: Mojtaba Abdolkhani
c Email: mojtababdolkhani@gmail.com
SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL,
1 DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP,
2 PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS,
3 COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER,
4 KSPT, KSTEP, KINC)
c
INCLUDE 'ABA_PARAM.INC'
c
CHARACTER*20 CMNAME
c
DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS, NTENS),
1 STRAN(NTENS), DSTRAN(NTENS),PROPS(NPROPS), DROT(3, 3),TIME(2),
2 EPLAS(NTENS), FLOW(NTENS)
c
PARAMETER(zero=0.D0, one=1.D0, two=2.D0, three=3.D0)
c
parameter (six = 6.d0, newton = 100, dsqrt23 = dsqrt(two/three),
1 toler = 1.0d-8, zeta = 1.d0)
c
DIMENSION one_hydro(ntens,ntens), dstran_trace(ntens),
1 one_dev(ntens,ntens), STRESS_MINUS(NTENS),
2 DDSDDE_MINUS(NTENS,NTENS)
c
REAL*8 EMOD,ENU,EBULK3,EBULK,EG,EG2,EG3,ELAM
C ======================================================================
c
c The energy decomposition is triggered via material PROPS(10) in input.
c Two options are available: ELASTOPLASTIC_ISO(41) or ELASTOPLASTIC_SD(40)
c
c ELASTOPLASTIC_ISO - no deformation energy decomposition
c ELASTOPLASTIC_SD - spherical-deviatoric energy decomposition
C ======================================================================
phi=temp+dtemp
c
c
C ======================================================================
!Material properties from input file
EMOD=PROPS(1)
ENU=PROPS(2)
xl=props(3) ! Phase field length scale
Gc=props(4) ! Toughness
type2 = PROPS(10) ! SD(40) OR ISO(41)
xK = 0.0001d0
c
!Degradation function
degfnc = zero
degfnc = (one - xK)*(one - phi)**2 + xK
dg=-2.d0*(1.d0-phi)
ddg=2.d0
c
!Computing lame constants
EBULK3=EMOD/(one-two*ENU)
EBULK = EBULK3/three
EG2=EMOD/(one+ENU)
EG=EG2/two
EG3=three*EG
ELAM=(EBULK3-EG2)/three
c
C ======================================================================
!Building material matrix
DO 20 K1=1,NTENS
DO 11 K2=1,NTENS
DDSDDE(K2,K1)=0.0
11 CONTINUE
20 CONTINUE
C
DO 40 K1=1,NDI
DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM
30 CONTINUE
DDSDDE(K1,K1)=EG2+ELAM
40 CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50 CONTINUE
C ======================================================================
!RECOVER PLASTIC STRAINS AND ROTATE FORWARD
!RECOVER EQUIVALENT PLASTIC STRAIN
c
eqplas = statev(1)
call rotsig(statev(1+1:ntens), drot, eplas, 2, ndi, nshr)
C ======================================================================
c
STRAN = STRAN + DSTRAN
STRESS = zero
c
!Stress predictor
DO K1=1, NTENS
DO K2=1, NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2, K1)*(STRAN(K1)-EPLAS(K1))
END DO
END DO
c
if (type2.eq.40) then
! DECOMPOSITION - VOLUMETRIC DEVIATORIC
ONE_hydro = zero
A_trace = zero
dstran_trace = zero
c
do i = 1, ndi
do j = 1, ndi
ONE_hydro(i,j) = one/three
end do
end do
c
ONE_Dev = zero
do i = 1, ntens
do j = 1, ntens
if (i .eq. j) ONE_dev (i,j) = ONE
end do
end do
c
ONE_dev = ONE_dev - ONE_hydro
dstran_trace = matmul(ONE_hydro, STRAN)
c
if (dstran_trace(1) .le. zero) then
A_trace = one
c
DDSDDE_MINUS = A_trace*matmul(ONE_hydro, DDSDDE)
DDSDDE = DDSDDE - DDSDDE_MINUS
c
STRESS_MINUS = A_trace*matmul(ONE_hydro, STRESS)
STRESS = STRESS - STRESS_MINUS
else
STRESS_MINUS = zero
DDSDDE_MINUS = zero
end if
c
else if (type2.eq.41) then
STRESS_MINUS = zero
DDSDDE_MINUS = zero
A_trace = zero
else
write(6,*)""
write(6,*)"WARNING: CHECK MATERIAL NAME!"
end if
C ======================================================================
c
!CALCULATE EQUIVALENT VON MISES STRESS
c
smises = (stress(1) - stress(2))**2 + (stress(2) - stress(3))**2
1 + (stress(3) - stress(1))**2
do k1 = ndi + 1, ntens
smises = smises + six*stress(k1)**2
enddo
c
smises = dsqrt(smises/two)*dsqrt23
C ======================================================================
!GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
c
SYIEL0 = PROPS(5)
EPS0 = PROPS(6)
xn = PROPS(7)
SYIELD_max = PROPS(8)
jeltype1 = PROPS(9)
c
SYIELD = SYIELD_max + SYIEL0*(EPS0 + eqplas)**xn
C ======================================================================
c
!DETERMINE ACTIVE YIELDING
c
if (smises .gt. (one + toler)*syield*dsqrt23) then
c
!ACTIVE YIELDING
!SEPARATE HYDROSTATIC STRESS FROM DEVIATORIC
!CALCULATE FLOW DIRECTION
c
shydro = (stress(1) + stress(2) + stress(3))/three
onesy = one/smises
do k1 = 1, ndi
flow(k1) = onesy*(stress(k1) - shydro)
enddo
c
do k1 = 1 + ndi, ntens
flow(k1) = stress(k1)*onesy
enddo
C ======================================================================
c
!SOLVE FOR EQUIVALENT VON MISES STRESS AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
c
syiel_0 = syield
dgama = zero
hard = zero
c
do 130 knewton = 1, newton
c
rhs = smises - two*EG*dgama - dsqrt23*syield
dgama = dgama + rhs/(two*EG + (two/three)*hard)
c
SYIELD = SYIELD_max+SYIEL0*(EPS0+eqplas+dsqrt23*dgama)**xn
HARD = xn*SYIEL0*(EPS0 + eqplas+dsqrt23*dgama)**(xn-one)
c
if (abs(rhs) .lt. toler*syiel0*dsqrt23) go to 10
c
130 CONTINUE
WRITE(6,2) NEWTON
2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1 'CONVERGE AFTER ',I3,' ITERATIONS')
c
PNEWDT = 0.65d0
c
10 continue
C ======================================================================
c
!UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND EQUIVALENT PLASTIC STRAIN
c
do k1 = 1, ndi
stress(K1) = stress(k1) - 2.d0*EG*dgama*FLOW(K1)
EPLAS(K1)=EPLAS(K1)+FLOW(K1)*dgama
enddo
c
do k1 = ndi + 1, ntens
stress(K1) = stress(k1) - 2.d0*EG*dgama*FLOW(K1)
EPLAS(K1)=EPLAS(K1)+FLOW(K1)*dgama*TWO
enddo
c
eqplas = eqplas + dsqrt23*dgama
c
!CALCULATE PLASTIC DISSIPATION
spd = spd + dsqrt23*dgama*(syiel_0 + syield)/two
C ======================================================================
c
!FORMULATE THE JACOBIAN (MATERIAL TANGENT)
c
DDSDDE = zero
c
theta = one - 2.d0*EG*dgama/(smises)
theta_bar = (one/(one + hard/EG3)) - (one-theta)
c
do k1 = 1, ndi
do k2 = 1, ndi
DDSDDE(k1,k2) = EBULK - 2.d0*EG*theta/three -
1 2.d0*EG*theta_bar*flow(k1)*flow(k2)
end do
DDSDDE(k1,k1) = DDSDDE(k1,k1) + 2.d0*EG*theta
end do
c
do k1 = ndi + 1, ntens
DDSDDE(k1,k1) = EG*theta
1 - EG*theta_bar*flow(k1)*flow(k1)
end do
end if
C ======================================================================
DDSDDE = degfnc*DDSDDE
STRESS = degfnc*STRESS
c
DDSDDE = DDSDDE + DDSDDE_MINUS
STRESS = STRESS + STRESS_MINUS
C ======================================================================
c
!STORE STRAINS IN STATE VARIABLE ARRAY
c
do k1 = 1, ntens
statev(1 + k1) = eplas(k1)
end do
statev(1) = eqplas
c
SSE = zero !elastic deformation energy array
c
SSE = 0.5d0*ELAM*(one-A_trace)*((stran(1)-eplas(1)) +
1 (stran(2)-eplas(2)) + (stran(3)-eplas(3)))**2 +
2 EG*((stran(1)-eplas(1))**2 + (stran(2)-eplas(2))**2 +
3 (stran(3)-eplas(3))**2) +
4 0.5d0*EG*((stran(4)-eplas(4))**2)
c
c ************************************************************************************************************
c
!SSE - Elastic energy
psiE_old = statev(2+ntens)
psiE_new = SSE
c
!SPD - Plastic energy
psiE_SPD = SPD
c
psiE_SSE = zero
psiE = zero
c
c ! AT-2 or TH model
c if (jeltype1.eq.30)then !AT-2 model
c psiE_crit = zero
c elseif (jeltype1.eq.31)then !TH model
psiE_crit = 0.2652d0*(Gc/xL)
c else
c write(6,*)"***ERROR - check element types"
c end if
c
! Elastic energy history field
if (psiE_new.gt.psiE_old) then
psiE_SSE = psiE_new
else
psiE_SSE = psiE_old
end if
c
if ((psiE_SSE + psiE_SPD).gt.psiE_crit) then
psiE = psiE_SSE + psiE_SPD - psiE_crit
psiE = zeta*psiE
else
psiE = zero
end if
c
statev(2+ntens) = psiE_SSE
c ***************************************************************************************
w=phi**2
dw=2.d0*phi
ddw=2.d0
cw=0.5d0
rpl=((psiE/psiE_crit)*((1.d0-phi)/(xl**2)))-(phi/(xl**2))
drpldt=-((1.d0+(psiE/psiE_crit))/(xl**2))
c
RETURN
END