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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps110c_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)

Function/Subroutine Documentation

◆ sigeps110c_lite_nice()

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

Definition at line 34 of file sigeps110c_lite_nice.F.

45 !=======================================================================
46 ! Modules
47 !=======================================================================
48 USE table_mod
49 USE interface_table_mod
50 !=======================================================================
51 ! Implicit types
52 !=======================================================================
53#include "implicit_f.inc"
54 !=======================================================================
55 ! Common
56 !=======================================================================
57#include "com04_c.inc"
58 !=======================================================================
59 ! Dummy arguments
60 !=======================================================================
61 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP,NPF(*),INLOC
62 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
63 my_real
64 . time,timestep,asrate,tf(*)
65 INTEGER :: VARTMP(NEL,NVARTMP)
66 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
67 . uparam
68 my_real,DIMENSION(NEL), INTENT(IN) ::
69 . rho,tempel,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
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
interface_table_mod::table_vinterp
Definition table_mod.F:251
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112