OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps110c_newton.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_newton ../engine/source/materials/mat/mat110/sigeps110c_newton.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_newton(
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,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
1095 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_newton
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)
Definition sigeps110c_newton.F:45
table_mod::ttable
Definition table_mod.F:125