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