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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ sigeps110c_lite_newton()

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

Definition at line 34 of file sigeps110c_lite_newton.F.

45 !=======================================================================
46 ! Modules
47 !=======================================================================
48 USE table_mod
49 USE interface_table_mod
50 !=======================================================================
51 ! Implicit types
52 !=======================================================================
53#include "implicit_f.inc"
54 !=======================================================================
55 ! Common
56 !=======================================================================
57#include "com04_c.inc"
58 !=======================================================================
59 ! Dummy arguments
60 !=======================================================================
61 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP,NPF(*),INLOC
62 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
63 my_real
64 . time,timestep,asrate,tf(*)
65 INTEGER :: VARTMP(NEL,NVARTMP)
66 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
67 . uparam
68 my_real,DIMENSION(NEL), INTENT(IN) ::
69 . rho,tempel,dplanl,
70 . depsxx,depsyy,depsxy,depsyz,depszx,
71 . epspxx,epspyy,epspxy,epspyz,epspzx ,
72 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
73 . gs,thkly,loff
74c
75 my_real ,DIMENSION(NEL), INTENT(OUT) ::
76 . soundsp,sigy,et,epsp,
77 . signxx,signyy,signxy,signyz,signzx
78c
79 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
80 . pla,dpla,off,thk,temp,seq
81 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
82 . uvar
83 my_real ,DIMENSION(NEL,3) ,INTENT(INOUT) ::
84 . siga
85c
86 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
87 !=======================================================================
88 ! Local Variables
89 !=======================================================================
90 INTEGER I,J,II,K,ITER,NITER,NINDX,INDEX(NEL),VFLAG,IPOS(NEL,2),NANGLE,
91 . IPOS0(NEL,2),ISMOOTH
92c
93 my_real
94 . young,g,g2,nu,nnu,a11,a12,yld0,dsigm,beta,omega,nexp,eps0,sigst,
95 . dg0,deps0,mexp,fbi(2),kboltz,fisokin,tini,xscale,yscale
96 my_real
97 . mohr_radius,mohr_center,normsig,tmp,sig_ratio,var_a,var_b,var_c,
98 . a(2),b(2),c(2),h,dpdt,dlam,ddep,dav,deve1,deve2,deve3,deve4,
99 . xvec(nel,4)
100 my_real
101 . fun_theta(nel,2),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
my_real
#define my_real
Definition cppsort.cpp:32
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