OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps110c_nice.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ sigeps110c_nice()

subroutine sigeps110c_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
integer, dimension(*) npf,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) gs,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(out) epsp,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
asrate,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(in) thkly,
intent(inout) thk,
intent(out) sigy,
intent(out) et,
intent(in) tempel,
intent(inout) temp,
intent(inout) seq,
tf,
integer numtabl,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table,
integer nvartmp,
integer, dimension(nel,nvartmp) vartmp,
intent(inout) siga,
integer inloc,
intent(in) dplanl,
intent(in) loff )

Definition at line 34 of file sigeps110c_nice.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
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
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112