44
45
46
49
50
51
52#include "implicit_f.inc"
53
54
55
56#include "tabsiz_c.inc"
57
58
59
60 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,NFUNC,NUMTABL,
61 INTEGER, INTENT(IN) :: NVARTMP
62 INTEGER, DIMENSION(NFUNC), INTENT(IN) :: IFUNC
63 my_real,
INTENT(IN) :: time,timestep
64 my_real,
DIMENSION(NUPARAM),
INTENT(IN) :: uparam
65 my_real,
DIMENSION(NEL),
INTENT(IN) :: gs,rho,dplanl,
66 . depsxx,depsyy,depsxy,depsyz,depszx,
67 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
68 TYPE(TABLE_4D_), DIMENSION(NUMTABL) ,INTENT(IN) :: TABLE
69 my_real,
DIMENSION(NEL),
INTENT(IN) :: loff
70
71
72
73 my_real,
DIMENSION(NEL),
INTENT(OUT) ::
74 . signxx,signyy,signxy,signyz,signzx,dezz,dpla
75
76
77
78 INTEGER, DIMENSION(NEL,NVARTMP), INTENT(INOUT) :: VARTMP
79 my_real,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) :: uvar
80 my_real,
DIMENSION(NEL),
INTENT(INOUT) :: off,yld,pla,epsd,soundsp
81
82
83
84 INTEGER :: (SNPC)
86
87
88
89 INTEGER I,II,ITER,NITER,ICONV,,NINDX
90 INTEGER ,PARAMETER :: FUNC_TRAC = 1
91 INTEGER ,PARAMETER :: FUNC_COMP = 2
92 INTEGER ,PARAMETER :: FUNC_SHEAR = 3
93 INTEGER, DIMENSION(NEL) :: INDX,IAD,ILEN
94 my_real :: lam,dlam,epsdot,da0,da1,da2,
95 . ca,cb,aa,bb,cc,a1s,a1c,a1t,a2s,a2c,a2t,e,nu,nu1
96 . yy,dx2,dfdsigdlam,yld_norm,epd_min,epd_max,
97 . normxx,normyy,normxy,normzz,normyz,normzx,alfa,alfi,dtinv,
98 . normgxx,normgyy,normgxy,normgzz,sig_dfdsig,
99 . epdt_min,epdt_max,epdc_min,epdc_max,epds_min,epds_max
100 my_real,
DIMENSION(NEL) :: ff,p,svm,svm2,ylds,sxx,syy,sxy,szz,
101 . dpxx,dpyy,dpxy,dpzz,sigt,sigc,sigs,dft,dfc,dfs,a11_2d,a12_2d,g,
102 . a0,a1,a2,nup,epspt,epspc,epsps,epdt,epdc,epds,dydx
103
104 my_real,
DIMENSION(NEL,2) :: xvec
105 my_real,
DIMENSION(NUMTABL) :: tfac
106 my_real,
DIMENSION(NFUNC) :: yfac
107 LOGICAL :: CONV(NEL)
108 my_real,
PARAMETER :: sfac = 1.05d0
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132 e = uparam(1)
133 nu = uparam(5)
134
135 nupc = uparam(9)
136
137 iconv = nint(uparam(15))
138 icas = uparam(17)
139 xfac = uparam(18)
140 alfa =
min(one, uparam(16)*timestep)
141 epdt_min = uparam(19)
142 epdt_max = uparam(20)
143 epdc_min = uparam(21)
144 epdc_max = uparam(22)
145 epds_min = uparam(23)
146 epds_max = uparam(24)
147 tfac(1) = uparam(25)
148 tfac(2) = uparam(26)
149 tfac(3) = uparam(27)
150 yfac(1) = uparam(28)
151 yfac(2) = uparam(29)
152 a11_2d(1:nel) = uparam(2)
153 a12_2d(1:nel) = uparam(3)
154 g(1:nel) = uparam(4)
155 nu1 = nu/(one - nu)
156 alfi = one-alfa
157 dtinv = one /
max(em20, timestep)
158
159
160 IF (time == zero) THEN
161 nup(1:nel) = nupc
162 uvar(1:nel,7) = nupc
163 ELSE
164 nup(1:nel) = uvar(1:nel,7)
165 ENDIF
166
167
168 epspt(1:nel) = uvar(1:nel,1)
169 epspc(1:nel) = uvar(1:nel,2)
170 epsps(1:nel) = uvar(1:nel,3)
171 epdt_f(1:nel) = uvar(1:nel,4)
172 epdc_f(1:nel) = uvar(1:nel,5)
173 epds_f(1:nel) = uvar(1:nel,6)
174 dplat(1:nel) = zero
175 dplac(1:nel) = zero
176 dplas(1:nel) = zero
177 DO i=1,nel
178 epdt(i) =
min(epdt_max,
max(epdt_f(i),epdt_min)) * xfac
179 epdc(i) =
min(epdc_max,
max(epdc_f(i),epdc_min)) * xfac
180 epds(i) =
min(epds_max,
max(epds_f(i),epds_min)) * xfac
181 ENDDO
182
183
184 xvec(1:nel,1) = epspt(1:nel)
185 xvec(1:nel,2) = epdt(1:nel)
187 sigt(1:nel) = sigt(1:nel) * tfac(1)
188 dft(1:nel) = dft(1:nel) * tfac(1)
190 xvec(1:nel,1) = epspc(1:nel)
191 xvec(1:nel,2) = epdc(1:nel)
193 sigc(1:nel) = sigc(1:nel) * tfac(2)
194 dfc(1:nel) = dfc(1:nel) * tfac(2)
195 END IF
196 IF (table(func_shear)%NOTABLE > 0) THEN
197 xvec(1:nel,1) = epsps(1:nel)
198 xvec(1:nel,2) = epds(1:nel)
200 sigs(1:nel) = sigs(1:nel) * tfac(3)
201 dfs(1:nel) = dfs(1:nel) * tfac(3)
202 END IF
203 IF (icas == 0) THEN
204 sigc(1:nel) = sigt(1:nel)
205 sigs(1:nel) = sigt(1:nel)/sqr3
206 ELSEIF (icas == 1) THEN
207 DO i=1,nel
208 sigs(i) = sqrt(sigt(i)*sigc(i)/three)
209 END DO
210 ENDIF
211
212 IF (iconv == 1) THEN
213 DO i = 1,nel
214 conv(i) = .false.
215 IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN
216 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
217 conv(i) = .true.
218 ENDIF
219 ENDDO
220 ENDIF
221 DO i=1,nel
222 aa = one/sigc(i)/sigt(i)
223 a0(i) = three*sigs(i)**2
224 a1(i) = nine*sigs(i)**2*(sigc(i) - sigt(i))*aa
225 a2(i) = nine*(sigc(i)*sigt(i) - three*sigs(i)**2)*aa
226 END DO
227
228
229
230
231 DO i=1,nel
232
233 signxx(i) = sigoxx(i) + a11_2d(i)*depsxx(i) + a12_2d(i)*depsyy(i)
234 signyy(i) = sigoyy(i) + a11_2d(i)*depsyy(i) + a12_2d(i)*depsxx(i)
235 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
236 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
237 signzx(i) = sigozx(i) + depszx(i)*gs(i)
238 p(i) = -(signxx(i) + signyy(i)) * third
239
240 sxx(i) = signxx(i) + p(i)
241 syy(i) = signyy(i) + p(i)
242 szz(i) = p(i)
243 dezz(i)= -nu1 * (depsxx(i) + depsyy(i))
244 soundsp(i) = sqrt(a11_2d(i)/rho(i))
245 ENDDO
246
247
248
249
250
251 DO i = 1,nel
252 alpha(i) = (nine/two)*((one-two*nup(i))/(one+nup(i)))
253 ENDDO
254 niter = 4
255 nindx = 0
256 DO i=1,nel
257 svm2(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*signxy(i)**2
258 svm(i) = sqrt(svm2(i))
259 ylds(i) = svm2(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
260 IF (ylds(i) > 0 .AND. off(i) == one) THEN
261 nindx = nindx + 1
262 indx(nindx) = i
263 ENDIF
264 ENDDO
265
266
267
268
269 IF (nindx > 0) THEN
270
271
272 DO iter = 1,niter
273
274
275 DO ii = 1,nindx
276 i = indx(ii)
277
278
279
280
281
282
283
284
285
286
287
288 gf(i) = sqrt(
max(em20,(svm(i)**2)+
alpha(i)*p(i)**2))
289
290
291 normgxx = (three_half*sxx(i) - third*
alpha(i)*p(i) ) /gf(i)
292 normgyy = (three_half*syy(i) - third*
alpha(i)*p(i) ) /gf(i)
293 normgzz = (three_half*szz(i) - third*
alpha(i)*p(i) ) /gf(i)
294 normgxy = three*signxy(i)/gf(i)
295
296
297 cb = a1(i) + two*a2(i)*p(i)
298 normxx = three * sxx(i) + cb /three
299 normyy = three * syy(i) + cb /three
300 normzz = three * szz(i) + cb /three
301 normxy = two *three * signxy(i)
302
303
304 dfdsigdlam = normxx * (a11_2d(i)*normgxx + a12_2d(i)*normgyy)
305 . + normyy * (a11_2d(i)*normgyy + a12_2d(i)*normgxx)
306 . + normxy * normgxy * g(i)
307
308 yld_norm = svm(i)/gf(i)
309 bb = three_half/(one + nup(i))
310 dft(i) = dft(i) * yld_norm * bb
311 IF (table(
func_comp)%NOTABLE > 0) dfc(i) = dfc(i) * yld_norm * bb
312 IF (table(func_shear)%NOTABLE > 0) dfs(i) = dfs(i) * yld_norm * sqr3/two
313 IF (icas == 0) THEN
314 dfc(i) = dft(i)
315 dfs(i) = (one/sqr3)*dft(i)
316 ELSEIF (icas == 1) THEN
317 dfs(i) = (one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
318 . (dfc(i)*sigt(i) + sigc(i)*dft(i))
319 ENDIF
320 IF (iconv == 1) THEN
321 IF (conv(i)) THEN
322 dfs(i) = sfac*(one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
323 . (dfc(i)*sigt(i) + sigc(i)*dft(i))
324 ENDIF
325 ENDIF
326
327
328 cc = sigs(i)/sigc(i)/sigt(i)
329
330 a1s = eighteen*(sigc(i) - sigt(i))*cc
331 a1c = nine*(sigs(i)/sigc(i))**2
332 a1t = -nine*(sigs(i)/sigt(i))**2
333
334 a2s = -cinquante4*cc
335 a2c = twenty7*cc*sigs(i)/sigc(i)
336 a2t = twenty7*cc*sigs(i)/sigt(i)
337
338 da0 = six*sigs(i)*dfs(i)
339 da1 = a1s*dfs(i) + a1t*dft(i) + a1c*dfc(i)
340 da2 = a2s*dfs(i) + a2t*dft(i) + a2c*dfc(i)
341
342 ff(i) = dfdsigdlam + da0 + p(i)*da1 + p(i)**2 * da2
343 ff(i) = sign(
max(abs(ff(i)),em20) ,ff(i))
344
345 dlam = ylds(i)/ff(i)
346
347 dpxx(i) = dlam * normgxx
348 dpyy(i) = dlam * normgyy
349 dpzz(i) = dlam * normgzz
350 dpxy(i) = dlam * normgxy
351
352
353 signxx(i) = signxx(i) - (a11_2d(i)*dpxx(i) + a12_2d(i)*dpyy(i))
354 signyy(i) = signyy(i) - (a11_2d(i)*dpyy(i) + a12_2d(i)*dpxx(i))
355 signxy(i) = signxy(i) - dpxy(i)*g(i)
356
357
358 epspt(i) = epspt(i) + dlam* yld_norm * bb
359 epspc(i) = epspt(i)
360 epsps(i) = epsps(i) + dlam* yld_norm*sqr3/two
361
362 pla(i) = pla(i) + dlam*yld_norm*off(i)
363 dpla(i) = dpla(i) + dlam *yld_norm
364 dplat(i) = dplat(i) + dlam* yld_norm*bb
365 dplac(i) = dplat(i)
366 dplas(i) = dplas(i) + dlam* yld_norm*sqr3/two
367
368 ENDDO
369
370
371 xvec(1:nel,1) = epspt(1:nel)
372 xvec(1:nel,2) = epdt(1:nel)
374 sigt(1:nel) = sigt(1:nel) * tfac(1)
375 dft(1:nel) = dft(1:nel) * tfac(1)
377 xvec(1:nel,1) = epspc(1:nel)
378 xvec(1:nel,2) = epdc(1:nel)
380 sigc(1:nel) = sigc(1:nel) * tfac(2)
381 dfc(1:nel) = dfc(1:nel) * tfac(2)
382 END IF
383 IF (table(func_shear)%NOTABLE > 0) THEN
384 xvec(1:nel,1) = epsps(1:nel)
385 xvec(1:nel,2) = epds(1:nel)
387 sigs(1:nel) = sigs(1:nel) * tfac(3)
388 dfs(1:nel) = dfs(1:nel) * tfac(3)
389 END IF
390 IF (icas == 0) THEN
391 DO ii = 1,nindx
392 i = indx(ii)
393 sigc(i) = sigt(i)
394 sigs(i) = sigt(i)/sqr3
395 ENDDO
396 ELSEIF (icas == 1) THEN
397 DO ii = 1,nindx
398 i = indx(ii)
399 sigs(i) = sqrt(sigt(i)*sigc(i)/three)
400 ENDDO
401 ENDIF
402 IF (iconv == 1) THEN
403 DO ii = 1,nindx
404 i = indx(ii)
405 conv(i) = .false.
406 IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN
407 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
408 conv(i) = .true.
409 ENDIF
410 END DO
411 ENDIF
412 DO ii = 1,nindx
413 i = indx(ii)
414 aa = one/sigc(i)/sigt(i)
415 a0(i) = three*sigs(i)**2
416 a1(i) = nine*sigs(i)**2*(sigc(i) - sigt(i))*aa
417 a2(i) = nine*(sigc(i)*sigt(i) - three*sigs(i)**2)*aa
418 END DO
419
420 ! update yield function value
421 DO ii = 1,nindx
422 i = indx(ii)
423 p(i) = -third*(signxx(i) + signyy(i) )
424 sxx(i) = signxx(i) + p(i)
425 syy(i) = signyy(i) + p(i)
426 szz(i) = p(i)
427 svm2(i)= three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*signxy(i)**2
428 svm(i) = sqrt(svm2(i))
429 ylds(i) = svm2(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
430 IF (inloc == 0) THEN
431 dezz(i) = dezz(i) + nu1*(dpxx(i) + dpyy(i)) + dpzz(i)
432 ENDIF
433 ENDDO
434 ENDDO
435 ENDIF
436
437
438
439
440 IF (ifunc(1) > 0) THEN
441 iad(1:nel) = npf(ifunc(1)) / 2 + 1
442 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - vartmp(1:nel,8)
443
444 CALL vinter(tf,iad,vartmp(1:nel,8),ilen,nel,pla,dydx,nup)
445
446 uvar(1:nel,7) = yfac(1) * nup(1:nel)
447 uvar(1:nel,7) =
max(zero,
min(nup(1:nel), half))
448 END IF
449
450
451
452
453 IF (inloc > 0) THEN
454 DO i = 1,nel
455 IF (loff(i) == one) THEN
456 alpha(i)= (nine/two)*((one-two*nup(i))/(one+nup(i)))
457 gf(i) = sqrt(
max(em20,(svm(i)**2)+
alpha(i)*p(i)**2))
458 normgxx = (three_half*sxx(i) - third*
alpha(i)*p(i))/gf(i)
459 normgyy = (three_half*syy(i) - third*
alpha(i)*p(i))/gf(i)
460 normgzz = (three_half*szz(i) - third*
alpha(i)*p(i))/gf(i)
461 yld_norm = svm(i)/gf(i)
462 IF (yld_norm /= zero) THEN
463 dlam_nl(i) = (one/yld_norm)*
max(dplanl(i),zero)
464 dezz(i) = dezz(i) + nu1*(dlam_nl(i)*normgxx)
465 . + nu1*(dlam_nl(i)*normgyy)
466 . + dlam_nl(i)*normgzz
467 ENDIF
468 ENDIF
469 ENDDO
470 ENDIF
471
472
473
474
475 uvar(1:nel,1) = epspt(1:nel)
476 uvar(1:nel,2) = epspc(1:nel)
477 uvar(1:nel,3) = epsps(1:nel)
478 uvar(1:nel,4) = alfa*dplat(1:nel)*dtinv + alfi*epdt_f(1:nel)
479 uvar(1:nel,5) = alfa*dplac(1:nel)*dtinv + alfi*epdc_f(1:nel)
480 uvar(1:nel,6) = alfa*dplas(1:nel)*dtinv + alfi*epds_f(1
481 epsd(1:nel) = alfa*dpla(1:nel)*dtinv + alfi*epsd(1:nel)
482
subroutine func_comp(table, ntable, nptmax, npt_trac, npt_shear, npt_comp, x_comp, y_comp, nup)
subroutine table_mat_vinterp(table, dimx, nel, ipos, xx, yy, dydx)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)