OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps110c_nice.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sigeps110c_nice ../engine/source/materials/mat/mat110/sigeps110c_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps110c ../engine/source/materials/mat/mat110/sigeps110c.F
27!||--- calls -----------------------------------------------------
28!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.F
29!|| table_vinterp ../engine/source/tools/curve/table_tools.F
30!||--- uses -----------------------------------------------------
31!|| interface_table_mod ../engine/share/modules/table_mod.F
32!|| table_mod ../engine/share/modules/table_mod.F
33!||====================================================================
34 SUBROUTINE sigeps110c_nice(
35 1 NEL ,NGL ,NUPARAM ,NUVAR ,NPF ,
36 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
37 3 GS ,RHO ,PLA ,DPLA ,EPSP ,SOUNDSP ,
38 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,ASRATE ,
39 5 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
40 6 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 7 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,THKLY ,
42 8 THK ,SIGY ,ET ,TEMPEL ,TEMP ,SEQ ,
43 9 TF ,NUMTABL ,ITABLE ,TABLE ,NVARTMP ,VARTMP ,
44 A SIGA ,INLOC ,DPLANL ,LOFF )
45 !=======================================================================
46 ! Modules
47 !=======================================================================
48 USE table_mod
49 USE interface_table_mod
50 !=======================================================================
51 ! Implicit types
52 !=======================================================================
53#include "implicit_f.inc"
54 !=======================================================================
55 ! Common
56 !=======================================================================
57#include "com04_c.inc"
58 !=======================================================================
59 ! Dummy arguments
60 !=======================================================================
61 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP,NPF(*),INLOC
62 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
63 my_real
64 . TIME,TIMESTEP,ASRATE,TF(*)
65 INTEGER :: VARTMP(NEL,NVARTMP)
66 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
67 . UPARAM
68 my_real,DIMENSION(NEL), INTENT(IN) ::
69 . RHO,TEMPEL,
70 . DEPSXX,DEPSYY,DEPSXY,DEPSYZ,DEPSZX,
71 . epspxx,epspyy,epspxy,epspyz,epspzx ,
72 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
73 . gs,thkly,dplanl,loff
74c
75 my_real ,DIMENSION(NEL), INTENT(OUT) ::
76 . soundsp,sigy,et,epsp,
77 . signxx,signyy,signxy,signyz,signzx
78c
79 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
80 . pla,dpla,off,thk,temp,seq
81 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
82 . uvar
83 my_real ,DIMENSION(NEL,3) ,INTENT(INOUT) ::
84 . siga
85c
86 TYPE(ttable), DIMENSION(NTABLE) :: TABLE
87 !=======================================================================
88 ! Local Variables
89 !=======================================================================
90 INTEGER I,J,II,K,ITER,NITER,NINDX,INDEX(NEL),VFLAG,IPOS(NEL,2),NANGLE,
91 . ipos0(nel,2),ismooth
92c
93 my_real
94 . young,g,g2,nu,nnu,a11,a12,yld0,dsigm,beta,omega,nexp,eps0,sigst,
95 . dg0,deps0,mexp,fbi(2),kboltz,fisokin,tini,xscale,yscale
96 my_real
97 . mohr_radius,mohr_center,normsig,tmp,sig_ratio,var_a,var_b,var_c,
98 . a(2),b(2),c(2),h,dpdt,dlam,ddep,dav,deve1,deve2,deve3,deve4,
99 . xvec(nel,4)
100 my_real
101 . fun_theta(nel,2),fsh_theta(nel,2),fps_theta(nel,2),
102 . hips_theta(nel,2),hiun_theta(nel,2),hish_theta(nel,2)
103 my_real
104 . f1,f2,df1_dmu,df2_dmu,normxx,normyy,normxy,
105 . denom,sig_dfdsig,dfdsig2,dphi_dsig1,dphi_dsig2,
106 . dsxx,dsyy,dsxy,dexx,deyy,dexy,alpha,da_dcos2(2),
107 . db_dcos2(2),dc_dcos2(2),df1_dcos2,df2_dcos2,
108 . dphi_dcos2,cos2(10,10),dphi
109c
110 my_real, DIMENSION(NEL) ::
111 . dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,sigvg,yld,hardp,sighard,sigrate,
112 . phi,dpla_dlam,dezz,dphi_dlam,dpxx,dpyy,dpxy,dpyz,dpzx,dpzz,
113 . sig1,sig2,cos2theta,sin2theta,deelzz,mu,dyld_dpla,yl0,sigexx,sigeyy,
114 . sigexy,hardr,yld_tref,dydx,yld_temp,tfac,phi0
115c
116 my_real, DIMENSION(:,:), ALLOCATABLE ::
117 . hips,hiun,hish,q_fsh,q_fps,q_hiun,q_hips,q_hish
118c
119 my_real, DIMENSION(:), ALLOCATABLE ::
120 . q_fun
121c
122 my_real, PARAMETER :: tol = 1.0d-6
123c
124 LOGICAL :: SIGN_CHANGED(NEL)
125c
126 DATA cos2/
127 1 1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
128 2 0. ,1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
129 3 -1. ,0. ,2. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
130 4 0. ,-3. ,0. ,4. ,0. ,0. ,0. ,0. ,0. ,0. ,
131 5 1. ,0. ,-8. ,0. ,8. ,0. ,0. ,0. ,0. ,0. ,
132 6 0. ,5. ,0. ,-20. ,0. ,16. ,0. ,0. ,0. ,0. ,
133 7 -1. ,0. ,18. ,0. ,-48. ,0. ,32. ,0. ,0. ,0. ,
134 8 0. ,-7. ,0. ,56. ,0. ,-112.,0. ,64. ,0. ,0. ,
135 9 1. ,0. ,-32. ,0. ,160. ,0. ,-256.,0. ,128. ,0. ,
136 a 0. ,9. ,0. ,-120.,0. ,432. ,0. ,-576 ,0. ,256. /
137 !=======================================================================
138 ! DRUCKER MATERIAL LAW WITH NO DAMAGE
139 !=======================================================================
140 !UVAR(1)
141 !UVAR(2) YLD YIELD STRESS
142 !UVAR(3) H HARDENING RATE
143 !UVAR(4) PHI YIELD FUNCTION VALUE
144 !-
145 !DEPIJ PLASTIC STRAIN TENSOR COMPONENT
146 !DEPSIJ TOTAL STRAIN TENSOR COMPONENT (EL+PL)
147 !=======================================================================
148c
149 !=======================================================================
150 ! - INITIALISATION OF COMPUTATION ON TIME STEP
151 !=======================================================================
152 ! Recovering model parameters
153 young = uparam(1)
154 nu = uparam(2)
155 g = uparam(3)
156 g2 = uparam(4)
157 nnu = uparam(5)
158 a11 = uparam(6)
159 a12 = uparam(7)
160 yld0 = uparam(8)
161 dsigm = uparam(9)
162 beta = uparam(10)
163 omega = uparam(11)
164 nexp = uparam(12)
165 eps0 = uparam(13)
166 sigst = uparam(14)
167 dg0 = uparam(15)
168 deps0 = uparam(16)
169 mexp = uparam(17)
170 kboltz = uparam(18)
171 xscale = uparam(19)
172 yscale = uparam(20)
173 fisokin = uparam(21)
174 vflag = nint(uparam(23))
175 tini = uparam(24)
176 fbi(1) = uparam(25)
177 fbi(2) = uparam(25)
178 nangle = nint(uparam(26))
179 ismooth = nint(uparam(28))
180c
181 ALLOCATE(hips(nangle,2),hiun(nangle,2),hish(nangle,2),
182 . q_fsh(nangle,2),q_fun(nangle),q_fps(nangle,2),
183 . q_hish(nangle,2),q_hiun(nangle,2),q_hips(nangle,2))
184c
185 DO i = 1, nangle
186 ! Hinge points
187 hips(i,1) = uparam(30+17*(i-1))
188 hips(i,2) = uparam(31+17*(i-1))
189 hiun(i,1) = uparam(32+17*(i-1))
190 hiun(i,2) = uparam(33+17*(i-1))
191 hish(i,1) = uparam(34+17*(i-1))
192 hish(i,2) = uparam(35+17*(i-1))
193 ! Interpolation factors
194 q_fsh(i,1) = uparam(36+17*(i-1))
195 q_fsh(i,2) = uparam(37+17*(i-1))
196 q_fun(i) = uparam(38+17*(i-1))
197 q_fps(i,1) = uparam(39+17*(i-1))
198 q_fps(i,2) = uparam(40+17*(i-1))
199 q_hish(i,1) = uparam(41+17*(i-1))
200 q_hish(i,2) = uparam(42+17*(i-1))
201 q_hiun(i,1) = uparam(43+17*(i-1))
202 q_hiun(i,2) = uparam(44+17*(i-1))
203 q_hips(i,1) = uparam(45+17*(i-1))
204 q_hips(i,2) = uparam(46+17*(i-1))
205 ENDDO
206c
207 ! Initial variable
208 IF (time == zero) THEN
209 IF (jthe == 0) THEN
210 temp(1:nel) = tini
211 ENDIF
212 IF (eps0 > zero) THEN
213 pla(1:nel) = eps0
214 ENDIF
215 ENDIF
216c
217 ! Recovering internal variables
218 DO i=1,nel
219 IF (off(i) < 0.1) off(i) = zero
220 IF (off(i) < one) off(i) = off(i)*four_over_5
221 ! User inputs
222 phi0(i) = uvar(i,2) ! Previous yield function value
223 ! Standard inputs
224 dpla(i) = zero ! Initialization of the plastic strain increment
225 et(i) = one ! Initialization of hourglass coefficient
226 hardp(i) = zero ! Initialization of hourglass coefficient
227 sign_changed(i) = .false.
228 dezz(i) = zero
229 ENDDO
230c
231 ! /HEAT/MAT temperature
232 IF (jthe > 0) THEN
233 temp(1:nel) = tempel(1:nel)
234 ENDIF
235c
236 ! Computation of the strain-rate depending on the flags
237 IF ((vflag == 1) .OR. (vflag == 3)) THEN
238 DO i = 1, nel
239 IF (vflag == 1) THEN
240 epsp(i) = uvar(i,1)
241 ELSEIF (vflag == 3) THEN
242 dav = (epspxx(i)+epspyy(i))*third
243 deve1 = epspxx(i) - dav
244 deve2 = epspyy(i) - dav
245 deve3 = - dav
246 deve4 = half*epspxy(i)
247 epsp(i) = half*(deve1**2 + deve2**2 + deve3**2) + deve4**2
248 epsp(i) = sqrt(three*epsp(i))/three_half
249 epsp(i) = asrate*epsp(i) + (one - asrate)*uvar(i,1)
250 uvar(i,1) = epsp(i)
251 ENDIF
252 ENDDO
253 ENDIF
254c
255 ! Initial yield stress for kinematic hardening
256 IF (fisokin > zero) THEN
257 IF (numtabl > 0) THEN
258 xvec(1:nel,1) = zero
259 xvec(1:nel,2) = epsp(1:nel) * xscale
260 ipos0(1:nel,1) = 1
261 ipos0(1:nel,2) = 1
262 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos0,xvec ,yl0 ,hardp,hardr)
263 yl0(1:nel) = yl0(1:nel) * yscale
264 ! Tabulation with Temperature
265 IF (numtabl == 2) THEN
266 ! Reference temperature factor
267 xvec(1:nel,2) = tini
268 ipos0(1:nel,1) = 1
269 ipos0(1:nel,2) = 1
270 CALL table_vinterp(table(itable(2)),nel,nel,ipos0,xvec,yld_tref,dydx)
271 ! Current temperature factor
272 xvec(1:nel,2) = temp(1:nel)
273 CALL table_vinterp(table(itable(2)),nel,nel,ipos0,xvec,yld_temp,dydx)
274 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
275 yl0(1:nel) = yl0(1:nel) * tfac(1:nel)
276 ENDIF
277 ELSE
278 yl0(1:nel) = yld0
279 ENDIF
280 ELSE
281 yl0(1:nel) = zero
282 ENDIF
283c
284 ! Computation of the yield stress
285 IF (numtabl > 0) THEN
286 ! Tabulation with strain-rate
287 xvec(1:nel,1) = pla(1:nel)
288 xvec(1:nel,2) = epsp(1:nel) * xscale
289 ipos(1:nel,1) = vartmp(1:nel,1)
290 ipos(1:nel,2) = 1
291 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,yld,hardp,hardr)
292 yld(1:nel) = yld(1:nel) * yscale
293 hardp(1:nel) = hardp(1:nel) * yscale
294 vartmp(1:nel,1) = ipos(1:nel,1)
295 ! Tabulation with Temperature
296 IF (numtabl == 2) THEN
297 ! Reference temperature factor
298 xvec(1:nel,2) = tini
299 ipos(1:nel,1) = vartmp(1:nel,2)
300 ipos(1:nel,2) = vartmp(1:nel,3)
301 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_tref,dydx)
302 vartmp(1:nel,2) = ipos(1:nel,1)
303 vartmp(1:nel,3) = ipos(1:nel,2)
304 ! Current temperature factor
305 xvec(1:nel,2) = temp(1:nel)
306 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_temp,dydx)
307 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
308 yld(1:nel) = yld(1:nel) * tfac(1:nel)
309 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
310 ELSE
311 tfac(1:nel) = one
312 ENDIF
313 ! Including kinematic hardening effect
314 yld(1:nel) = (one-fisokin)*yld(1:nel)+fisokin*yl0(1:nel)
315 ELSE
316 DO i = 1,nel
317 ! a) - Hardening law
318 ! Continuous law
319 sighard(i) = yld0 + dsigm*(beta*(pla(i))+(one-exp(-omega*(pla(i))))**nexp)
320 ! b) - Strain-rate dependent hardening stress
321 sigrate(i) = sigst*(one + (kboltz*temp(i)/dg0)*log(one + (epsp(i)/deps0)))**mexp
322 ! c) - Computation of the yield function and value check
323 yld(i) = sighard(i) + sigrate(i)
324 ! d) - Including kinematic hardening
325 IF (fisokin > zero) THEN
326 yl0(i) = yl0(i) + sigrate(i)
327 ENDIF
328 yld(i) = (one-fisokin)*yld(i)+fisokin*yl0(i)
329 ! d) - Checking values
330 yld(i) = max(em10, yld(i))
331 ENDDO
332 ENDIF
333c
334 !========================================================================
335 ! - COMPUTATION OF TRIAL VALUES
336 !========================================================================
337 DO i=1,nel
338c
339 ! computation of the trial stress tensor
340 IF (fisokin > zero) THEN
341 signxx(i) = sigoxx(i) - siga(i,1) + a11*depsxx(i) + a12*depsyy(i)
342 signyy(i) = sigoyy(i) - siga(i,2) + a12*depsxx(i) + a11*depsyy(i)
343 signxy(i) = sigoxy(i) - siga(i,3) + g*depsxy(i)
344 sigexx(i) = signxx(i)
345 sigeyy(i) = signyy(i)
346 sigexy(i) = signxy(i)
347 ELSE
348 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
349 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
350 signxy(i) = sigoxy(i) + depsxy(i)*g
351 ENDIF
352 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
353 signzx(i) = sigozx(i) + depszx(i)*gs(i)
354c
355 ! Computation of the equivalent stress
356 normsig = sqrt(signxx(i)*signxx(i)
357 . + signyy(i)*signyy(i)
358 . + two*signxy(i)*signxy(i))
359 IF (normsig < em20) THEN
360 sigvg(i) = zero
361 ELSE
362c
363 ! Computation of the principal stresses
364 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
365 mohr_center = (signxx(i)+signyy(i))/two
366 sig1(i) = mohr_center + mohr_radius
367 sig2(i) = mohr_center - mohr_radius
368 IF (mohr_radius>em20) THEN
369 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
370 sin2theta(i) = signxy(i)/mohr_radius
371 ELSE
372 cos2theta(i) = one
373 sin2theta(i) = zero
374 ENDIF
375c
376 ! Computation of scale factors for shear
377 fsh_theta(i,1:2) = zero
378 fps_theta(i,1:2) = zero
379 DO j = 1,nangle
380 DO k = 1,j
381 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
382 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
383 ENDDO
384 ENDDO
385c
386 ! Computation of the equivalent stress of Vegter
387 IF (sig1(i)<zero .OR. ((sig2(i)<zero) .AND. (sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2)))) THEN
388 cos2theta(i) = -cos2theta(i)
389 sin2theta(i) = -sin2theta(i)
390 tmp = sig1(i)
391 sig1(i) = -sig2(i)
392 sig2(i) = -tmp
393 sign_changed(i) = .true.
394 ! Computation of scale factors
395 fsh_theta(i,1:2) = zero
396 fps_theta(i,1:2) = zero
397 DO j = 1,nangle
398 DO k = 1,j
399 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
400 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
401 ENDDO
402 ENDDO
403 ELSE
404 sign_changed(i) = .false.
405 ENDIF
406 ! Between shear and uniaxial point
407 IF (sig2(i)<zero) THEN
408 ! Interpolation of factors
409 fun_theta(i,1:2) = zero
410 hish_theta(i,1:2) = zero
411 DO j = 1,nangle
412 DO k = 1,j
413 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
414 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
415 ENDDO
416 ENDDO
417 ! Filling the table
418 a(1:2) = fsh_theta(i,1:2)
419 b(1:2) = hish_theta(i,1:2)
420 c(1:2) = fun_theta(i,1:2)
421 ! Between plane strain and uniaxial point
422 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
423 ! Interpolation of factors
424 fun_theta(i,1:2) = zero
425 hiun_theta(i,1:2) = zero
426 DO j = 1,nangle
427 DO k = 1,j
428 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
429 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
430 ENDDO
431 ENDDO
432 ! Filling the table
433 a(1:2) = fun_theta(i,1:2)
434 b(1:2) = hiun_theta(i,1:2)
435 c(1:2) = fps_theta(i,1:2)
436 ! Between plane strain and equibiaxial point
437 ELSE
438 ! Interpolation of factors
439 hips_theta(i,1:2) = zero
440 DO j = 1,nangle
441 DO k = 1,j
442 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
443 ENDDO
444 ENDDO
445 ! Filling the table
446 a(1:2) = fps_theta(i,1:2)
447 b(1:2) = hips_theta(i,1:2)
448 c(1:2) = fbi(1:2)
449 ENDIF
450 ! Determine MU and SIGVG
451 IF (sig1(i)<em20) THEN
452 mu(i) = zero
453 sigvg(i) = zero
454 ELSE
455 sig_ratio = sig2(i)/sig1(i)
456 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
457 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
458 var_c = a(2) - sig_ratio*a(1)
459 IF (abs(var_a)<em08) THEN
460 mu(i) = -var_c/var_b
461 ELSE
462 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
463 ENDIF
464 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
465 END IF
466 ENDIF
467c
468 ENDDO
469c
470 !========================================================================
471 ! - COMPUTATION OF YIELD FONCTION
472 !========================================================================
473 phi(1:nel) = sigvg(1:nel) - yld(1:nel)
474c
475 ! Checking plastic behavior for all elements
476 nindx = 0
477 DO i=1,nel
478 IF (phi(i) > zero .AND. off(i) == one) THEN
479 nindx=nindx+1
480 index(nindx)=i
481 ENDIF
482 ENDDO
483c
484 !=======================================================
485 ! - PLASTIC CORRECTION WITH NICE EXPLICIT RETURN MAPPING
486 !=======================================================
487c
488 ! Loop over yielding elements
489#include "vectorize.inc"
490 DO ii=1,nindx
491c
492 ! Number of the element with plastic behaviour
493 i = index(ii)
494c
495 ! Computation of the trial stress increment
496 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
497 dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
498 dsigxy(i) = depsxy(i)*g
499 dsigyz(i) = depsyz(i)*gs(i)
500 dsigzx(i) = depszx(i)*gs(i)
501 dlam = zero
502c
503 ! Computation of the old principal stresses
504 mohr_radius = sqrt(((sigoxx(i)-sigoyy(i))/two)**2 + sigoxy(i)**2)
505 mohr_center = (sigoxx(i)+sigoyy(i))/two
506 sig1(i) = mohr_center + mohr_radius
507 sig2(i) = mohr_center - mohr_radius
508 IF (mohr_radius>em20) THEN
509 cos2theta(i) = ((sigoxx(i)-sigoyy(i))/two)/mohr_radius
510 sin2theta(i) = sigoxy(i)/mohr_radius
511 ELSE
512 cos2theta(i) = one
513 sin2theta(i) = zero
514 ENDIF
515c
516 ! Computation of scale factors for shear
517 fsh_theta(i,1:2) = zero
518 fps_theta(i,1:2) = zero
519 DO j = 1,nangle
520 DO k = 1,j
521 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
522 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
523 ENDDO
524 ENDDO
525c
526 ! computation of the equivalent stress of vegter
527 IF (sig1(i)<zero .OR. ((sig2(i)<zero) .AND. (sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2)))) THEN
528 cos2theta(i) = -cos2theta(i)
529 sin2theta(i) = -sin2theta(i)
530 tmp = sig1(i)
531 sig1(i) = -sig2(i)
532 sig2(i) = -tmp
533 sign_changed(i) = .true.
534 ! Computation of scale factors
535 fsh_theta(i,1:2) = zero
536 fps_theta(i,1:2) = zero
537 DO j = 1,nangle
538 DO k = 1,j
539 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
540 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
541 ENDDO
542 ENDDO
543 ELSE
544 sign_changed(i) = .false.
545 ENDIF
546 ! Between shear and uniaxial point
547 IF (sig2(i)<zero) THEN
548 ! Interpolation of factors
549 fun_theta(i,1:2) = zero
550 hish_theta(i,1:2) = zero
551 DO j = 1,nangle
552 DO k = 1,j
553 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
554 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
555 ENDDO
556 ENDDO
557 ! Filling the table
558 a(1:2) = fsh_theta(i,1:2)
559 b(1:2) = hish_theta(i,1:2)
560 c(1:2) = fun_theta(i,1:2)
561 ! Derivatives of PHI with respect to THETA
562 da_dcos2(1:2) = zero
563 db_dcos2(1:2) = zero
564 dc_dcos2(1:2) = zero
565 ! If anisotropic, compute derivatives with respect to COS2THET
566 IF (nangle > 1) THEN
567 ! Computation of their first derivative with respect to COS2THET
568 DO j = 2, nangle
569 DO k = 2, j
570 da_dcos2(1:2) = da_dcos2(1:2) + q_fsh(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
571 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
572 dc_dcos2(1) = dc_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
573 ENDDO
574 ENDDO
575 ENDIF
576 ! Between plane strain and uniaxial point
577 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
578 ! Interpolation of factors
579 fun_theta(i,1:2) = zero
580 hiun_theta(i,1:2) = zero
581 DO j = 1,nangle
582 DO k = 1,j
583 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
584 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
585 ENDDO
586 ENDDO
587 ! Filling the table
588 a(1:2) = fun_theta(i,1:2)
589 b(1:2) = hiun_theta(i,1:2)
590 c(1:2) = fps_theta(i,1:2)
591 ! Derivatives of PHI with respect to THETA
592 da_dcos2(1:2) = zero
593 db_dcos2(1:2) = zero
594 dc_dcos2(1:2) = zero
595 ! If anisotropic, compute derivatives with respect to COS2THET
596 IF (nangle > 1) THEN
597 ! Computation of their first derivative with respect to COS2THET
598 DO j = 2, nangle
599 DO k = 2,j
600 da_dcos2(1) = da_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
601 db_dcos2(1:2) = db_dcos2(1:2) + q_hiun(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
602 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
603 ENDDO
604 ENDDO
605 ENDIF
606 ! Between plane strain and equibiaxial point
607 ELSE
608 ! Interpolation of factors
609 hips_theta(i,1:2) = zero
610 DO j = 1,nangle
611 DO k = 1,j
612 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
613 ENDDO
614 ENDDO
615 ! Filling the table
616 a(1:2) = fps_theta(i,1:2)
617 b(1:2) = hips_theta(i,1:2)
618 c(1:2) = fbi(1:2)
619 ! derivatives of phi with respect to theta
620 da_dcos2(1:2) = zero
621 db_dcos2(1:2) = zero
622 dc_dcos2(1:2) = zero
623 ! If anisotropic, compute derivatives with respect to COS2THET
624 IF (nangle > 1) THEN
625 ! Computation of their first derivative with respect to COS2THET
626 DO j = 2, nangle
627 DO k = 2,j
628 da_dcos2(1:2) = da_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
629 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
630 ENDDO
631 ENDDO
632 ENDIF
633 ENDIF
634 ! Determine MU and SIGVG
635 IF (sig1(i)<em20) THEN
636 mu(i) = zero
637 sigvg(i) = zero
638 ELSE
639 sig_ratio = sig2(i)/sig1(i)
640 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
641 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
642 var_c = a(2) - sig_ratio*a(1)
643 IF (abs(var_a)<em08) THEN
644 mu(i) = -var_c/var_b
645 ELSE
646 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
647 ENDIF
648 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
649 END IF
650c
651 ! Note : in this part, the purpose is to compute for each iteration
652 ! a plastic multiplier allowing to update internal variables to satisfy
653 ! the consistency condition using the backward Euler implicit method
654 ! with a Newton-Raphson iterative procedure
655 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
656 ! -> PHI : current value of yield function (known)
657 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
658 ! into account of internal variables kinetic :
659 ! plasticity, temperature and damage (to compute)
660c
661 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
662 !-------------------------------------------------------------
663c
664 ! Derivatives of PHI with respect to COS2THETA
665 df1_dcos2 = da_dcos2(1) + two*mu(i)*(db_dcos2(1)-da_dcos2(1)) +
666 . mu(i)*mu(i)*(da_dcos2(1)+dc_dcos2(1)-two*db_dcos2(1))
667 df2_dcos2 = da_dcos2(2) + two*mu(i)*(db_dcos2(2)-da_dcos2(2)) +
668 . mu(i)*mu(i)*(da_dcos2(2)+dc_dcos2(2)-two*db_dcos2(2))
669c
670 ! Computation of the derivatives
671 f1 = mu(i)*mu(i)*(a(1)+c(1)-two*b(1))+two*mu(i)*(b(1)-a(1))+a(1)
672 f2 = mu(i)*mu(i)*(a(2)+c(2)-two*b(2))+two*mu(i)*(b(2)-a(2))+a(2)
673c
674 ! Derivatives of Phi with respect to MU
675 df1_dmu = two*(b(1)-a(1)) + two*mu(i)*(a(1)+c(1)-two*b(1))
676 df2_dmu = two*(b(2)-a(2)) + two*mu(i)*(a(2)+c(2)-two*b(2))
677c
678 ! Derivatives of PHI with respect to SIG1, SIG2 and MU
679 IF ((f1*df2_dmu - f2*df1_dmu)/=zero) THEN
680 dphi_dsig1 = df2_dmu/(f1*df2_dmu - f2*df1_dmu)
681 dphi_dsig2 = -df1_dmu/(f1*df2_dmu - f2*df1_dmu)
682 IF (abs(sig1(i)-sig2(i))>tol) THEN
683 dphi_dcos2 = (sigvg(i)/(sig1(i) - sig2(i)))*(df2_dcos2*df1_dmu -
684 . df1_dcos2*df2_dmu)/(f1*df2_dmu - f2*df1_dmu)
685 ELSE
686 dphi_dcos2 = (two*((one-mu(i))*da_dcos2(2)+two*mu(i)*db_dcos2(2))*df1_dmu -
687 . two*((one-mu(i))*da_dcos2(1)+two*mu(i)*db_dcos2(1))*df2_dmu)/
688 . ((one-mu(i))*(a(1)-a(2)) + two*mu(i)*(b(1)-b(2)))
689 ENDIF
690 ELSE
691 dphi_dsig1 = zero
692 dphi_dsig2 = zero
693 dphi_dcos2 = zero
694 ENDIF
695 normxx = half*(one + cos2theta(i))*dphi_dsig1 + half*(one-cos2theta(i))*dphi_dsig2 +
696 . (sin2theta(i)**2)*dphi_dcos2
697 normyy = half*(one - cos2theta(i))*dphi_dsig1 + half*(one+cos2theta(i))*dphi_dsig2 -
698 . (sin2theta(i)**2)*dphi_dcos2
699 normxy = sin2theta(i)*dphi_dsig1 - sin2theta(i)*dphi_dsig2 -
700 . (two*sin2theta(i)*cos2theta(i))*dphi_dcos2
701 IF (sign_changed(i)) THEN
702 normxx = -normxx
703 normyy = -normyy
704 normxy = -normxy
705 ENDIF
706c
707 ! Restoring previous value of the yield function
708 phi(i) = phi0(i)
709c
710 ! Computation of yield surface trial increment DPHI
711 dphi = normxx * dsigxx(i)
712 . + normyy * dsigyy(i)
713 . + normxy * dsigxy(i)
714c
715 ! 2 - Computation of DPHI_DLAMBDA
716 !---------------------------------------------------------
717c
718 ! a) Derivative with respect stress increments tensor DSIG
719 ! --------------------------------------------------------
720 dfdsig2 = normxx * (a11*normxx + a12*normyy)
721 . + normyy * (a11*normyy + a12*normxx)
722 . + normxy * normxy * g
723c
724 ! b) Derivatives with respect to plastic strain P
725 ! ------------------------------------------------
726c
727 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
728 ! ----------------------------------------------------------------------------
729 IF (numtabl == 0) THEN
730 ! Continuous hardening law
731 hardp(i) = dsigm*beta
732 IF (pla(i)>zero) THEN
733 hardp(i) = hardp(i) + dsigm*(nexp*(one-exp(-omega*(pla(i))))**(nexp-one))*
734 . (omega*exp(-omega*(pla(i))))
735 ENDIF
736 ENDIF
737 ! Accounting for kinematic hardening
738 dyld_dpla(i) = (one-fisokin)*hardp(i)
739c
740 ! ii) Derivative of dPLA with respect to DLAM
741 ! -------------------------------------------
742 sig_dfdsig = sigoxx(i) * normxx
743 . + sigoyy(i) * normyy
744 . + sigoxy(i) * normxy
745 dpla_dlam(i) = sig_dfdsig/yld(i)
746c
747 ! 3 - Computation of plastic multiplier and variables update
748 !----------------------------------------------------------
749c
750 ! Derivative of PHI with respect to DLAM
751 dphi_dlam(i) = - dfdsig2 - dyld_dpla(i)*dpla_dlam(i)
752 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
753c
754 ! Computation of the plastic multiplier
755 dlam = - (dphi + phi(i)) / dphi_dlam(i)
756 dlam = max(dlam, zero)
757c
758 ! Plastic strains tensor update
759 dpxx(i) = dlam * normxx
760 dpyy(i) = dlam * normyy
761 dpxy(i) = dlam * normxy
762c
763 ! Elasto-plastic stresses update
764 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
765 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
766 signxy(i) = signxy(i) - dpxy(i)*g
767c
768 ! Cumulated plastic strain and strain rate update
769 ddep = dlam*sig_dfdsig/yld(i)
770 dpla(i) = max(zero, dpla(i) + ddep)
771 pla(i) = pla(i) + ddep
772c
773 ! Computation of the equivalent stress
774 normsig = sqrt(signxx(i)*signxx(i)
775 . + signyy(i)*signyy(i)
776 . + signxy(i)*signxy(i))
777 IF (normsig < em20) THEN
778 sigvg(i) = zero
779 ELSE
780 ! Computation of the principal stresses
781 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
782 mohr_center = (signxx(i)+signyy(i))/two
783 sig1(i) = mohr_center + mohr_radius
784 sig2(i) = mohr_center - mohr_radius
785 IF (mohr_radius>em20) THEN
786 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
787 sin2theta(i) = signxy(i)/mohr_radius
788 ELSE
789 cos2theta(i) = one
790 sin2theta(i) = zero
791 ENDIF
792 ! Computation of scale factors
793 fsh_theta(i,1:2) = zero
794 fps_theta(i,1:2) = zero
795 DO j = 1,nangle
796 DO k = 1,j
797 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
798 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
799 ENDDO
800 ENDDO
801 ! Computation of the equivalent stress of Vegter
802 IF (sig1(i)<zero .OR. (sig2(i)<zero .AND. sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2))) THEN
803 cos2theta(i) = -cos2theta(i)
804 sin2theta(i) = -sin2theta(i)
805 tmp = sig1(i)
806 sig1(i) = -sig2(i)
807 sig2(i) = -tmp
808 sign_changed(i) = .true.
809 fsh_theta(i,1:2) = zero
810 fps_theta(i,1:2) = zero
811 DO j = 1,nangle
812 DO k = 1,j
813 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
814 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
815 ENDDO
816 ENDDO
817 ELSE
818 sign_changed(i) = .false.
819 ENDIF
820 ! Between shear and uniaxial point
821 IF (sig2(i)<zero) THEN
822 ! Interpolation of factors
823 fun_theta(i,1:2) = zero
824 hish_theta(i,1:2) = zero
825 DO j = 1,nangle
826 DO k = 1,j
827 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
828 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
829 ENDDO
830 ENDDO
831 ! Filling the tables
832 a(1:2) = fsh_theta(i,1:2)
833 b(1:2) = hish_theta(i,1:2)
834 c(1:2) = fun_theta(i,1:2)
835 ! Between plane strain and uniaxial point
836 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
837 ! Interpolation of factors
838 fun_theta(i,1:2) = zero
839 hiun_theta(i,1:2) = zero
840 DO j = 1,nangle
841 DO k = 1,j
842 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
843 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
844 ENDDO
845 ENDDO
846 ! Filling the tables
847 a(1:2) = fun_theta(i,1:2)
848 b(1:2) = hiun_theta(i,1:2)
849 c(1:2) = fps_theta(i,1:2)
850 ! Between plane strain and equibiaxial point
851 ELSE
852 ! Interpolation of factors
853 hips_theta(i,1:2) = zero
854 DO j = 1,nangle
855 DO k = 1,j
856 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
857 ENDDO
858 ENDDO
859 ! Filling the table
860 a(1:2) = fps_theta(i,1:2)
861 b(1:2) = hips_theta(i,1:2)
862 c(1:2) = fbi(1:2)
863 ENDIF
864 ! Determine MU and SIGVG
865 IF (sig1(i)<em20) THEN
866 mu(i) = zero
867 sigvg(i) = zero
868 ELSE
869 sig_ratio = sig2(i)/sig1(i)
870 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
871 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
872 var_c = a(2) - sig_ratio*a(1)
873 IF (abs(var_a)<em08) THEN
874 mu(i) = -var_c/var_b
875 ELSE
876 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
877 ENDIF
878 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
879 END IF
880 ENDIF
881c
882 IF (numtabl == 0) THEN
883 ! Continuous hardening law update
884 sighard(i) = yld0 + dsigm*(beta*(pla(i))+(one-exp(-omega*(pla(i))))**nexp)
885 ! Yield stress update
886 yld(i) = sighard(i) + sigrate(i)
887 yld(i) = (one-fisokin)*yld(i)+fisokin*yl0(i)
888 yld(i) = max(yld(i), em10)
889 ! Yield function value update
890 phi(i) = sigvg(i) - yld(i)
891 ENDIF
892c
893 ! Transverse strain update
894 dezz(i) = dezz(i) - (dpxx(i)+dpyy(i))
895c
896 ENDDO
897 ! End of the loop over yielding elements
898 !==============================================================
899 ! - END OF PLASTIC CORRECTION WITH NICE EXPLICIT RETURN MAPPING
900 !==============================================================
901c
902 ! If tabulated yield function, update of the yield stress for all element
903 IF ((numtabl > 0).AND.(nindx > 0)) THEN
904 xvec(1:nel,1) = pla(1:nel)
905 xvec(1:nel,2) = epsp(1:nel) * xscale
906 ipos(1:nel,1) = vartmp(1:nel,1)
907 ipos(1:nel,2) = 1
908 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec ,yld ,hardp,hardr)
909 yld(1:nel) = yld(1:nel) * yscale
910 hardp(1:nel) = hardp(1:nel) * yscale
911 vartmp(1:nel,1) = ipos(1:nel,1)
912 ! Tabulation with Temperature
913 IF (numtabl == 2) THEN
914 ! Reference temperature factor
915 xvec(1:nel,2) = tini
916 ipos(1:nel,1) = vartmp(1:nel,2)
917 ipos(1:nel,2) = vartmp(1:nel,3)
918 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_tref,dydx)
919 vartmp(1:nel,2) = ipos(1:nel,1)
920 vartmp(1:nel,3) = ipos(1:nel,2)
921 ! Current temperature factor
922 xvec(1:nel,2) = temp(1:nel)
923 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_temp,dydx)
924 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
925 yld(1:nel) = yld(1:nel) * tfac(1:nel)
926 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
927 ELSE
928 tfac(1:nel) = one
929 ENDIF
930 ! Including kinematic hardening effect
931 yld(1:nel) = (one-fisokin)*yld(1:nel)+fisokin*yl0(1:nel)
932 ! Yield function value update
933 phi(1:nel) = sigvg(1:nel) - yld(1:nel)
934 ENDIF
935c
936 ! Kinematic hardening
937 IF (fisokin > zero) THEN
938 DO i=1,nel
939 dsxx = sigexx(i) - signxx(i)
940 dsyy = sigeyy(i) - signyy(i)
941 dsxy = sigexy(i) - signxy(i)
942 dexx = (dsxx - nu*dsyy)
943 deyy = (dsyy - nu*dsxx)
944 dexy = two*(one+nu)*dsxy
945 alpha = fisokin*hardp(i)/(young+hardp(i))*third
946 signxx(i) = signxx(i) + siga(i,1)
947 signyy(i) = signyy(i) + siga(i,2)
948 signxy(i) = signxy(i) + siga(i,3)
949 siga(i,1) = siga(i,1) + alpha*(four*dexx+two*deyy)
950 siga(i,2) = siga(i,2) + alpha*(four*deyy+two*dexx)
951 siga(i,3) = siga(i,3) + alpha*dexy
952 ENDDO
953 ENDIF
954c
955 ! Storing new values
956 DO i=1,nel
957 ! Computation of the plastic strain rate
958 IF (vflag == 1) THEN
959 dpdt = dpla(i)/max(em20,timestep)
960 uvar(i,1) = asrate * dpdt + (one - asrate) * uvar(i,1)
961 epsp(i) = uvar(i,1)
962 ENDIF
963 ! USR Outputs & coefficient for hourglass
964 IF (dpla(i) > zero) THEN
965 uvar(i,2) = phi(i) ! Yield function value
966 et(i) = hardp(i) / (hardp(i) + young)
967 ELSE
968 uvar(i,2) = zero
969 et(i) = one
970 ENDIF
971 seq(i) = sigvg(i) ! SIGEQ
972 ! Computation of the sound speed
973 soundsp(i) = sqrt(a11/rho(i))
974 ! Storing the yield stress
975 sigy(i) = yld(i)
976 ! Computation of the thickness variation
977 deelzz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young
978 ! Computation of the non-local thickness variation
979 IF ((inloc > 0).AND.(loff(i) == one)) THEN
980 ! Computation of the principal stresses
981 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
982 mohr_center = (signxx(i)+signyy(i))/two
983 sig1(i) = mohr_center + mohr_radius
984 sig2(i) = mohr_center - mohr_radius
985 IF (mohr_radius>em20) THEN
986 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
987 sin2theta(i) = signxy(i)/mohr_radius
988 ELSE
989 cos2theta(i) = one
990 sin2theta(i) = zero
991 ENDIF
992 ! Computation of scale factors for shear
993 fsh_theta(i,1:2) = zero
994 fps_theta(i,1:2) = zero
995 DO j = 1,nangle
996 DO k = 1,j
997 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
998 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
999 ENDDO
1000 ENDDO
1001 ! Computation of the equivalent stress of Vegter
1002 IF (sig1(i)<zero .OR. ((sig2(i)<zero) .AND. (sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2)))) THEN
1003 cos2theta(i) = -cos2theta(i)
1004 sin2theta(i) = -sin2theta(i)
1005 tmp = sig1(i)
1006 sig1(i) = -sig2(i)
1007 sig2(i) = -tmp
1008 sign_changed(i) = .true.
1009 ! Computation of scale factors
1010 fsh_theta(i,1:2) = zero
1011 fps_theta(i,1:2) = zero
1012 DO j = 1,nangle
1013 DO k = 1,j
1014 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1015 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1016 ENDDO
1017 ENDDO
1018 ELSE
1019 sign_changed(i) = .false.
1020 ENDIF
1021 ! between shear and uniaxial point
1022 IF (sig2(i)<zero) THEN
1023 ! Interpolation of factors
1024 fun_theta(i,1:2) = zero
1025 hish_theta(i,1:2) = zero
1026 DO j = 1,nangle
1027 DO k = 1,j
1028 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
1029 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1030 ENDDO
1031 ENDDO
1032 ! Filling the table
1033 a(1:2) = fsh_theta(i,1:2)
1034 b(1:2) = hish_theta(i,1:2)
1035 c(1:2) = fun_theta(i,1:2)
1036 ! Derivatives of PHI with respect to THETA
1037 da_dcos2(1:2) = zero
1038 db_dcos2(1:2) = zero
1039 dc_dcos2(1:2) = zero
1040 ! If anisotropic, compute derivatives with respect to COS2THET
1041 IF (nangle > 1) THEN
1042 ! Computation of their first derivative with respect to COS2THET
1043 DO j = 2, nangle
1044 DO k = 2, j
1045 da_dcos2(1:2) = da_dcos2(1:2) + q_fsh(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1046 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1047 dc_dcos2(1) = dc_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1048 ENDDO
1049 ENDDO
1050 ENDIF
1051 ! Between plane strain and uniaxial point
1052 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
1053 ! Interpolation of factors
1054 fun_theta(i,1:2) = zero
1055 hiun_theta(i,1:2) = zero
1056 DO j = 1,nangle
1057 DO k = 1,j
1058 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
1059 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1060 ENDDO
1061 ENDDO
1062 ! Filling the table
1063 a(1:2) = fun_theta(i,1:2)
1064 b(1:2) = hiun_theta(i,1:2)
1065 c(1:2) = fps_theta(i,1:2)
1066 ! Derivatives of PHI with respect to THETA
1067 da_dcos2(1:2) = zero
1068 db_dcos2(1:2) = zero
1069 dc_dcos2(1:2) = zero
1070 ! If anisotropic, compute derivatives with respect to COS2THET
1071 IF (nangle > 1) THEN
1072 ! Computation of their first derivative with respect to COS2THET
1073 DO j = 2, nangle
1074 DO k = 2,j
1075 da_dcos2(1) = da_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1076 db_dcos2(1:2) = db_dcos2(1:2) + q_hiun(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1077 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1078 ENDDO
1079 ENDDO
1080 ENDIF
1081 ! Between plane strain and equibiaxial point
1082 ELSE
1083 ! Interpolation of factors
1084 hips_theta(i,1:2) = zero
1085 DO j = 1,nangle
1086 DO k = 1,j
1087 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1088 ENDDO
1089 ENDDO
1090 ! Filling the table
1091 a(1:2) = fps_theta(i,1:2)
1092 b(1:2) = hips_theta(i,1:2)
1093 c(1:2) = fbi(1:2)
1094 ! Derivatives of PHI with respect to THETA
1095 da_dcos2(1:2) = zero
1096 db_dcos2(1:2) = zero
1097 dc_dcos2(1:2) = zero
1098 ! If anisotropic, compute derivatives with respect to COS2THET
1099 IF (nangle > 1) THEN
1100 ! Computation of their first derivative with respect to COS2THET
1101 DO j = 2, nangle
1102 DO k = 2,j
1103 da_dcos2(1:2) = da_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1104 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1105 ENDDO
1106 ENDDO
1107 ENDIF
1108 ENDIF
1109 ! Determine MU and SIGVG
1110 IF (sig1(i)<em20) THEN
1111 mu(i) = zero
1112 sigvg(i) = zero
1113 ELSE
1114 sig_ratio = sig2(i)/sig1(i)
1115 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
1116 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
1117 var_c = a(2) - sig_ratio*a(1)
1118 IF (abs(var_a)<em08) THEN
1119 mu(i) = -var_c/var_b
1120 ELSE
1121 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
1122 ENDIF
1123 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
1124 END IF
1125 ! Derivatives of PHI with respect to COS2THETA
1126 df1_dcos2 = da_dcos2(1) + two*mu(i)*(db_dcos2(1)-da_dcos2(1)) +
1127 . mu(i)*mu(i)*(da_dcos2(1)+dc_dcos2(1)-two*db_dcos2(1))
1128 df2_dcos2 = da_dcos2(2) + two*mu(i)*(db_dcos2(2)-da_dcos2(2)) +
1129 . mu(i)*mu(i)*(da_dcos2(2)+dc_dcos2(2)-two*db_dcos2(2))
1130 ! Computation of the derivatives
1131 f1 = mu(i)*mu(i)*(a(1)+c(1)-two*b(1))+two*mu(i)*(b(1)-a(1))+a(1)
1132 f2 = mu(i)*mu(i)*(a(2)+c(2)-two*b(2))+two*mu(i)*(b(2)-a(2))+a(2)
1133 ! Derivatives of Phi with respect to MU
1134 df1_dmu = two*(b(1)-a(1)) + two*mu(i)*(a(1)+c(1)-two*b(1))
1135 df2_dmu = two*(b(2)-a(2)) + two*mu(i)*(a(2)+c(2)-two*b(2))
1136 ! Derivatives of PHI with respect to SIG1, SIG2 and MU
1137 IF ((f1*df2_dmu - f2*df1_dmu)/=zero) THEN
1138 dphi_dsig1 = df2_dmu/(f1*df2_dmu - f2*df1_dmu)
1139 dphi_dsig2 = -df1_dmu/(f1*df2_dmu - f2*df1_dmu)
1140 IF (abs(sig1(i)-sig2(i))>tol) THEN
1141 dphi_dcos2 = (sigvg(i)/(sig1(i) - sig2(i)))*(df2_dcos2*df1_dmu -
1142 . df1_dcos2*df2_dmu)/(f1*df2_dmu - f2*df1_dmu)
1143 ELSE
1144 dphi_dcos2 = (two*((one-mu(i))*da_dcos2(2)+two*mu(i)*db_dcos2(2))*df1_dmu -
1145 . two*((one-mu(i))*da_dcos2(1)+two*mu(i)*db_dcos2(1))*df2_dmu)/
1146 . ((one-mu(i))*(a(1)-a(2)) + two*mu(i)*(b(1)-b(2)))
1147 ENDIF
1148 ELSE
1149 dphi_dsig1 = zero
1150 dphi_dsig2 = zero
1151 dphi_dcos2 = zero
1152 ENDIF
1153 normxx = half*(one + cos2theta(i))*dphi_dsig1 + half*(one-cos2theta(i))*dphi_dsig2 +
1154 . (sin2theta(i)**2)*dphi_dcos2
1155 normyy = half*(one - cos2theta(i))*dphi_dsig1 + half*(one+cos2theta(i))*dphi_dsig2 -
1156 . (sin2theta(i)**2)*dphi_dcos2
1157 normxy = sin2theta(i)*dphi_dsig1 - sin2theta(i)*dphi_dsig2 -
1158 . (two*sin2theta(i)*cos2theta(i))*dphi_dcos2
1159 IF (sign_changed(i)) THEN
1160 normxx = -normxx
1161 normyy = -normyy
1162 normxy = -normxy
1163 ENDIF
1164 sig_dfdsig = signxx(i) * normxx
1165 . + signyy(i) * normyy
1166 . + signxy(i) * normxy
1167 IF (sig_dfdsig > zero) THEN
1168 dezz(i) = -max(dplanl(i),zero)*(yld(i)/sig_dfdsig)*(normxx+normyy)
1169 ELSE
1170 dezz(i) = zero
1171 ENDIF
1172 ENDIF
1173 dezz(i) = deelzz(i) + dezz(i)
1174 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
1175 ENDDO
1176c
1177 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
alpha
#define alpha
Definition eval.h:35
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
interface_table_mod::table_vinterp
Definition table_mod.F:251
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
sigeps110c_nice
subroutine sigeps110c_nice(nel, ngl, nuparam, nuvar, npf, time, timestep, uparam, uvar, jthe, off, gs, rho, pla, dpla, epsp, soundsp, depsxx, depsyy, depsxy, depsyz, depszx, asrate, epspxx, epspyy, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thkly, thk, sigy, et, tempel, temp, seq, tf, numtabl, itable, table, nvartmp, vartmp, siga, inloc, dplanl, loff)
Definition sigeps110c_nice.F:45
table_mod::ttable
Definition table_mod.F:125