OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112_xia_nice.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| mat112_xia_nice ../engine/source/materials/mat/mat112/mat112_xia_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps112 ../engine/source/materials/mat/mat112/sigeps112.F
27!||--- calls -----------------------------------------------------
28!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.F
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.F
31!|| table_mod ../engine/share/modules/table_mod.F
32!||====================================================================
33 SUBROUTINE mat112_xia_nice(
34 1 NEL ,NGL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP,
35 2 UPARAM ,UVAR ,JTHE ,OFF ,RHO0 ,RHO ,
36 3 PLA ,DPLA ,EPSD ,SOUNDSP ,EPSZZ ,
37 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 7 SIGY ,ET ,
41 8 NVARTMP ,NUMTABL ,VARTMP ,ITABLE ,TABLE )
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
913 END
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
mat112_xia_nice
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)
Definition mat112_xia_nice.F:42
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
table_mod::ttable
Definition table_mod.F:125