OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112_xia_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 mat112_xia_nice (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, epszz, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, nvartmp, numtabl, vartmp, itable, table)

Function/Subroutine Documentation

◆ mat112_xia_nice()

subroutine mat112_xia_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) rho0,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) epszz,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
integer nvartmp,
integer numtabl,
integer, dimension(nel,nvartmp) vartmp,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table )

Definition at line 33 of file mat112_xia_nice.F.

42 !=======================================================================
43 ! Modules
44 !=======================================================================
45 USE table_mod
46 USE interface_table_mod
47 !=======================================================================
48 ! Implicit types
49 !=======================================================================
50#include "implicit_f.inc"
51 !=======================================================================
52 ! Common
53 !=======================================================================
54#include "com04_c.inc"
55 !=======================================================================
56 ! Dummy arguments
57 !=======================================================================
58 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP
59 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
60 my_real
61 . time,timestep
62 INTEGER :: VARTMP(NEL,NVARTMP)
63 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
64 . uparam
65 my_real,DIMENSION(NEL), INTENT(IN) ::
66 . rho,rho0,epszz,
67 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL), INTENT(OUT) ::
71 . soundsp,sigy,et,
72 . signxx,signyy,signzz,signxy,signyz,signzx
73c
74 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
75 . dpla,off
76 my_real ,DIMENSION(NEL,4),INTENT(INOUT) ::
77 . pla,epsd
78 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
79 . uvar
80c
81 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
82 !=======================================================================
83 ! Local Variables
84 !=======================================================================
85 INTEGER I,K,II,NINDX,INDEX(NEL),NINDX2,INDEX2(NEL),NINDX3,INDEX3(NEL),
86 . ITAB,ISMOOTH,IPOS(NEL,2)
87c
88 my_real
89 . young1,young2,young3,nu12,nu21,
90 . a11,a12,a21,a22,
91 . g12,g23,g31,deuxk,e3c,cc,c1,
92 . nu1p,nu2p,nu4p,nu5p,s01,a01,b01,c01,
93 . s02,a02,b02,c02,s03,a03,b03,c03,
94 . s04,a04,b04,c04,s05,a05,b05,c05,
95 . asig,bsig,csig,tau0,atau,btau,n3,n6,
96 . bulk,nu13,nu31,nu23,nu32
97 my_real
98 . normsig,h,dpdt,dlam,ddep,depxx,depyy
99 my_real
100 . normxx,normyy,normxy,normpxx,normpyy,normpxy,
101 . dphi_dlam,dphi,dfdsig2,dphi_dpla,dpxx,dpyy,dpxy,
102 . dphi_dsigy1,dphi_dsigy2,dphi_dsigy3,dphi_dsigy4,
103 . dphi_dsigy5
104 my_real
105 . normzz,dtheta_dlam,dtheta,dtheta_dpla,
106 . dtheta_dsigyo
107 my_real
108 . normyz,normzx,dpsi_dlam,dpsi,dpsi_dpla,
109 . dpyz,dpzx,dpsi_dsigys
110 my_real
111 . xscale1,yscale1,xscale2,yscale2,xscale3,yscale3,xscale4,yscale4,
112 . xscale5,yscale5,xscalec,yscalec,xscales,yscales,xvec(nel,2),asrate
113c
114 my_real, DIMENSION(2) ::
115 . n1,n2,n4,n5
116c
117 my_real, DIMENSION(NEL) ::
118 . khi1,khi2,khi3,khi4,khi5,khi6,sigy1,sigy2,sigy3,sigy4,sigy5,dpla2,dpla3,
119 . dpla4,sigys,sigyo,dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,hardp,phi,phi0,psi,psi0,
120 . theta,theta0,dsigzz,eezz,gam1,gam2,gam3,gam4,gam5,gam6,dsigy1_dp,dsigy2_dp,
121 . dsigy3_dp,dsigy4_dp,dsigy5_dp,dsigyo_dp,dsigys_dp,hardr
122c
123 !=======================================================================
124 ! - INITIALISATION OF COMPUTATION ON TIME STEP
125 !=======================================================================
126 ! Recovering model parameters
127 ! Elastic parameters
128 young1 = uparam(1) ! Young modulus in direction 1 (MD)
129 young2 = uparam(2) ! Young modulus in direction 2 (CD)
130 young3 = uparam(3) ! Young modulus in direction 3 (Out of plane)
131 nu12 = uparam(4) ! Poisson's ratio in 12
132 nu21 = uparam(5) ! Poisson's ratio in 21
133 a11 = uparam(6) ! Component 11 of orthotropic shell elasticity matrix
134 a12 = uparam(7) ! Component 12 of orthotropic shell elasticity matrix
135 a21 = uparam(8) ! Component 21 of orthotropic shell elasticity matrix
136 a22 = uparam(9) ! Component 22 of orthotropic shell elasticity matrix
137 g12 = uparam(10) ! Shear modulus in 12
138 g23 = uparam(11) ! Shear modulus in 23
139 g31 = uparam(12) ! Shear modulus in 31
140 itab = int(uparam(14)) ! Tabulated yield stress flag
141 deuxk = uparam(15) ! Yield criterion exponent
142 e3c = uparam(16) ! Compressive elasticity parameter
143 cc = uparam(17) ! Exponent compressive elasticity parameter
144 nu1p = uparam(18) ! Tensile plastic Poisson ratio in direction 1 (MD)
145 nu2p = uparam(19) ! Tensile plastic Poisson ratio in direction 2 (CD)
146 nu4p = uparam(20) ! Compressive plastic Poisson ratio in direction 1 (MD)
147 nu5p = uparam(21) ! Compressive plastic Poisson ratio in direction 2 (CD)
148 IF (itab == 0) THEN
149 s01 = uparam(22) ! 1st Tensile plasticity parameter direction 1 (MD)
150 a01 = uparam(23) ! 2nd Tensile plasticity parameter direction 1 (MD)
151 b01 = uparam(24) ! 3rd Tensile plasticity parameter direction 1 (MD)
152 c01 = uparam(25) ! 4th Tensile plasticity parameter direction 1 (MD)
153 s02 = uparam(26) ! 1st Tensile plasticity parameter direction 2 (CD)
154 a02 = uparam(27) ! 2nd Tensile plasticity parameter direction 2 (CD)
155 b02 = uparam(28) ! 3rd Tensile plasticity parameter direction 2 (CD)
156 c02 = uparam(29) ! 4th Tensile plasticity parameter direction 2 (CD)
157 s03 = uparam(30) ! 1st Positive shear plasticity parameter
158 a03 = uparam(31) ! 2nd Positive shear plasticity parameter
159 b03 = uparam(32) ! 3rd Positive shear plasticity parameter
160 c03 = uparam(33) ! 4th Positive shear plasticity parameter
161 s04 = uparam(34) ! 1st Compressive plasticity parameter direction 1 (MD)
162 a04 = uparam(35) ! 2nd Compressive plasticity parameter direction 1 (MD)
163 b04 = uparam(36) ! 3rd Compressive plasticity parameter direction 1 (MD)
164 c04 = uparam(37) ! 4th Compressive plasticity parameter direction 1 (MD)
165 s05 = uparam(38) ! 1st Compressive plasticity parameter direction 2 (CD)
166 a05 = uparam(39) ! 2nd Compressive plasticity parameter direction 2 (CD)
167 b05 = uparam(40) ! 3rd Compressive plasticity parameter direction 2 (CD)
168 c05 = uparam(41) ! 4th Compressive plasticity parameter direction 2 (CD)
169 asig = uparam(42) ! 1st Out-of-plane plasticity parameter
170 bsig = uparam(43) ! 2nd Out-of-plane plasticity parameter
171 csig = uparam(44) ! 3rd Out-of-plane plasticity parameter
172 tau0 = uparam(45) ! 1st Transverse shear plasticity parameter
173 atau = uparam(46) ! 2nd Transverse shear plasticity parameter
174 btau = uparam(47) ! 3rd Transverse shear plasticity parameter
175 ELSE
176 xscale1 = uparam(22)
177 yscale1 = uparam(23)
178 xscale2 = uparam(24)
179 yscale2 = uparam(25)
180 xscale3 = uparam(26)
181 yscale3 = uparam(27)
182 xscale4 = uparam(28)
183 yscale4 = uparam(29)
184 xscale5 = uparam(30)
185 yscale5 = uparam(31)
186 xscalec = uparam(32)
187 yscalec = uparam(33)
188 xscales = uparam(34)
189 yscales = uparam(35)
190 asrate = uparam(36)
191 asrate = (two*pi*asrate*timestep)/(two*pi*asrate*timestep + one)
192 ismooth = int(uparam(37))
193 ENDIF
194c
195 ! Yield planes normal
196 n1(1) = one/(sqrt(one+nu1p**2))
197 n1(2) = -nu1p/(sqrt(one+nu1p**2))
198 n2(1) = -nu2p/(sqrt(one+nu2p**2))
199 n2(2) = one/(sqrt(one+nu2p**2))
200 n3 = one
201 n4(1) = -one/(sqrt(one+nu4p**2))
202 n4(2) = nu4p/(sqrt(one+nu4p**2))
203 n5(1) = nu5p/(sqrt(one+nu5p**2))
204 n5(2) = -one/(sqrt(one+nu5p**2))
205 n6 = -one
206c
207 ! Recovering internal variables
208 DO i=1,nel
209 ! OFF parameter for element deletion
210 IF (off(i) < em01) off(i) = zero
211 IF (off(i) < one) off(i) = off(i)*four_over_5
212 ! User inputs
213 phi0(i) = uvar(i,1) ! Previous in-plane yield function value
214 theta0(i) = uvar(i,2) ! Previous out-of-plane yield function value
215 psi0(i) = uvar(i,3) ! Previous transverse shear yield function value
216 eezz(i) = epszz(i) + pla(i,3) ! Out-of-plane elastic strain
217 ! Standard inputs
218 dpla(i) = zero ! Initialization of the "global" plastic strain increment
219 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
220 dpla3(i) = zero ! Initialization of the out-of-plane plastic strain increment
221 dpla4(i) = zero ! Initialization of the transverse shear plastic strain increment
222 et(i) = one ! Initialization of hourglass coefficient
223 hardp(i) = zero ! Initialization of hourglass coefficient
224 ENDDO
225c
226 ! Computation of the initial yield stress
227 IF (itab > 0) THEN
228 ! In-plane tabulation with strain-rate
229 xvec(1:nel,1) = pla(1:nel,2)
230 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
231 ! -> Tensile yield stress in direction 1 (MD)
232 ipos(1:nel,1) = vartmp(1:nel,1)
233 ipos(1:nel,2) = 1
234 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
235 sigy1(1:nel) = sigy1(1:nel) * yscale1
236 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
237 vartmp(1:nel,1) = ipos(1:nel,1)
238 ! -> Tensile yield stress in direction 2 (CD)
239 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
240 ipos(1:nel,1) = vartmp(1:nel,2)
241 ipos(1:nel,2) = 1
242 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
243 sigy2(1:nel) = sigy2(1:nel) * yscale2
244 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
245 vartmp(1:nel,2) = ipos(1:nel,1)
246 ! -> Positive shear yield stress
247 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
248 ipos(1:nel,1) = vartmp(1:nel,3)
249 ipos(1:nel,2) = 1
250 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
251 sigy3(1:nel) = sigy3(1:nel) * yscale3
252 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
253 vartmp(1:nel,3) = ipos(1:nel,1)
254 ! -> Compressive yield stress in direction 1 (MD)
255 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
256 ipos(1:nel,1) = vartmp(1:nel,4)
257 ipos(1:nel,2) = 1
258 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
259 sigy4(1:nel) = sigy4(1:nel) * yscale4
260 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
261 vartmp(1:nel,4) = ipos(1:nel,1)
262 ! -> Compressive yield stress in direction 2 (CD)
263 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
264 ipos(1:nel,1) = vartmp(1:nel,5)
265 ipos(1:nel,2) = 1
266 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
267 sigy5(1:nel) = sigy5(1:nel) * yscale5
268 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
269 vartmp(1:nel,5) = ipos(1:nel,1)
270 ! -> Out-of-plane tabulation with strain-rate
271 xvec(1:nel,1) = pla(1:nel,3)
272 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
273 ! -> Transverse shear yield stress
274 ipos(1:nel,1) = vartmp(1:nel,6)
275 ipos(1:nel,2) = 1
276 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
277 sigyo(1:nel) = sigyo(1:nel) * yscalec
278 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
279 vartmp(1:nel,6) = ipos(1:nel,1)
280 ! -> Transverse shear tabulation with strain-rate
281 xvec(1:nel,1) = pla(1:nel,4)
282 xvec(1:nel,2) = epsd(1:nel,4) * xscales
283 ! -> Transverse shear yield stress
284 ipos(1:nel,1) = vartmp(1:nel,7)
285 ipos(1:nel,2) = 1
286 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
287 sigys(1:nel) = sigys(1:nel) * yscales
288 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
289 vartmp(1:nel,7) = ipos(1:nel,1)
290 ELSE
291 DO i = 1,nel
292 ! Tensile yield stress in direction 1 (MD)
293 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
294 ! Tensile yield stress in direction 2 (CD)
295 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
296 ! Positive shear yield stress
297 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
298 ! Compressive yield stress in direction 1 (MD)
299 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
300 ! Compressive yield stress in direction 2 (CD)
301 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
302 ! Out-of-plane yield stress
303 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
304 ! transverse shear yield stress
305 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
306 ENDDO
307 ENDIF
308c
309 !========================================================================
310 ! - COMPUTATION OF TRIAL VALUES
311 !========================================================================
312 DO i=1,nel
313c
314 ! Computation of the trial stress tensor
315 ! - in-plane stresses
316 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
317 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
318 signxy(i) = sigoxy(i) + depsxy(i)*g12
319 ! - out-of-plane stress
320 IF (eezz(i) >= zero) THEN
321 signzz(i) = sigozz(i) + young3*depszz(i)
322 ELSE
323 signzz(i) = sigozz(i) + e3c*cc*exp(-cc*eezz(i))*depszz(i)
324 ENDIF
325 ! - transverse shear stresses
326 signyz(i) = sigoyz(i) + depsyz(i)*g23
327 signzx(i) = sigozx(i) + depszx(i)*g31
328c
329 ! Computation of trial khi coefficients
330 khi1(i) = zero
331 khi2(i) = zero
332 khi3(i) = zero
333 khi4(i) = zero
334 khi5(i) = zero
335 khi6(i) = zero
336 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
337 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
338 IF (n3*signxy(i) > zero) khi3(i) = one
339 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
340 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
341 IF (n6*signxy(i) > zero) khi6(i) = one
342c
343 ! Computation of the in-plane yield function
344 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
345 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
346 gam3(i) = n3*signxy(i)/sigy3(i)
347 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
348 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
349 gam6(i) = n6*signxy(i)/sigy3(i)
350 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
351 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
352 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
353c
354 ! Computation of the out-of-plane shear yield function
355 theta(i) = - signzz(i) - sigyo(i)
356c
357 ! Computation of the transverse shear yield function
358 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/(sigys(i))) - one
359c
360 ENDDO
361c
362 !========================================================================
363 ! - CHECKING THE PLASTIC BEHAVIOR
364 !========================================================================
365c
366 ! checking plastic behavior
367 nindx = 0
368 nindx2 = 0
369 nindx3 = 0
370 DO i=1,nel
371 ! In plane
372 IF (phi(i) > zero .AND. off(i) == one) THEN
373 nindx=nindx+1
374 index(nindx)=i
375 ENDIF
376 ! Out-of-plane
377 IF (theta(i) > zero .AND. off(i) == one) THEN
378 nindx2=nindx2+1
379 index2(nindx2)=i
380 ENDIF
381 ! Transverse shear
382 IF (psi(i) > zero .AND. off(i) == one) THEN
383 nindx3=nindx3+1
384 index3(nindx3)=i
385 ENDIF
386 ENDDO
387c
388 !====================================================================
389 ! - I IN-PLANE PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
390 !====================================================================
391#include "vectorize.inc"
392 DO ii=1,nindx
393c
394 ! Number of the element with plastic behaviour
395 i = index(ii)
396c
397 ! Note : in this part, the purpose is to compute for each iteration
398 ! a plastic multiplier allowing to update internal variables to satisfy
399 ! the consistency condition.
400 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
401 ! -> PHI_OLD : old value of yield function (known)
402 ! -> DPHI : yield function prediction (to compute)
403 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
404 ! into account of internal variables kinetic :
405 ! plasticity, temperature and damage (to compute)
406c
407 ! 1 - Computation of DPHI_DSIG the normal to the yield criterion
408 !-------------------------------------------------------------
409 ! Computation of derivatives to the yield criterion
410 ! - in-plane normal
411 gam1(i) = (n1(1)*sigoxx(i)+n1(2)*sigoyy(i))/sigy1(i)
412 gam2(i) = (n2(1)*sigoxx(i)+n2(2)*sigoyy(i))/sigy2(i)
413 gam3(i) = n3*sigoxy(i)/sigy3(i)
414 gam4(i) = (n4(1)*sigoxx(i)+n4(2)*sigoyy(i))/sigy4(i)
415 gam5(i) = (n5(1)*sigoxx(i)+n5(2)*sigoyy(i))/sigy5(i)
416 gam6(i) = n6*sigoxy(i)/sigy3(i)
417 normxx = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(1)/sigy1(i)) +
418 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(1)/sigy2(i)) +
419 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(1)/sigy4(i)) +
420 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(1)/sigy5(i))
421 normyy = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(2)/sigy1(i)) +
422 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(2)/sigy2(i)) +
423 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(2)/sigy4(i)) +
424 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(2)/sigy5(i))
425 normxy = deuxk/two*khi3(i)*(gam3(i)**(deuxk-1))*(n3/sigy3(i)) +
426 . deuxk/two*khi6(i)*(gam6(i)**(deuxk-1))*(n6/sigy3(i))
427c
428 ! Normalization of the derivatives
429 ! - in-plane normal
430 normsig = sqrt(normxx**2 + normyy**2 + two*normxy**2)
431 normsig = max(normsig,em20)
432 normpxx = normxx/normsig
433 normpyy = normyy/normsig
434 normpxy = normxy/normsig
435c
436 ! 3 - Computation of DPHI_DLAMBDA
437 !---------------------------------------------------------
438c
439 ! a) Derivative with respect stress increments tensor DSIG
440 ! --------------------------------------------------------
441 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
442 . + normyy * (a21*normpxx + a22*normpyy)
443 . + two*normxy * two*normpxy * g12
444c
445 ! b) Derivatives with respect to hardening parameters
446 ! ---------------------------------------------------
447c
448 ! i) Derivatives of PHI with respect to the yield stresses
449 dphi_dsigy1 = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(-(n1(1)*sigoxx(i)+n1(2)*sigoyy(i))/(sigy1(i)**2))
450 dphi_dsigy2 = deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(-(n2(1)*sigoxx(i)+n2(2)*sigoyy(i))/(sigy2(i)**2))
451 dphi_dsigy3 = deuxk*khi3(i)*(gam3(i)**(deuxk-1))*(-n3*sigoxy(i)/(sigy3(i)**2)) +
452 . deuxk*khi6(i)*(gam6(i)**(deuxk-1))*(-n6*sigoxy(i)/(sigy3(i)**2))
453 dphi_dsigy4 = deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(-(n4(1)*sigoxx(i)+n4(2)*sigoyy(i))/(sigy4(i)**2))
454 dphi_dsigy5 = deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(-(n5(1)*sigoxx(i)+n5(2)*sigoyy(i))/(sigy5(i)**2))
455 ! ii) Derivatives of yield stresses with respect to hardening parameters
456 IF (itab == 0) THEN
457 dsigy1_dp(i) = a01*b01*(one-(tanh(b01*pla(i,2)))**2) + c01
458 dsigy2_dp(i) = a02*b02*(one-(tanh(b02*pla(i,2)))**2) + c02
459 dsigy3_dp(i) = a03*b03*(one-(tanh(b03*pla(i,2)))**2) + c03
460 dsigy4_dp(i) = a04*b04*(one-(tanh(b04*pla(i,2)))**2) + c04
461 dsigy5_dp(i) = a05*b05*(one-(tanh(b05*pla(i,2)))**2) + c05
462 ENDIF
463 ! iii) Assembling derivatives of PHI with respect to hardening parameter
464 hardp(i) = sqrt(khi1(i)*dsigy1_dp(i)**2 + khi2(i)*dsigy2_dp(i)**2 +
465 . khi3(i)*two*dsigy3_dp(i)**2 + khi4(i)*dsigy4_dp(i)**2 +
466 . khi5(i)*dsigy5_dp(i)**2)
467 dphi_dpla = dphi_dsigy1*dsigy1_dp(i) + dphi_dsigy2*dsigy2_dp(i) +
468 . dphi_dsigy3*dsigy3_dp(i) + dphi_dsigy4*dsigy4_dp(i) +
469 . dphi_dsigy5*dsigy5_dp(i)
470c
471 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
472 dphi_dlam = -dfdsig2 + dphi_dpla
473 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
474c
475 ! 4 - Computation of plastic multiplier and variables update
476 !----------------------------------------------------------
477c
478 ! a) Restoring previous value of the yield function
479 phi(i) = phi0(i)
480c
481 ! b) Computation of the trial stress increment
482 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
483 dsigyy(i) = a21*depsxx(i) + a22*depsyy(i)
484 dsigxy(i) = depsxy(i)*g12
485c
486 ! c) Computation of yield surface trial increment DPHI
487 dphi = normxx * dsigxx(i)
488 . + normyy * dsigyy(i)
489 . + two * normxy * dsigxy(i)
490c
491 ! d) Computation of the plastic multiplier
492 dlam = -(phi(i) + dphi) / dphi_dlam
493 dlam = max(dlam, zero)
494c
495 ! e) Plastic strains tensor update
496 dpxx = dlam * normpxx
497 dpyy = dlam * normpyy
498 dpxy = two * dlam * normpxy
499c
500 ! f) Elasto-plastic stresses update
501 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
502 signyy(i) = sigoyy(i) + dsigyy(i) - (a21*dpxx + a22*dpyy)
503 signxy(i) = sigoxy(i) + dsigxy(i) - g12*dpxy
504c
505 ! g) Cumulated plastic strain and hardening parameter update
506 ddep = dlam
507 dpla2(i) = max(zero, dpla2(i) + ddep)
508 pla(i,2) = pla(i,2) + ddep
509c
510 ! h) Yield stresses update
511 IF (itab == 0) THEN
512 ! Tensile yield stress in direction 1 (MD)
513 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
514 ! Tensile yield stress in direction 2 (CD)
515 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
516 ! Positive shear yield stress
517 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
518 ! Compressive yield stress in direction 1 (MD)
519 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
520 ! Compressive yield stress in direction 2 (CD)
521 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
522 ENDIF
523c
524 ! i) Update of trial alpha coefficients
525 khi1(i) = zero
526 khi2(i) = zero
527 khi3(i) = zero
528 khi4(i) = zero
529 khi5(i) = zero
530 khi6(i) = zero
531 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
532 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
533 IF (n3*signxy(i) > zero) khi3(i) = one
534 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
535 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
536 IF (n6*signxy(i) > zero) khi6(i) = one
537c
538 ! j) Yield function value update
539 IF (itab == 0) THEN
540 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
541 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
542 gam3(i) = n3*signxy(i)/sigy3(i)
543 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
544 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
545 gam6(i) = n6*signxy(i)/sigy3(i)
546 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
547 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
548 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
549 ENDIF
550c
551 ENDDO
552 !===================================================================
553 ! - END OF IN-PLANE PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
554 !===================================================================
555c
556 !====================================================================
557 ! - II OUT-OF-PLANE PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
558 !====================================================================
559#include "vectorize.inc"
560 DO ii=1,nindx2
561c
562 ! Number of the element with plastic behaviour
563 i = index2(ii)
564c
565 ! Note : in this part, the purpose is to compute for each iteration
566 ! a plastic multiplier allowing to update internal variables to satisfy
567 ! the consistency condition.
568 ! Its expression is : DLAMBDA = - (THETA + DTHETA) / DPHI_DLAMBDA
569 ! -> THETA : current value of yield function (known)
570 ! -> DTHETA : yield function prediction (to compute)
571 ! -> DTHETA_DLAMBDA : derivative of THETA with respect to DLAMBDA by taking
572 ! into account of internal variables kinetic :
573 ! plasticity, temperature and damage (to compute)
574c
575 ! 1 - Computation of DTHETA_DSIG the normal to the yield criterion
576 !-------------------------------------------------------------
577 ! Computation of derivatives to the yield criterion
578 ! - in-plane normal
579 normzz = -one
580c
581 ! 3 - Computation of DTHETA_DLAMBDA
582 !---------------------------------------------------------
583c
584 ! a) Derivative with respect stress increments tensor DSIG
585 ! --------------------------------------------------------
586 IF (eezz(i) >= zero) THEN
587 dfdsig2 = normzz * young3 * normzz
588 ELSE
589 dfdsig2 = normzz * e3c*cc*exp(-cc*eezz(i)) * normzz
590 ENDIF
591c
592 ! b) Derivatives with respect to hardening parameters
593 ! ---------------------------------------------------
594c
595 ! i) Derivatives of THETA with respect to the yield stresses
596 dtheta_dsigyo = -one
597 ! ii) Derivatives of yield stresses with respect to hardening parameters
598 IF (itab == 0) dsigyo_dp(i) = csig*bsig*exp(csig*pla(i,3))
599 ! iii) Assembling derivatives of THETA with respect to hardening parameter
600 dtheta_dpla = dtheta_dsigyo*dsigyo_dp(i)
601c
602 ! iv) Derivative of THETA with respect to DLAM ( = -DENOM)
603 dtheta_dlam = -dfdsig2 + dtheta_dpla
604 dtheta_dlam = sign(max(abs(dtheta_dlam),em20) ,dtheta_dlam)
605c
606 ! 4 - Computation of plastic multiplier and variables update
607 !----------------------------------------------------------
608c
609 ! a) Restoring previous value of the yield function
610 theta(i) = theta0(i)
611c
612 ! b) Computation of the trial stress increment
613 IF (eezz(i) >= zero) THEN
614 dsigzz(i) = young3*depszz(i)
615 ELSE
616 dsigzz(i) = e3c*cc*exp(-cc*eezz(i))*depszz(i)
617 ENDIF
618c
619 ! c) Computation of yield surface trial increment DTHETA
620 dtheta = normzz * dsigzz(i)
621c
622 ! d) Computation of the plastic multiplier
623 dlam = -(theta(i) + dtheta) / dtheta_dlam
624 dlam = max(dlam, zero)
625c
626 ! e) Elasto-plastic stresses update
627 IF (eezz(i)>=zero) THEN
628 signzz(i) = sigozz(i) + dsigzz(i) - young3*dlam*normzz
629 ELSE
630 signzz(i) = sigozz(i) + dsigzz(i) - e3c*cc*exp(-cc*eezz(i))*dlam*normzz
631 ENDIF
632c
633 ! f) Cumulated plastic strain and hardening parameter update
634 ddep = dlam
635 dpla3(i) = max(zero,dpla3(i) + ddep)
636 pla(i,3) = pla(i,3) + ddep
637 eezz(i) = epszz(i) + pla(i,3)
638c
639 ! g) Yield stresses update
640 IF (itab == 0) THEN
641 ! Out-of-plane yield stress
642 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
643 ! Yield function value update
644 theta(i) = - signzz(i) - sigyo(i)
645 ENDIF
646c
647 ENDDO
648 !=======================================================================
649 ! - END OF OUT-OF-PLANE PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
650 !=======================================================================
651c
652 !========================================================================
653 ! - III TRANSVERSE SHEAR PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
654 !========================================================================
655#include "vectorize.inc"
656 DO ii=1,nindx3
657c
658 ! Number of the element with plastic behaviour
659 i = index3(ii)
660c
661 ! Note : in this part, the purpose is to compute for each iteration
662 ! a plastic multiplier allowing to update internal variables to satisfy
663 ! the consistency condition.
664 ! Its expression is : DLAMBDA = - (PSI + DPSI) / DPSI_DLAMBDA
665 ! -> PSI : current value of yield function (known)
666 ! -> DPSI : yield function prediction (to compute)
667 ! -> DPSI_DLAMBDA : derivative of PSI with respect to DLAMBDA by taking
668 ! into account of internal variables kinetic :
669 ! plasticity, temperature and damage (to compute)
670c
671 ! 1 - Computation of DPSI_DSIG the normal to the yield criterion
672 !-------------------------------------------------------------
673c
674 ! Computation of derivatives to the yield criterion
675 ! - transverse shear normal
676 normyz = sigoyz(i)/max((sqrt(sigoyz(i)**2 + sigozx(i)**2)*sigys(i)),em20)
677 normzx = sigozx(i)/max((sqrt(sigoyz(i)**2 + sigozx(i)**2)*sigys(i)),em20)
678c
679 ! 2 - Computation of DPSI_DLAMBDA
680 !---------------------------------------------------------
681c
682 ! a) Derivative with respect stress increments tensor DSIG
683 ! --------------------------------------------------------
684 dfdsig2 = two*normyz * g23 * two*normyz +
685 . two*normzx * g31 * two*normzx
686c
687 ! b) Derivatives with respect to hardening parameters
688 ! ---------------------------------------------------
689c
690 ! i) Derivatives of PSI with respect to the yield stresses
691 dpsi_dsigys = -sqrt(sigoyz(i)**2 + sigozx(i)**2)/(sigys(i)**2)
692 ! ii) Derivatives of yield stresses with respect to hardening parameters
693 IF (itab == 0) dsigys_dp(i) = atau - min(zero,sigozz(i))*btau
694 ! iii) Assembling derivatives of PSI with respect to hardening parameter
695 dpsi_dpla = dpsi_dsigys*dsigys_dp(i)
696c
697 ! iv) Derivative of PSI with respect to DLAM ( = -DENOM)
698 dpsi_dlam = -dfdsig2 + dpsi_dpla
699 dpsi_dlam = sign(max(abs(dpsi_dlam),em20),dpsi_dlam)
700c
701 ! 4 - Computation of plastic multiplier and variables update
702 !----------------------------------------------------------
703c
704 ! a) Restoring previous value of the yield function
705 psi(i) = psi0(i)
706c
707 ! b) Computation of the trial stress increment
708 dsigyz(i) = depsyz(i)*g23
709 dsigzx(i) = depszx(i)*g31
710c
711 ! c) Computation of yield surface trial increment DPSI
712 dpsi = two*normyz * dsigyz(i)
713 . + two*normzx * dsigzx(i)
714c
715 ! d) Computation of the plastic multiplier
716 dlam = -(psi(i) + dpsi) / dpsi_dlam
717 dlam = max(dlam, zero)
718c
719 ! e) Plastic strains tensor update
720 dpyz = two*dlam * normyz
721 dpzx = two*dlam * normzx
722c
723 ! f) Elasto-plastic stresses update
724 signyz(i) = sigoyz(i) + dsigyz(i) - g23*dpyz
725 signzx(i) = sigozx(i) + dsigzx(i) - g31*dpzx
726c
727 ! g) Cumulated plastic strain and hardening parameter update
728 ddep = dlam
729 dpla4(i) = max(zero, dpla4(i) + ddep)
730 pla(i,4) = pla(i,4) + ddep
731c
732 ! h) Yield stresses update
733 IF (itab == 0) THEN
734 ! Transverse shear yield stress
735 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
736 ! i) Yield function value update
737 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
738 ENDIF
739c
740 ENDDO
741 !===========================================================================
742 ! - END OF TRANSVERSE SHEAR PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
743 !===========================================================================
744c
745 ! If tabulated yield function, update of the yield stress for all element
746 IF (itab > 0) THEN
747 IF (nindx > 0) THEN
748 ! In-plane tabulation with strain-rate
749 xvec(1:nel,1) = pla(1:nel,2)
750 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
751 ! -> Tensile yield stress in direction 1 (MD)
752 ipos(1:nel,1) = vartmp(1:nel,1)
753 ipos(1:nel,2) = 1
754 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
755 sigy1(1:nel) = sigy1(1:nel) * yscale1
756 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
757 vartmp(1:nel,1) = ipos(1:nel,1)
758 ! -> Tensile yield stress in direction 2 (CD)
759 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
760 ipos(1:nel,1) = vartmp(1:nel,2)
761 ipos(1:nel,2) = 1
762 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
763 sigy2(1:nel) = sigy2(1:nel) * yscale2
764 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
765 vartmp(1:nel,2) = ipos(1:nel,1)
766 ! -> Positive shear yield stress
767 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
768 ipos(1:nel,1) = vartmp(1:nel,3)
769 ipos(1:nel,2) = 1
770 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
771 sigy3(1:nel) = sigy3(1:nel) * yscale3
772 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
773 vartmp(1:nel,3) = ipos(1:nel,1)
774 ! -> Compressive yield stress in direction 1 (MD)
775 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
776 ipos(1:nel,1) = vartmp(1:nel,4)
777 ipos(1:nel,2) = 1
778 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
779 sigy4(1:nel) = sigy4(1:nel) * yscale4
780 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
781 vartmp(1:nel,4) = ipos(1:nel,1)
782 ! -> Compressive yield stress in direction 2 (CD)
783 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
784 ipos(1:nel,1) = vartmp(1:nel,5)
785 ipos(1:nel,2) = 1
786 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
787 sigy5(1:nel) = sigy5(1:nel) * yscale5
788 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
789 vartmp(1:nel,5) = ipos(1:nel,1)
790 ! Yield function value update
791#include "vectorize.inc"
792 DO ii=1,nindx
793 i = index(ii)
794 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
795 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
796 gam3(i) = n3*signxy(i)/sigy3(i)
797 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
798 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
799 gam6(i) = n6*signxy(i)/sigy3(i)
800 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
801 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
802 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
803 ENDDO
804 ENDIF
805 IF (nindx2 > 0) THEN
806 ! -> Transverse shear tabulation with strain-rate
807 xvec(1:nel,1) = pla(1:nel,3)
808 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
809 ! -> Transverse shear yield stress
810 ipos(1:nel,1) = vartmp(1:nel,6)
811 ipos(1:nel,2) = 1
812 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
813 sigyo(1:nel) = sigyo(1:nel) * yscalec
814 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
815 vartmp(1:nel,6) = ipos(1:nel,1)
816 ! Yield function value update
817#include "vectorize.inc"
818 DO ii=1,nindx2
819 i = index2(ii)
820 theta(i) = - signzz(i) - sigyo(i)
821 ENDDO
822 ENDIF
823 IF (nindx3 > 0) THEN
824 ! Transverse shear tabulation with strain-rate
825 xvec(1:nel,1) = pla(1:nel,4)
826 xvec(1:nel,2) = epsd(1:nel,4) * xscales
827 ! Transverse shear yield stress
828 ipos(1:nel,1) = vartmp(1:nel,7)
829 ipos(1:nel,2) = 1
830 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
831 sigys(1:nel) = sigys(1:nel) * yscales
832 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
833 vartmp(1:nel,7) = ipos(1:nel,1)
834 ! Yield function value update
835#include "vectorize.inc"
836 DO ii=1,nindx3
837 i = index3(ii)
838 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
839 ENDDO
840 ENDIF
841 ENDIF
842c
843 ! Storing new values
844 DO i=1,nel
845 ! USR Outputs
846 uvar(i,1) = zero ! reset in-plane yield function value
847 uvar(i,2) = zero ! Reset out-of-plane yield function value
848 uvar(i,3) = zero ! Reset transverse shear yield function value
849 ! Plastic strain
850 pla(i,1) = sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2)
851 dpla(i) = (pla(i,2)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla2(i) +
852 . (pla(i,3)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla3(i) +
853 . (pla(i,4)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla4(i)
854 ! Plastic strain-rate
855 IF (itab == 0) THEN
856 epsd(i,1) = dpla(i) / max(em20,timestep)
857 epsd(i,2) = dpla2(i) / max(em20,timestep)
858 epsd(i,3) = dpla3(i) / max(em20,timestep)
859 epsd(i,4) = dpla4(i) / max(em20,timestep)
860 ELSE
861 dpdt = dpla(i)/max(em20,timestep)
862 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
863 dpdt = dpla2(i)/max(em20,timestep)
864 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
865 dpdt = dpla3(i)/max(em20,timestep)
866 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
867 dpdt = dpla4(i)/max(em20,timestep)
868 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
869 ENDIF
870 ! Coefficient for hourglass and soundspeed
871 IF (eezz(i) >= zero) THEN
872 IF (dpla(i) > zero) THEN
873 et(i) = hardp(i) / (hardp(i) + max(young1,young2,young3))
874 ELSE
875 et(i) = one
876 ENDIF
877 soundsp(i) = sqrt(max(a11,a12,a21,a22,young3,g12,g23,g31)/ rho(i))
878 ELSE
879 IF (dpla(i) > zero) THEN
880 et(i) = hardp(i) / (hardp(i) + max(young1,young2,e3c*cc*exp(-cc*eezz(i))))
881 ELSE
882 et(i) = one
883 ENDIF
884 soundsp(i) = sqrt(max(a11,a12,a21,a22,e3c*cc*exp(-cc*eezz(i)),g12,g23,g31)/ rho(i))
885 ENDIF
886 ! Storing the yield stress
887 sigy(i) = sqrt(khi1(i)*sigy1(i)**2 + khi2(i)*sigy2(i)**2 + khi4(i)*sigy4(i)**2 +
888 . khi5(i)*sigy5(i)**2 + two*khi3(i)*sigy3(i)**2 + sigyo(i)**2
889 . + two*sigys(i)**2)
890 ENDDO
891c
892 ! Save in-plane yield function remaining value
893#include "vectorize.inc"
894 DO ii=1,nindx
895 i = index(ii)
896 uvar(i,1) = phi(i)
897 ENDDO
898c
899 ! Save out-of-plane yield function remaining value
900#include "vectorize.inc"
901 DO ii=1,nindx2
902 i = index2(ii)
903 uvar(i,2) = theta(i)
904 ENDDO
905c
906 ! Save transverse shear yield function remaining value
907#include "vectorize.inc"
908 DO ii=1,nindx3
909 i = index3(ii)
910 uvar(i,3) = psi(i)
911 ENDDO
912c
my_real
#define my_real
Definition cppsort.cpp:32
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112