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