OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat107c_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!|| mat107c_nice ../engine/source/materials/mat/mat107/mat107c_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps107c ../engine/source/materials/mat/mat107/sigeps107c.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 mat107c_nice(
34 1 NEL ,NGL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP,
35 2 UPARAM ,UVAR ,JTHE ,OFF ,RHO ,
36 3 PLA ,DPLA ,EPSD ,SOUNDSP ,SHF ,
37 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 7 SIGY ,ET ,
41 8 NUMTABL ,ITABLE ,TABLE ,NVARTMP ,VARTMP )
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,shf,
67 . depsxx,depsyy,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL), INTENT(OUT) ::
71 . soundsp,sigy,et,
72 . signxx,signyy,signxy,signyz,signzx
73c
74 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
75 . dpla,off
76 my_real ,DIMENSION(NEL,6),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),ITAB,ISMOOTH,IPOS(NEL,2)
86c
87 my_real
88 . young1,young2,nu12,nu21,g12,a11,a12,a21,a22,c1,
89 . xi1,xi2,k1,k2,k3,k4,k5,k6,g1c,d1,d2,sigy1,
90 . cini1,s1,sigy2,cini2,s2,sigy1c,cini1c,s1c,
91 . sigy2c,cini2c,s2c,sigyt,cinit,st,g23,g31
92 my_real
93 . normsig,beta1,beta2,
94 . dpdt,dlam,ddep
95 my_real
96 . normxx,normyy,normxy,
97 . normpxx,normpyy,normpxy,
98 . dphi_dlam,dphi,dfdsig2,dphi_dpla,dpxx,dpyy,dpxy,
99 . dphi_dr1,dphi_dr2,dphi_dr1c,dphi_dr2c,dphi_drt
100 my_real
101 . xscale1,yscale1,xscale1c,yscale1c,xscale2,yscale2,xscale2c,yscale2c,
102 . xscalet,yscalet,xvec(nel,2),asrate
103c
104 my_real, DIMENSION(NEL) ::
105 . alpha1,alpha2,alpha3,alpha4,alpha5,r1,r2,r1c,r2c,rt,
106 . dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,hardp,phi,phi0,eta1,eta2,
107 . dr1_dp,dr2_dp,dr1c_dp,dr2c_dp,drt_dp,hardr,dpla2,dpla3,dpla4,
108 . dpla5,dpla6
109c
110 !=======================================================================
111 ! - INITIALISATION OF COMPUTATION ON TIME STEP
112 !=======================================================================
113 ! Recovering model parameters
114 ! Elastic parameters
115 young1 = uparam(1) ! Young modulus in direction 1 (MD)
116 young2 = uparam(2) ! Young modulus in direction 2 (CD)
117 nu12 = uparam(4) ! Poisson's ratio in 12
118 nu21 = uparam(5) ! Poisson's ratio in 21
119 a11 = uparam(6) ! Component 11 of orthotropic shell elasticity matrix
120 a12 = uparam(7) ! Component 12 of orthotropic shell elasticity matrix
121 a21 = uparam(8) ! Component 21 of orthotropic shell elasticity matrix
122 a22 = uparam(9) ! Component 22 of orthotropic shell elasticity matrix
123 g12 = uparam(10) ! Shear modulus in 12
124 g23 = uparam(11) ! Shear modulus in 23
125 g31 = uparam(12) ! Shear modulus in 31
126 itab = int(uparam(14)) ! Tabulated yield stress flag
127 xi1 = uparam(15) ! 1st coupling parameter
128 xi2 = uparam(16) ! 2nd coupling parameter
129 k1 = uparam(17) ! 1st plastic potential parameter
130 k2 = uparam(18) ! 2nd plastic potential parameter
131 k3 = uparam(19) ! 3rd plastic potential parameter
132 k4 = uparam(20) ! 4th plastic potential parameter
133 k5 = uparam(21) ! 5th plastic potential parameter
134 k6 = uparam(22) ! 6th plastic potential parameter
135 g1c = uparam(23) ! Correction factor for R1C
136 IF (itab == 0) THEN
137 d1 = uparam(24) ! 1st Shear yield stress parameter
138 d2 = uparam(25) ! 2nd Shear yield stress parameter
139 sigy1 = uparam(26) ! Initial yield stress in tension in direction 1 (MD)
140 cini1 = uparam(27) ! Yield stress intersection with ordinate axis in tension in direction 1 (MD)
141 s1 = uparam(28) ! Yield stress slope in tension in direction 1 (MD)
142 sigy2 = uparam(29) ! Initial yield stress in tension in direction 2 (CD)
143 cini2 = uparam(30) ! Yield stress intersection with ordinate axis in tension in direction 2 (CD)
144 s2 = uparam(31) ! Yield stress slope in tension in direction 2 (CD)
145 sigy1c = uparam(32) ! Initial yield stress in compression in direction 1 (MD)
146 cini1c = uparam(33) ! Yield stress intersection with ordinate axis in compression in direction 1 (MD)
147 s1c = uparam(34) ! Yield stress slope in compression in direction 1 (MD)
148 sigy2c = uparam(35) ! Initial yield stress in compression in direction 2 (CD)
149 cini2c = uparam(36) ! Yield stress intersection with ordinate axis in compression in direction 2 (CD)
150 s2c = uparam(37) ! Yield stress slope in compression in direction 2 (CD)
151 sigyt = uparam(38) ! Initial yield stress in shear
152 cinit = uparam(39) ! Yield stress intersection with ordinate axis in shear
153 st = uparam(40) ! Yield stress slope in shear
154 ELSE
155 xscale1 = uparam(24)
156 yscale1 = uparam(25)
157 xscale2 = uparam(26)
158 yscale2 = uparam(27)
159 xscale1c= uparam(28)
160 yscale1c= uparam(29)
161 xscale2c= uparam(30)
162 yscale2c= uparam(31)
163 xscalet = uparam(32)
164 yscalet = uparam(33)
165 asrate = uparam(34)
166 asrate = (two*pi*asrate*timestep)/(two*pi*asrate*timestep + one)
167 ismooth = int(uparam(35))
168 ENDIF
169c
170 ! Recovering internal variables
171 DO i=1,nel
172 ! OFF parameter for element deletion
173 IF (off(i) < 0.1) off(i) = zero
174 IF (off(i) < one) off(i) = off(i)*four_over_5
175 ! User inputs
176 phi0(i) = uvar(i,1) ! Previous yield function value
177 ! Standard inputs
178 dpla(i) = zero ! Initialization of the "global" plastic strain increment
179 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
180 dpla3(i) = zero ! Initialization of the transverse shear plastic strain increment
181 dpla4(i) = zero ! Initialization of the in-plane plastic strain increment
182 dpla5(i) = zero ! Initialization of the transverse shear plastic strain increment
183 dpla6(i) = zero ! Initialization of the in-plane plastic strain increment
184 et(i) = one ! Initialization of hourglass coefficient
185 hardp(i) = zero ! Initialization of hourglass coefficient
186 ENDDO
187c
188 ! Computation of the initial yield stress
189 IF (itab > 0) THEN
190 ! -> Tensile yield stress in direction 1 (MD)
191 xvec(1:nel,1) = pla(1:nel,2)
192 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
193 ipos(1:nel,1) = vartmp(1:nel,1)
194 ipos(1:nel,2) = 1
195 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
196 r1(1:nel) = r1(1:nel) * yscale1
197 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
198 vartmp(1:nel,1) = ipos(1:nel,1)
199 ! -> Compressive yield stress in direction 1 (MD)
200 xvec(1:nel,1) = pla(1:nel,3)
201 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
202 ipos(1:nel,1) = vartmp(1:nel,2)
203 ipos(1:nel,2) = 1
204 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
205 r2(1:nel) = r2(1:nel) * yscale2
206 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
207 vartmp(1:nel,2) = ipos(1:nel,1)
208 ! -> Positive shear yield stress
209 xvec(1:nel,1) = pla(1:nel,4)
210 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
211 ipos(1:nel,1) = vartmp(1:nel,3)
212 ipos(1:nel,2) = 1
213 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
214 r1c(1:nel) = r1c(1:nel) * yscale1c
215 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
216 vartmp(1:nel,3) = ipos(1:nel,1)
217 ! -> Tensile yield stress in direction 2 (CD)
218 xvec(1:nel,1) = pla(1:nel,5)
219 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
220 ipos(1:nel,1) = vartmp(1:nel,4)
221 ipos(1:nel,2) = 1
222 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
223 r2c(1:nel) = r2c(1:nel) * yscale2c
224 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
225 vartmp(1:nel,4) = ipos(1:nel,1)
226 ! Compressive yield stress in direction 2 (CD)
227 xvec(1:nel,1) = pla(1:nel,6)
228 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
229 ipos(1:nel,1) = vartmp(1:nel,5)
230 ipos(1:nel,2) = 1
231 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
232 rt(1:nel) = rt(1:nel) * yscalet
233 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
234 vartmp(1:nel,5) = ipos(1:nel,1)
235 ENDIF
236c
237 !========================================================================
238 ! - COMPUTATION OF TRIAL VALUES
239 !========================================================================
240 DO i=1,nel
241c
242 ! Computation of the trial stress tensor
243 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
244 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
245 signxy(i) = sigoxy(i) + depsxy(i)*g12
246 signyz(i) = sigoyz(i) + depsyz(i)*g23*shf(i)
247 signzx(i) = sigozx(i) + depszx(i)*g31*shf(i)
248c
249 ! Computation of trial alpha coefficients
250 alpha1(i) = zero
251 alpha2(i) = zero
252 alpha3(i) = zero
253 alpha4(i) = zero
254 alpha5(i) = zero
255 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
256 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
257 IF (-signxx(i) > zero) alpha3(i) = one
258 IF (-signyy(i) > zero) alpha4(i) = one
259 IF (abs(signxy(i)) > zero) alpha5(i) = one
260c
261 IF (itab == 0) THEN
262 ! Computation of yield stresses
263 ! Tensile yield stresses
264 ! - in direction 1 (MD)
265 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
266 ! - in direction 2 (CD)
267 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
268 ! Compressive yield stresses
269 ! - in direction 1 (MD) (including correction)
270 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
271 r1c(i) = r1c(i)/sqrt(one - g1c)
272 ! - in direction 2 (CD)
273 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
274 ! Shear yield stress
275 eta1(i) = one + alpha3(i)*alpha4(i)*d1
276 eta2(i) = one + alpha3(i)*alpha4(i)*d2
277 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
278 ENDIF
279c
280 ! Computation of the yield function
281 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
282 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
283 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
284 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
285 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
286 ENDDO
287c
288 !========================================================================
289 ! - CHECKING THE PLASTIC BEHAVIOR
290 !========================================================================
291c
292 ! Checking plastic behavior for all elements
293 nindx = 0
294 DO i=1,nel
295 IF (phi(i) > zero .AND. off(i) == one) THEN
296 nindx=nindx+1
297 index(nindx)=i
298 ENDIF
299 ENDDO
300c
301 !====================================================================
302 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
303 !====================================================================
304#include "vectorize.inc"
305 DO ii=1,nindx
306c
307 ! Number of the element with plastic behaviour
308 i = index(ii)
309c
310 ! Note : in this part, the purpose is to compute for each iteration
311 ! a plastic multiplier allowing to update internal variables to satisfy
312 ! the consistency condition.
313 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
314 ! -> PHI : current value of yield function (known)
315 ! -> DPHI : yield function prediction (to compute)
316 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
317 ! into account of internal variables kinetic :
318 ! plasticity, temperature and damage (to compute)
319c
320 ! 1 - Computation of DPHI_DSIG the normal to the yield criterion
321 !-------------------------------------------------------------
322c
323 normxx = two*(alpha1(i)/(r1(i)**2))*(sigoxx(i)-xi1*sigoyy(i))
324 . - two*(alpha2(i)*xi2/(r2(i)**2))*(sigoyy(i)-xi2*sigoxx(i))
325 . + two*(alpha3(i)/(r1c(i)**2))*sigoxx(i)
326 normyy = -two*(alpha1(i)*xi1/(r1(i)**2))*(sigoxx(i)-xi1*sigoyy(i))
327 . + two*(alpha2(i)/(r2(i)**2))*(sigoyy(i)-xi2*sigoxx(i))
328 . + two*(alpha4(i)/(r2c(i)**2))*sigoyy(i)
329 normxy = (alpha5(i)/(rt(i)**2))*abs(sigoxy(i))*sign(one,sigoxy(i))
330c
331 ! 2 - Computation of DFP_DSIG the normal to the non-associated plastic potential
332 !-------------------------------------------------------------------------------
333 ! a) Computation of BETA switching coefficient
334 normsig = sqrt(sigoxx(i)*sigoxx(i) + sigoyy(i)*sigoyy(i) + two*sigoxy(i)*sigoxy(i))
335 normsig = max(normsig,em20)
336 beta1 = zero
337 beta2 = zero
338 IF ((sigoxx(i)/normsig > em01).OR.(sigoyy(i)/normsig > em01)) beta1 = one
339 IF ((sigoxx(i)/normsig < -em01).OR.(sigoyy(i)/normsig < -em01)) beta2 = one
340 IF (((abs(sigoxx(i))/normsig)<em01).AND.((abs(sigoyy(i))/normsig)<em01)) THEN
341 beta1 = one
342 beta2 = one
343 ENDIF
344 ! b) Computation of derivatives
345 normpxx = beta1*(two*sigoxx(i) + k2*sigoyy(i)) + beta2*(two*sigoxx(i) + k4*sigoyy(i))
346 normpyy = beta1*(two*k1*sigoyy(i) + k2*sigoxx(i)) + beta2*(two*k3*sigoyy(i) + k4*sigoxx(i))
347 normpxy = (beta1*k5 + beta2*k6)*sigoxy(i)
348 ! c) Computation of the norm of the derivative
349 normsig = sqrt(normpxx*normpxx + normpyy*normpyy + two*normpxy*normpxy)
350 normsig = max(normsig,em20)
351 ! d) Computation of the normal to the plastic potential
352 normpxx = normpxx/normsig
353 normpyy = normpyy/normsig
354 normpxy = normpxy/normsig
355c
356 ! 3 - Computation of DPHI_DLAMBDA
357 !---------------------------------------------------------
358c
359 ! a) Derivative with respect stress increments tensor DSIG
360 ! --------------------------------------------------------
361 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
362 . + normyy * (a21*normpxx + a22*normpyy)
363 . + two*normxy * two*normpxy * g12
364c
365 ! b) Derivatives with respect to hardening parameters
366 ! ---------------------------------------------------
367c
368 ! i) Derivatives of PHI with respect to the yield stresses
369 dphi_dr1 = -two*alpha1(i)*(((sigoxx(i)-xi1*sigoyy(i))**2)/(r1(i)**3))
370 dphi_dr2 = -two*alpha2(i)*(((sigoyy(i)-xi2*sigoxx(i))**2)/(r2(i)**3))
371 dphi_dr1c = -two*alpha3(i)*(sigoxx(i)**2)/(r1c(i)**3)
372 dphi_dr2c = -two*alpha4(i)*(sigoyy(i)**2)/(r2c(i)**3)
373 dphi_drt = -two*alpha5(i)*(sigoxy(i)**2)/(rt(i)**3)
374 ! ii) Derivatives of yield stresses with respect to hardening parameters
375 IF (itab == 0) THEN
376 dr1_dp(i) = (one/(cini1 + s1*pla(i,2))) * (one - (s1*pla(i,2) /(cini1 + s1*pla(i,2))))
377 dr2_dp(i) = (one/(cini2 + s2*pla(i,3))) * (one - (s2*pla(i,3) /(cini2 + s2*pla(i,3))))
378 dr1c_dp(i) = (one/(cini1c + s1c*pla(i,4))) * (one - (s1c*pla(i,4)/(cini1c + s1c*pla(i,4))))
379 dr1c_dp(i) = dr1c_dp(i)/sqrt(one - g1c)
380 dr2c_dp(i) = (one/(cini2c + s2c*pla(i,5))) * (one - (s2c*pla(i,5)/(cini2c + s2c*pla(i,5))))
381 drt_dp(i) = eta2(i)*(one/(cinit + st*pla(i,6))) * (one - (st*pla(i,6) /(cinit + st*pla(i,6))))
382 ENDIF
383 ! iii) Assembling derivatives of PHI with respect to hardening parameter
384 hardp(i) = sqrt(alpha1(i)*dr1_dp(i)**2 + alpha2(i)*dr2_dp(i)**2 + alpha3(i)*dr1c_dp(i)**2
385 . + alpha4(i)*dr2c_dp(i)**2 + two*alpha5(i)*drt_dp(i)**2)
386 dphi_dpla = dphi_dr1*dr1_dp(i)*alpha1(i) + dphi_dr2*dr2_dp(i)*alpha2(i) +
387 . dphi_dr1c*dr1c_dp(i)*alpha3(i) + dphi_dr2c*dr2c_dp(i)*alpha4(i) +
388 . dphi_drt*drt_dp(i)*alpha5(i)
389c
390 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
391 dphi_dlam = -dfdsig2 + dphi_dpla
392 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
393c
394 ! 4 - Computation of plastic multiplier and variables update
395 !----------------------------------------------------------
396c
397 ! a) Restoring previous value of the yield function
398 phi(i) = phi0(i)
399c
400 ! b) Computation of the trial stress increment
401 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
402 dsigyy(i) = a21*depsxx(i) + a22*depsyy(i)
403 dsigxy(i) = depsxy(i)*g12
404c
405 ! c) Computation of yield surface trial increment DPHI
406 dphi = normxx * dsigxx(i)
407 . + normyy * dsigyy(i)
408 . + two * normxy * dsigxy(i)
409c
410 ! d) Computation of the plastic multiplier
411 dlam = -(phi(i) + dphi) / dphi_dlam
412 dlam = max(dlam, zero)
413c
414 ! e) Plastic strains tensor update
415 dpxx = dlam * normpxx
416 dpyy = dlam * normpyy
417 dpxy = two * dlam * normpxy
418c
419 ! f) Elasto-plastic stresses update
420 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
421 signyy(i) = sigoyy(i) + dsigyy(i) - (a21*dpxx + a22*dpyy)
422 signxy(i) = sigoxy(i) + dsigxy(i) - g12*dpxy
423c
424 ! g) Cumulated plastic strain and hardening parameter update
425 ddep = dlam
426 dpla(i) = max(zero, dpla(i) + ddep)
427 dpla2(i) = max(zero, dpla2(i) + alpha1(i)*ddep)
428 dpla3(i) = max(zero, dpla3(i) + alpha2(i)*ddep)
429 dpla4(i) = max(zero, dpla4(i) + alpha3(i)*ddep)
430 dpla5(i) = max(zero, dpla5(i) + alpha4(i)*ddep)
431 dpla6(i) = max(zero, dpla6(i) + alpha5(i)*ddep)
432 pla(i,1) = pla(i,1) + ddep
433 pla(i,2) = pla(i,2) + alpha1(i)*ddep
434 pla(i,3) = pla(i,3) + alpha2(i)*ddep
435 pla(i,4) = pla(i,4) + alpha3(i)*ddep
436 pla(i,5) = pla(i,5) + alpha4(i)*ddep
437 pla(i,6) = pla(i,6) + alpha5(i)*ddep
438c
439 ! h) Update of Alpha coefficient
440 alpha1(i) = zero
441 alpha2(i) = zero
442 alpha3(i) = zero
443 alpha4(i) = zero
444 alpha5(i) = zero
445 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
446 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
447 IF (-signxx(i) > zero) alpha3(i) = one
448 IF (-signyy(i) > zero) alpha4(i) = one
449 IF (abs(signxy(i)) > zero) alpha5(i) = one
450c
451 ! i) Yield stresses update
452 IF (itab == 0) THEN
453 ! i) Tensile yield stresses update
454 ! - in direction 1 (MD)
455 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
456 ! - in direction 2 (CD)
457 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
458 ! ii) Compressive yield stresses update
459 ! - in direction 1 (MD)
460 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
461 r1c(i) = r1c(i)/sqrt(one - g1c)
462 ! - in direction 2 (CD)
463 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
464 ! iii) Shear yield stress update
465 eta1(i) = one + alpha3(i)*alpha4(i)*d1
466 eta2(i) = one + alpha3(i)*alpha4(i)*d2
467 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
468c
469 ! i) Yield function value update
470 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
471 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
472 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
473 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
474 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
475 ENDIF
476c
477 ENDDO
478 !===================================================================
479 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
480 !===================================================================
481c
482 ! If tabulated yield function, update of the yield stress for all element
483 IF (itab > 0) THEN
484 IF (nindx > 0) THEN
485 ! -> Tensile yield stress in direction 1 (MD)
486 xvec(1:nel,1) = pla(1:nel,2)
487 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
488 ipos(1:nel,1) = vartmp(1:nel,1)
489 ipos(1:nel,2) = 1
490 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
491 r1(1:nel) = r1(1:nel) * yscale1
492 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
493 vartmp(1:nel,1) = ipos(1:nel,1)
494 ! -> Compressive yield stress in direction 1 (MD)
495 xvec(1:nel,1) = pla(1:nel,3)
496 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
497 ipos(1:nel,1) = vartmp(1:nel,2)
498 ipos(1:nel,2) = 1
499 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
500 r2(1:nel) = r2(1:nel) * yscale2
501 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
502 vartmp(1:nel,2) = ipos(1:nel,1)
503 ! -> Positive shear yield stress
504 xvec(1:nel,1) = pla(1:nel,4)
505 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
506 ipos(1:nel,1) = vartmp(1:nel,3)
507 ipos(1:nel,2) = 1
508 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
509 r1c(1:nel) = r1c(1:nel) * yscale1c
510 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
511 vartmp(1:nel,3) = ipos(1:nel,1)
512 ! -> Tensile yield stress in direction 2 (CD)
513 xvec(1:nel,1) = pla(1:nel,5)
514 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
515 ipos(1:nel,1) = vartmp(1:nel,4)
516 ipos(1:nel,2) = 1
517 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
518 r2c(1:nel) = r2c(1:nel) * yscale2c
519 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
520 vartmp(1:nel,4) = ipos(1:nel,1)
521 ! Compressive yield stress in direction 2 (CD)
522 xvec(1:nel,1) = pla(1:nel,6)
523 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
524 ipos(1:nel,1) = vartmp(1:nel,5)
525 ipos(1:nel,2) = 1
526 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
527 rt(1:nel) = rt(1:nel) * yscalet
528 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
529 vartmp(1:nel,5) = ipos(1:nel,1)
530 ! Yield function value update
531 DO i = 1, nel
532 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
533 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
534 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
535 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
536 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
537 ENDDO
538 ENDIF
539 ENDIF
540c
541 ! Storing new values
542 DO i=1,nel
543 ! USR Outputs
544 uvar(i,1) = phi(i) ! Yield function value
545 ! Plastic strain-rate
546 IF (itab == 0) THEN
547 epsd(i,1) = dpla(i) / max(em20,timestep)
548 epsd(i,2) = dpla2(i) / max(em20,timestep)
549 epsd(i,3) = dpla3(i) / max(em20,timestep)
550 epsd(i,4) = dpla4(i) / max(em20,timestep)
551 epsd(i,5) = dpla5(i) / max(em20,timestep)
552 epsd(i,6) = dpla6(i) / max(em20,timestep)
553 ELSE
554 dpdt = dpla(i)/max(em20,timestep)
555 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
556 dpdt = dpla2(i)/max(em20,timestep)
557 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
558 dpdt = dpla3(i)/max(em20,timestep)
559 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
560 dpdt = dpla4(i)/max(em20,timestep)
561 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
562 dpdt = dpla5(i)/max(em20,timestep)
563 epsd(i,5) = asrate * dpdt + (one - asrate) * epsd(i,5)
564 dpdt = dpla6(i)/max(em20,timestep)
565 epsd(i,6) = asrate * dpdt + (one - asrate) * epsd(i,6)
566 ENDIF
567 ! Coefficient for hourglass
568 et(i) = hardp(i) / (hardp(i) + max(young1,young2))
569 ! Computation of the sound speed
570 soundsp(i) = sqrt(max(a11,a12,a21,a22,g12,g23,g31)/ rho(i))
571 ! Storing the yield stress
572 sigy(i) = sqrt(r1(i)**2 + r2(i)**2 + r1c(i)**2 +
573 . r2c(i)**2 + two*rt(i)**2)
574 ENDDO
575c
576 END
alpha2
#define alpha2
Definition eval.h:48
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
max
#define max(a, b)
Definition macros.h:21
mat107c_nice
subroutine mat107c_nice(nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho, pla, dpla, epsd, soundsp, shf, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigy, et, numtabl, itable, table, nvartmp, vartmp)
Definition mat107c_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