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