OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps110c_newton.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_newton (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_newton()

subroutine sigeps110c_newton ( 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_newton.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,dplanl,
70 . depsxx,depsyy,depsxy,depsyz,depszx,
71 . epspxx,epspyy,epspxy,epspyz,epspzx ,
72 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
73 . gs,thkly,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)
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
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 ! Standard inputs
222 dpla(i) = zero ! Initialization of the plastic strain increment
223 et(i) = one ! Initialization of hourglass coefficient
224 hardp(i) = zero ! Initialization of hourglass coefficient
225 sign_changed(i) = .false.
226 dezz(i) = zero
227 ENDDO
228c
229 ! /HEAT/MAT temperature
230 IF (jthe > 0) THEN
231 temp(1:nel) = tempel(1:nel)
232 ENDIF
233c
234 ! Computation of the strain-rate depending on the flags
235 IF ((vflag == 1) .OR. (vflag == 3)) THEN
236 DO i = 1, nel
237 IF (vflag == 1) THEN
238 epsp(i) = uvar(i,1)
239 ELSEIF (vflag == 3) THEN
240 dav = (epspxx(i)+epspyy(i))*third
241 deve1 = epspxx(i) - dav
242 deve2 = epspyy(i) - dav
243 deve3 = - dav
244 deve4 = half*epspxy(i)
245 epsp(i) = half*(deve1**2 + deve2**2 + deve3**2) + deve4**2
246 epsp(i) = sqrt(three*epsp(i))/three_half
247 epsp(i) = asrate*epsp(i) + (one - asrate)*uvar(i,1)
248 uvar(i,1) = epsp(i)
249 ENDIF
250 ENDDO
251 ENDIF
252c
253 ! Initial yield stress for kinematic hardening
254 IF (fisokin > zero) THEN
255 IF (numtabl > 0) THEN
256 xvec(1:nel,1) = zero
257 xvec(1:nel,2) = epsp(1:nel) * xscale
258 ipos0(1:nel,1) = 1
259 ipos0(1:nel,2) = 1
260 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos0,xvec ,yl0 ,hardp,hardr)
261 yl0(1:nel) = yl0(1:nel) * yscale
262 ! Tabulation with Temperature
263 IF (numtabl == 2) THEN
264 ! Reference temperature factor
265 xvec(1:nel,2) = tini
266 ipos0(1:nel,1) = 1
267 ipos0(1:nel,2) = 1
268 CALL table_vinterp(table(itable(2)),nel,nel,ipos0,xvec,yld_tref,dydx)
269 ! Current temperature factor
270 xvec(1:nel,2) = temp(1:nel)
271 CALL table_vinterp(table(itable(2)),nel,nel,ipos0,xvec,yld_temp,dydx)
272 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
273 yl0(1:nel) = yl0(1:nel) * tfac(1:nel)
274 ENDIF
275 ELSE
276 yl0(1:nel) = yld0
277 ENDIF
278 ELSE
279 yl0(1:nel) = zero
280 ENDIF
281c
282 ! Computation of the yield stress
283 IF (numtabl > 0) THEN
284 ! Tabulation with strain-rate
285 xvec(1:nel,1) = pla(1:nel)
286 xvec(1:nel,2) = epsp(1:nel) * xscale
287 ipos(1:nel,1) = vartmp(1:nel,1)
288 ipos(1:nel,2) = 1
289 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,yld,hardp,hardr)
290 yld(1:nel) = yld(1:nel) * yscale
291 hardp(1:nel) = hardp(1:nel) * yscale
292 vartmp(1:nel,1) = ipos(1:nel,1)
293 ! Tabulation with Temperature
294 IF (numtabl == 2) THEN
295 ! Reference temperature factor
296 xvec(1:nel,2) = tini
297 ipos(1:nel,1) = vartmp(1:nel,2)
298 ipos(1:nel,2) = vartmp(1:nel,3)
299 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_tref,dydx)
300 vartmp(1:nel,2) = ipos(1:nel,1)
301 vartmp(1:nel,3) = ipos(1:nel,2)
302 ! Current temperature factor
303 xvec(1:nel,2) = temp(1:nel)
304 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_temp,dydx)
305 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
306 yld(1:nel) = yld(1:nel) * tfac(1:nel)
307 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
308 ELSE
309 tfac(1:nel) = one
310 ENDIF
311 ! Including kinematic hardening effect
312 yld(1:nel) = (one-fisokin)*yld(1:nel)+fisokin*yl0(1:nel)
313 ELSE
314 DO i = 1,nel
315 ! a) - Hardening law
316 ! Continuous law
317 sighard(i) = yld0 + dsigm*(beta*(pla(i))+(one-exp(-omega*(pla(i))))**nexp)
318 ! b) - Strain-rate dependent hardening stress
319 sigrate(i) = sigst*(one + (kboltz*temp(i)/dg0)*log(one + (epsp(i)/deps0)))**mexp
320 ! c) - Computation of the yield function and value check
321 yld(i) = sighard(i) + sigrate(i)
322 ! d) - Including kinematic hardening
323 IF (fisokin > zero) THEN
324 yl0(i) = yl0(i) + sigrate(i)
325 ENDIF
326 yld(i) = (one-fisokin)*yld(i)+fisokin*yl0(i)
327 ! d) - Checking values
328 yld(i) = max(em10, yld(i))
329 ENDDO
330 ENDIF
331c
332 !========================================================================
333 ! - COMPUTATION OF TRIAL VALUES
334 !========================================================================
335 DO i=1,nel
336c
337 ! Computation of the trial stress tensor
338 IF (fisokin > zero) THEN
339 signxx(i) = sigoxx(i) - siga(i,1) + a11*depsxx(i) + a12*depsyy(i)
340 signyy(i) = sigoyy(i) - siga(i,2) + a12*depsxx(i) + a11*depsyy(i)
341 signxy(i) = sigoxy(i) - siga(i,3) + g*depsxy(i)
342 sigexx(i) = signxx(i)
343 sigeyy(i) = signyy(i)
344 sigexy(i) = signxy(i)
345 ELSE
346 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
347 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
348 signxy(i) = sigoxy(i) + depsxy(i)*g
349 ENDIF
350 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
351 signzx(i) = sigozx(i) + depszx(i)*gs(i)
352c
353 ! Computation of the equivalent stress
354 normsig = sqrt(signxx(i)*signxx(i)
355 . + signyy(i)*signyy(i)
356 . + two*signxy(i)*signxy(i))
357 IF (normsig < em20) THEN
358 sigvg(i) = zero
359 ELSE
360c
361 ! Computation of the principal stresses
362 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
363 mohr_center = (signxx(i)+signyy(i))/two
364 sig1(i) = mohr_center + mohr_radius
365 sig2(i) = mohr_center - mohr_radius
366 IF (mohr_radius>em20) THEN
367 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
368 sin2theta(i) = signxy(i)/mohr_radius
369 ELSE
370 cos2theta(i) = one
371 sin2theta(i) = zero
372 ENDIF
373c
374 ! Computation of scale factors for shear
375 fsh_theta(i,1:2) = zero
376 fps_theta(i,1:2) = zero
377 DO j = 1,nangle
378 DO k = 1,j
379 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
380 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
381 ENDDO
382 ENDDO
383c
384 ! Computation of the equivalent stress of Vegter
385 IF (sig1(i)<zero .OR. ((sig2(i)<zero) .AND. (sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2)))) THEN
386 cos2theta(i) = -cos2theta(i)
387 sin2theta(i) = -sin2theta(i)
388 tmp = sig1(i)
389 sig1(i) = -sig2(i)
390 sig2(i) = -tmp
391 sign_changed(i) = .true.
392 ! Computation of scale factors
393 fsh_theta(i,1:2) = zero
394 fps_theta(i,1:2) = zero
395 DO j = 1,nangle
396 DO k = 1,j
397 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
398 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
399 ENDDO
400 ENDDO
401 ELSE
402 sign_changed(i) = .false.
403 ENDIF
404 ! Between shear and uniaxial point
405 IF (sig2(i)<zero) THEN
406 ! Interpolation of factors
407 fun_theta(i,1:2) = zero
408 hish_theta(i,1:2) = zero
409 DO j = 1,nangle
410 DO k = 1,j
411 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
412 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
413 ENDDO
414 ENDDO
415 ! Filling the table
416 a(1:2) = fsh_theta(i,1:2)
417 b(1:2) = hish_theta(i,1:2)
418 c(1:2) = fun_theta(i,1:2)
419 ! Between plane strain and uniaxial point
420 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
421 ! Interpolation of factors
422 fun_theta(i,1:2) = zero
423 hiun_theta(i,1:2) = zero
424 DO j = 1,nangle
425 DO k = 1,j
426 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
427 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
428 ENDDO
429 ENDDO
430 ! Filling the table
431 a(1:2) = fun_theta(i,1:2)
432 b(1:2) = hiun_theta(i,1:2)
433 c(1:2) = fps_theta(i,1:2)
434 ! Between plane strain and equibiaxial point
435 ELSE
436 ! Interpolation of factors
437 hips_theta(i,1:2) = zero
438 DO j = 1,nangle
439 DO k = 1,j
440 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
441 ENDDO
442 ENDDO
443 ! Filling the table
444 a(1:2) = fps_theta(i,1:2)
445 b(1:2) = hips_theta(i,1:2)
446 c(1:2) = fbi(1:2)
447 ENDIF
448 ! Determine MU and SIGVG
449 IF (sig1(i)<em20) THEN
450 mu(i) = zero
451 sigvg(i) = zero
452 ELSE
453 sig_ratio = sig2(i)/sig1(i)
454 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
455 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
456 var_c = a(2) - sig_ratio*a(1)
457 IF (abs(var_a)<em08) THEN
458 mu(i) = -var_c/var_b
459 ELSE
460 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
461 ENDIF
462 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
463 END IF
464 ENDIF
465c
466 ENDDO
467c
468 !========================================================================
469 ! - COMPUTATION OF YIELD FONCTION
470 !========================================================================
471 phi(1:nel) = sigvg(1:nel) - yld(1:nel)
472c
473 ! Checking plastic behavior for all elements
474 nindx = 0
475 DO i=1,nel
476 IF (phi(i) > zero .AND. off(i) == one) THEN
477 nindx=nindx+1
478 index(nindx)=i
479 ENDIF
480 ENDDO
481c
482 !====================================================================
483 ! - PLASTIC CORRECTION WITH BACKWARD EULER METHOD (NEWTON RESOLUTION)
484 !====================================================================
485c
486 ! Number of Newton iterations
487 niter = 3
488c
489 ! Loop over yielding elements
490#include "vectorize.inc"
491 DO ii=1,nindx
492c
493 ! Number of the element with plastic behaviour
494 i = index(ii)
495c
496 ! Initialization of the iterative Newton procedure
497 dpxx(i) = zero
498 dpyy(i) = zero
499 dpzz(i) = zero
500 dpxy(i) = zero
501 dpyz(i) = zero
502 dpzx(i) = zero
503 ENDDO
504c
505 ! Loop over the iterations
506 DO iter = 1,niter
507#include "vectorize.inc"
508 ! Loop over yielding elements
509 DO ii=1,nindx
510 i = index(ii)
511c
512 ! Note : in this part, the purpose is to compute for each iteration
513 ! a plastic multiplier allowing to update internal variables to satisfy
514 ! the consistency condition using the backward Euler implicit method
515 ! with a Newton-Raphson iterative procedure
516 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
517 ! -> PHI : current value of yield function (known)
518 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
519 ! into account of internal variables kinetic :
520 ! plasticity, temperature and damage (to compute)
521c
522 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
523 !-------------------------------------------------------------
524c
525 ! Choice of loading conditions
526 ! Between shear and uniaxial point
527 IF (sig2(i)<zero) THEN
528 a(1:2) = fsh_theta(i,1:2)
529 b(1:2) = hish_theta(i,1:2)
530 c(1:2) = fun_theta(i,1:2)
531 ! Derivatives of PHI with respect to THETA
532 da_dcos2(1:2) = zero
533 db_dcos2(1:2) = zero
534 dc_dcos2(1:2) = zero
535 ! If anisotropic, compute derivatives with respect to COS2THET
536 IF (nangle > 1) THEN
537 ! Computation of their first derivative with respect to COS2THET
538 DO j = 2, nangle
539 DO k = 2, j
540 da_dcos2(1:2) = da_dcos2(1:2) + q_fsh(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
541 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
542 dc_dcos2(1) = dc_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
543 ENDDO
544 ENDDO
545 ENDIF
546 ! Between plane strain and uniaxial point
547 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
548 a(1:2) = fun_theta(i,1:2)
549 b(1:2) = hiun_theta(i,1:2)
550 c(1:2) = fps_theta(i,1:2)
551 ! Derivatives of PHI with respect to THETA
552 da_dcos2(1:2) = zero
553 db_dcos2(1:2) = zero
554 dc_dcos2(1:2) = zero
555 ! If anisotropic, compute derivatives with respect to COS2THET
556 IF (nangle > 1) THEN
557 ! Computation of their first derivative with respect to COS2THET
558 DO j = 2, nangle
559 DO k = 2,j
560 da_dcos2(1) = da_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
561 db_dcos2(1:2) = db_dcos2(1:2) + q_hiun(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
562 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
563 ENDDO
564 ENDDO
565 ENDIF
566 ! Between plane strain and equibiaxial point
567 ELSE
568 a(1:2) = fps_theta(i,1:2)
569 b(1:2) = hips_theta(i,1:2)
570 c(1:2) = fbi(1:2)
571 ! Derivatives of PHI with respect to THETA
572 da_dcos2(1:2) = zero
573 db_dcos2(1:2) = zero
574 dc_dcos2(1:2) = zero
575 ! If anisotropic, compute derivatives with respect to COS2THET
576 IF (nangle > 1) THEN
577 ! Computation of their first derivative with respect to COS2THET
578 DO j = 2, nangle
579 DO k = 2,j
580 da_dcos2(1:2) = da_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
581 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
582 ENDDO
583 ENDDO
584 ENDIF
585 ENDIF
586c
587 ! Derivatives of PHI with respect to COS2THETA
588 df1_dcos2 = da_dcos2(1) + two*mu(i)*(db_dcos2(1)-da_dcos2(1)) +
589 . mu(i)*mu(i)*(da_dcos2(1)+dc_dcos2(1)-two*db_dcos2(1))
590 df2_dcos2 = da_dcos2(2) + two*mu(i)*(db_dcos2(2)-da_dcos2(2)) +
591 . mu(i)*mu(i)*(da_dcos2(2)+dc_dcos2(2)-two*db_dcos2(2))
592c
593 ! Computation of the derivatives
594 f1 = mu(i)*mu(i)*(a(1)+c(1)-two*b(1))+two*mu(i)*(b(1)-a(1))+a(1)
595 f2 = mu(i)*mu(i)*(a(2)+c(2)-two*b(2))+two*mu(i)*(b(2)-a(2))+a(2)
596c
597 ! Derivatives of Phi with respect to MU
598 df1_dmu = two*(b(1)-a(1)) + two*mu(i)*(a(1)+c(1)-two*b(1))
599 df2_dmu = two*(b(2)-a(2)) + two*mu(i)*(a(2)+c(2)-two*b(2))
600c
601 ! Derivatives of PHI with respect to SIG1, SIG2 and MU
602 IF ((f1*df2_dmu - f2*df1_dmu)/=zero) THEN
603 dphi_dsig1 = df2_dmu/(f1*df2_dmu - f2*df1_dmu)
604 dphi_dsig2 = -df1_dmu/(f1*df2_dmu - f2*df1_dmu)
605 IF (abs(sig1(i)-sig2(i))>tol) THEN
606 dphi_dcos2 = (sigvg(i)/(sig1(i) - sig2(i)))*(df2_dcos2*df1_dmu -
607 . df1_dcos2*df2_dmu)/(f1*df2_dmu - f2*df1_dmu)
608 ELSE
609 dphi_dcos2 = (two*((one-mu(i))*da_dcos2(2)+two*mu(i)*db_dcos2(2))*df1_dmu -
610 . two*((one-mu(i))*da_dcos2(1)+two*mu(i)*db_dcos2(1))*df2_dmu)/
611 . ((one-mu(i))*(a(1)-a(2)) + two*mu(i)*(b(1)-b(2)))
612 ENDIF
613 ELSE
614 dphi_dsig1 = zero
615 dphi_dsig2 = zero
616 dphi_dcos2 = zero
617 ENDIF
618 normxx = half*(one + cos2theta(i))*dphi_dsig1 + half*(one-cos2theta(i))*dphi_dsig2 +
619 . (sin2theta(i)**2)*dphi_dcos2
620 normyy = half*(one - cos2theta(i))*dphi_dsig1 + half*(one+cos2theta(i))*dphi_dsig2 -
621 . (sin2theta(i)**2)*dphi_dcos2
622 normxy = sin2theta(i)*dphi_dsig1 - sin2theta(i)*dphi_dsig2 -
623 . (two*sin2theta(i)*cos2theta(i))*dphi_dcos2
624 IF (sign_changed(i)) THEN
625 normxx = -normxx
626 normyy = -normyy
627 normxy = -normxy
628 ENDIF
629c
630 ! 2 - Computation of DPHI_DLAMBDA
631 !---------------------------------------------------------
632c
633 ! a) Derivative with respect stress increments tensor DSIG
634 ! --------------------------------------------------------
635 dfdsig2 = normxx * (a11*normxx + a12*normyy)
636 . + normyy * (a11*normyy + a12*normxx)
637 . + normxy * normxy * g
638c
639 ! b) Derivatives with respect to plastic strain P
640 ! ------------------------------------------------
641c
642 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
643 ! ----------------------------------------------------------------------------
644 IF (numtabl == 0) THEN
645 ! Continuous hardening law
646 hardp(i) = dsigm*beta
647 IF (pla(i)>zero) THEN
648 hardp(i) = hardp(i) + dsigm*(nexp*(one-exp(-omega*(pla(i))))**(nexp-one))*
649 . (omega*exp(-omega*(pla(i))))
650 ENDIF
651 ENDIF
652 ! Accounting for kinematic hardening
653 dyld_dpla(i) = (one-fisokin)*hardp(i)
654c
655 ! ii) Derivative of dPLA with respect to DLAM
656 ! -------------------------------------------
657 sig_dfdsig = signxx(i) * normxx
658 . + signyy(i) * normyy
659 . + signxy(i) * normxy
660 dpla_dlam(i) = sig_dfdsig/yld(i)
661c
662 ! 3 - Computation of plastic multiplier and variables update
663 !----------------------------------------------------------
664c
665 ! Derivative of PHI with respect to DLAM
666 dphi_dlam(i) = - dfdsig2 - dyld_dpla(i)*dpla_dlam(i)
667 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
668c
669 ! Computation of the plastic multiplier
670 dlam = -phi(i)/dphi_dlam(i)
671c
672 ! Plastic strains tensor update
673 dpxx(i) = dlam * normxx
674 dpyy(i) = dlam * normyy
675 dpxy(i) = dlam * normxy
676c
677 ! Elasto-plastic stresses update
678 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
679 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
680 signxy(i) = signxy(i) - dpxy(i)*g
681c
682 ! Cumulated plastic strain and strain rate update
683 ddep = dlam*sig_dfdsig/yld(i)
684 dpla(i) = max(zero, dpla(i) + ddep)
685 pla(i) = pla(i) + ddep
686c
687 ! Computation of the equivalent stress
688 normsig = sqrt(signxx(i)*signxx(i)
689 . + signyy(i)*signyy(i)
690 . + two*signxy(i)*signxy(i))
691 IF (normsig < em20) THEN
692 sigvg(i) = zero
693 ELSE
694 ! Computation of the principal stresses
695 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
696 mohr_center = (signxx(i)+signyy(i))/two
697 sig1(i) = mohr_center + mohr_radius
698 sig2(i) = mohr_center - mohr_radius
699 IF (mohr_radius>em20) THEN
700 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
701 sin2theta(i) = signxy(i)/mohr_radius
702 ELSE
703 cos2theta(i) = one
704 sin2theta(i) = zero
705 ENDIF
706 ! Computation of scale factors
707 fsh_theta(i,1:2) = zero
708 fps_theta(i,1:2) = zero
709 DO j = 1,nangle
710 DO k = 1,j
711 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
712 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
713 ENDDO
714 ENDDO
715 ! Computation of the equivalent stress of Vegter
716 IF (sig1(i)<zero .OR. (sig2(i)<zero .AND. sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2))) THEN
717 cos2theta(i) = -cos2theta(i)
718 sin2theta(i) = -sin2theta(i)
719 tmp = sig1(i)
720 sig1(i) = -sig2(i)
721 sig2(i) = -tmp
722 sign_changed(i) = .true.
723 fsh_theta(i,1:2) = zero
724 fps_theta(i,1:2) = zero
725 DO j = 1,nangle
726 DO k = 1,j
727 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
728 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
729 ENDDO
730 ENDDO
731 ELSE
732 sign_changed(i) = .false.
733 ENDIF
734 ! Between shear and uniaxial point
735 IF (sig2(i)<zero) THEN
736 ! Interpolation of factors
737 fun_theta(i,1:2) = zero
738 hish_theta(i,1:2) = zero
739 DO j = 1,nangle
740 DO k = 1,j
741 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
742 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
743 ENDDO
744 ENDDO
745 ! Filling the tables
746 a(1:2) = fsh_theta(i,1:2)
747 b(1:2) = hish_theta(i,1:2)
748 c(1:2) = fun_theta(i,1:2)
749 ! Between plane strain and uniaxial point
750 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
751 ! Interpolation of factors
752 fun_theta(i,1:2) = zero
753 hiun_theta(i,1:2) = zero
754 DO j = 1,nangle
755 DO k = 1,j
756 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
757 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
758 ENDDO
759 ENDDO
760 ! Filling the tables
761 a(1:2) = fun_theta(i,1:2)
762 b(1:2) = hiun_theta(i,1:2)
763 c(1:2) = fps_theta(i,1:2)
764 ! Between plane strain and equibiaxial point
765 ELSE
766 ! Interpolation of factors
767 hips_theta(i,1:2) = zero
768 DO j = 1,nangle
769 DO k = 1,j
770 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
771 ENDDO
772 ENDDO
773 ! Filling the table
774 a(1:2) = fps_theta(i,1:2)
775 b(1:2) = hips_theta(i,1:2)
776 c(1:2) = fbi(1:2)
777 ENDIF
778 ! Determine MU and SIGVG
779 IF (sig1(i)<em20) THEN
780 mu(i) = zero
781 sigvg(i) = zero
782 ELSE
783 sig_ratio = sig2(i)/sig1(i)
784 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
785 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
786 var_c = a(2) - sig_ratio*a(1)
787 IF (abs(var_a)<em08) THEN
788 mu(i) = -var_c/var_b
789 ELSE
790 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
791 ENDIF
792 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
793 END IF
794 ENDIF
795c
796 IF (numtabl == 0) THEN
797 ! Continuous hardening law update
798 sighard(i) = yld0 + dsigm*(beta*(pla(i))+(one-exp(-omega*(pla(i))))**nexp)
799 ! Yield stress update
800 yld(i) = sighard(i) + sigrate(i)
801 yld(i) = (one-fisokin)*yld(i)+fisokin*yl0(i)
802 yld(i) = max(yld(i), em10)
803 ! Yield function value update
804 phi(i) = sigvg(i) - yld(i)
805 ENDIF
806c
807 ! Transverse strain update
808 IF (inloc == 0) THEN
809 dezz(i) = dezz(i) - (dpxx(i)+dpyy(i))
810 ENDIF
811c
812 ENDDO
813 ! End of the loop over yielding elements
814c
815 ! If tabulated yield function, update of the yield stress for all element
816 IF ((numtabl > 0).AND.(nindx > 0)) THEN
817 xvec(1:nel,1) = pla(1:nel)
818 xvec(1:nel,2) = epsp(1:nel) * xscale
819 ipos(1:nel,1) = vartmp(1:nel,1)
820 ipos(1:nel,2) = 1
821 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec ,yld ,hardp,hardr)
822 yld(1:nel) = yld(1:nel) * yscale
823 hardp(1:nel) = hardp(1:nel) * yscale
824 vartmp(1:nel,1) = ipos(1:nel,1)
825 ! Tabulation with Temperature
826 IF (numtabl == 2) THEN
827 ! Reference temperature factor
828 xvec(1:nel,2) = tini
829 ipos(1:nel,1) = vartmp(1:nel,2)
830 ipos(1:nel,2) = vartmp(1:nel,3)
831 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_tref,dydx)
832 vartmp(1:nel,2) = ipos(1:nel,1)
833 vartmp(1:nel,3) = ipos(1:nel,2)
834 ! Current temperature factor
835 xvec(1:nel,2) = temp(1:nel)
836 CALL table_vinterp(table(itable(2)),nel,nel,ipos,xvec,yld_temp,dydx)
837 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
838 yld(1:nel) = yld(1:nel) * tfac(1:nel)
839 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
840 ELSE
841 tfac(1:nel) = one
842 ENDIF
843 ! Including kinematic hardening effect
844 yld(1:nel) = (one-fisokin)*yld(1:nel)+fisokin*yl0(1:nel)
845 ! Yield function value update
846 phi(1:nel) = sigvg(1:nel) - yld(1:nel)
847 ENDIF
848c
849 ENDDO
850 ! End of the loop over the iterations
851 !===================================================================
852 ! - END OF PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
853 !===================================================================
854c
855 ! Kinematic hardening
856 IF (fisokin > zero) THEN
857 DO i=1,nel
858 dsxx = sigexx(i) - signxx(i)
859 dsyy = sigeyy(i) - signyy(i)
860 dsxy = sigexy(i) - signxy(i)
861 dexx = (dsxx - nu*dsyy)
862 deyy = (dsyy - nu*dsxx)
863 dexy = two*(one+nu)*dsxy
864 alpha = fisokin*hardp(i)/(young+hardp(i))*third
865 signxx(i) = signxx(i) + siga(i,1)
866 signyy(i) = signyy(i) + siga(i,2)
867 signxy(i) = signxy(i) + siga(i,3)
868 siga(i,1) = siga(i,1) + alpha*(four*dexx+two*deyy)
869 siga(i,2) = siga(i,2) + alpha*(four*deyy+two*dexx)
870 siga(i,3) = siga(i,3) + alpha*dexy
871 ENDDO
872 ENDIF
873c
874 ! Storing new values
875 DO i=1,nel
876 ! Computation of the plastic strain rate
877 IF (vflag == 1) THEN
878 dpdt = dpla(i)/max(em20,timestep)
879 uvar(i,1) = asrate * dpdt + (one - asrate) * uvar(i,1)
880 epsp(i) = uvar(i,1)
881 ENDIF
882 ! USR Outputs
883 seq(i) = sigvg(i) ! SIGEQ
884 ! Coefficient for hourglass
885 IF (dpla(i) > zero) THEN
886 et(i) = hardp(i) / (hardp(i) + young)
887 ELSE
888 et(i) = one
889 ENDIF
890 ! Computation of the sound speed
891 soundsp(i) = sqrt(a11/rho(i))
892 ! Storing the yield stress
893 sigy(i) = yld(i)
894 ! Computation of the thickness variation
895 deelzz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young
896 ! Computation of the non-local thickness variation
897 IF ((inloc > 0).AND.(loff(i) == one)) THEN
898 ! Computation of the old principal stresses
899 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
900 mohr_center = (signxx(i)+signyy(i))/two
901 sig1(i) = mohr_center + mohr_radius
902 sig2(i) = mohr_center - mohr_radius
903 IF (mohr_radius>em20) THEN
904 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
905 sin2theta(i) = signxy(i)/mohr_radius
906 ELSE
907 cos2theta(i) = one
908 sin2theta(i) = zero
909 ENDIF
910 ! Computation of scale factors for shear
911 fsh_theta(i,1:2) = zero
912 fps_theta(i,1:2) = zero
913 DO j = 1,nangle
914 DO k = 1,j
915 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
916 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
917 ENDDO
918 ENDDO
919 ! Computation of the equivalent stress of Vegter
920 IF (sig1(i)<zero .OR. ((sig2(i)<zero) .AND. (sig2(i)*fsh_theta(i,1)<sig1(i)*fsh_theta(i,2)))) THEN
921 cos2theta(i) = -cos2theta(i)
922 sin2theta(i) = -sin2theta(i)
923 tmp = sig1(i)
924 sig1(i) = -sig2(i)
925 sig2(i) = -tmp
926 sign_changed(i) = .true.
927 ! Computation of scale factors
928 fsh_theta(i,1:2) = zero
929 fps_theta(i,1:2) = zero
930 DO j = 1,nangle
931 DO k = 1,j
932 fsh_theta(i,1:2) = fsh_theta(i,1:2) + q_fsh(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
933 fps_theta(i,1:2) = fps_theta(i,1:2) + q_fps(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
934 ENDDO
935 ENDDO
936 ELSE
937 sign_changed(i) = .false.
938 ENDIF
939 ! Between shear and uniaxial point
940 IF (sig2(i)<zero) THEN
941 ! Interpolation of factors
942 fun_theta(i,1:2) = zero
943 hish_theta(i,1:2) = zero
944 DO j = 1,nangle
945 DO k = 1,j
946 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
947 hish_theta(i,1:2) = hish_theta(i,1:2) + q_hish(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
948 ENDDO
949 ENDDO
950 ! Filling the table
951 a(1:2) = fsh_theta(i,1:2)
952 b(1:2) = hish_theta(i,1:2)
953 c(1:2) = fun_theta(i,1:2)
954 ! Derivatives of PHI with respect to THETA
955 da_dcos2(1:2) = zero
956 db_dcos2(1:2) = zero
957 dc_dcos2(1:2) = zero
958 ! If anisotropic, compute derivatives with respect to COS2THET
959 IF (nangle > 1) THEN
960 ! Computation of their first derivative with respect to COS2THET
961 DO j = 2, nangle
962 DO k = 2, j
963 da_dcos2(1:2) = da_dcos2(1:2) + q_fsh(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
964 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
965 dc_dcos2(1) = dc_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
966 ENDDO
967 ENDDO
968 ENDIF
969 ! Between plane strain and uniaxial point
970 ELSEIF (sig2(i)/sig1(i)<fps_theta(i,2)/fps_theta(i,1)) THEN
971 ! Interpolation of factors
972 fun_theta(i,1:2) = zero
973 hiun_theta(i,1:2) = zero
974 DO j = 1,nangle
975 DO k = 1,j
976 fun_theta(i,1) = fun_theta(i,1) + q_fun(j)*cos2(k,j)*(cos2theta(i))**(k-1)
977 hiun_theta(i,1:2) = hiun_theta(i,1:2) + q_hiun(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
978 ENDDO
979 ENDDO
980 ! Filling the table
981 a(1:2) = fun_theta(i,1:2)
982 b(1:2) = hiun_theta(i,1:2)
983 c(1:2) = fps_theta(i,1:2)
984 ! Derivatives of PHI with respect to THETA
985 da_dcos2(1:2) = zero
986 db_dcos2(1:2) = zero
987 dc_dcos2(1:2) = zero
988 ! If anisotropic, compute derivatives with respect to COS2THET
989 IF (nangle > 1) THEN
990 ! Computation of their first derivative with respect to COS2THET
991 DO j = 2, nangle
992 DO k = 2,j
993 da_dcos2(1) = da_dcos2(1) + q_fun(j)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
994 db_dcos2(1:2) = db_dcos2(1:2) + q_hiun(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
995 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
996 ENDDO
997 ENDDO
998 ENDIF
999 ! Between plane strain and equibiaxial point
1000 ELSE
1001 ! Interpolation of factors
1002 hips_theta(i,1:2) = zero
1003 DO j = 1,nangle
1004 DO k = 1,j
1005 hips_theta(i,1:2) = hips_theta(i,1:2) + q_hips(j,1:2)*cos2(k,j)*(cos2theta(i))**(k-1)
1006 ENDDO
1007 ENDDO
1008 ! Filling the table
1009 a(1:2) = fps_theta(i,1:2)
1010 b(1:2) = hips_theta(i,1:2)
1011 c(1:2) = fbi(1:2)
1012 ! Derivatives of PHI with respect to THETA
1013 da_dcos2(1:2) = zero
1014 db_dcos2(1:2) = zero
1015 dc_dcos2(1:2) = zero
1016 ! If anisotropic, compute derivatives with respect to COS2THET
1017 IF (nangle > 1) THEN
1018 ! Computation of their first derivative with respect to COS2THET
1019 DO j = 2, nangle
1020 DO k = 2,j
1021 da_dcos2(1:2) = da_dcos2(1:2) + q_fps(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1022 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(j,1:2)*(k-1)*cos2(k,j)*(cos2theta(i))**(k-2)
1023 ENDDO
1024 ENDDO
1025 ENDIF
1026 ENDIF
1027 ! Determine MU and SIGVG
1028 IF (sig1(i)<em20) THEN
1029 mu(i) = zero
1030 sigvg(i) = zero
1031 ELSE
1032 sig_ratio = sig2(i)/sig1(i)
1033 var_a = (c(2)+a(2)-two*b(2)) - sig_ratio*(c(1)+a(1)-two*b(1))
1034 var_b = two*((b(2)-a(2)) - sig_ratio*(b(1)-a(1)))
1035 var_c = a(2) - sig_ratio*a(1)
1036 IF (abs(var_a)<em08) THEN
1037 mu(i) = -var_c/var_b
1038 ELSE
1039 mu(i) = (-var_b+sqrt(var_b*var_b - four*var_a*var_c))/(two*var_a)
1040 ENDIF
1041 sigvg(i) = sig1(i)/(a(1)+two*mu(i)*(b(1)-a(1))+mu(i)*mu(i)*(c(1)+a(1)-two*b(1)))
1042 END IF
1043 ! Derivatives of PHI with respect to COS2THETA
1044 df1_dcos2 = da_dcos2(1) + two*mu(i)*(db_dcos2(1)-da_dcos2(1)) +
1045 . mu(i)*mu(i)*(da_dcos2(1)+dc_dcos2(1)-two*db_dcos2(1))
1046 df2_dcos2 = da_dcos2(2) + two*mu(i)*(db_dcos2(2)-da_dcos2(2)) +
1047 . mu(i)*mu(i)*(da_dcos2(2)+dc_dcos2(2)-two*db_dcos2(2))
1048 ! Computation of the derivatives
1049 f1 = mu(i)*mu(i)*(a(1)+c(1)-two*b(1))+two*mu(i)*(b(1)-a(1))+a(1)
1050 f2 = mu(i)*mu(i)*(a(2)+c(2)-two*b(2))+two*mu(i)*(b(2)-a(2))+a(2)
1051 ! Derivatives of Phi with respect to MU
1052 df1_dmu = two*(b(1)-a(1)) + two*mu(i)*(a(1)+c(1)-two*b(1))
1053 df2_dmu = two*(b(2)-a(2)) + two*mu(i)*(a(2)+c(2)-two*b(2))
1054 ! Derivatives of PHI with respect to SIG1, SIG2 and MU
1055 IF ((f1*df2_dmu - f2*df1_dmu)/=zero) THEN
1056 dphi_dsig1 = df2_dmu/(f1*df2_dmu - f2*df1_dmu)
1057 dphi_dsig2 = -df1_dmu/(f1*df2_dmu - f2*df1_dmu)
1058 IF (abs(sig1(i)-sig2(i))>tol) THEN
1059 dphi_dcos2 = (sigvg(i)/(sig1(i) - sig2(i)))*(df2_dcos2*df1_dmu -
1060 . df1_dcos2*df2_dmu)/(f1*df2_dmu - f2*df1_dmu)
1061 ELSE
1062 dphi_dcos2 = (two*((one-mu(i))*da_dcos2(2)+two*mu(i)*db_dcos2(2))*df1_dmu -
1063 . two*((one-mu(i))*da_dcos2(1)+two*mu(i)*db_dcos2(1))*df2_dmu)/
1064 . ((one-mu(i))*(a(1)-a(2)) + two*mu(i)*(b(1)-b(2)))
1065 ENDIF
1066 ELSE
1067 dphi_dsig1 = zero
1068 dphi_dsig2 = zero
1069 dphi_dcos2 = zero
1070 ENDIF
1071 normxx = half*(one + cos2theta(i))*dphi_dsig1 + half*(one-cos2theta(i))*dphi_dsig2 +
1072 . (sin2theta(i)**2)*dphi_dcos2
1073 normyy = half*(one - cos2theta(i))*dphi_dsig1 + half*(one+cos2theta(i))*dphi_dsig2 -
1074 . (sin2theta(i)**2)*dphi_dcos2
1075 normxy = sin2theta(i)*dphi_dsig1 - sin2theta(i)*dphi_dsig2 -
1076 . (two*sin2theta(i)*cos2theta(i))*dphi_dcos2
1077 IF (sign_changed(i)) THEN
1078 normxx = -normxx
1079 normyy = -normyy
1080 normxy = -normxy
1081 ENDIF
1082 sig_dfdsig = signxx(i) * normxx
1083 . + signyy(i) * normyy
1084 . + signxy(i) * normxy
1085 IF (sig_dfdsig > zero) THEN
1086 dezz(i) = -max(dplanl(i),zero)*(yld(i)/sig_dfdsig)*(normxx+normyy)
1087 ELSE
1088 dezz(i) = zero
1089 ENDIF
1090 ENDIF
1091 dezz(i) = deelzz(i) + dezz(i)
1092 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
1093 ENDDO
1094c
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
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112