OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat122_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!|| mat122_nice ../engine/source/materials/mat/mat122/mat122_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps122 ../engine/source/materials/mat/mat122/sigeps122.F
27!||--- calls -----------------------------------------------------
28!|| vinter ../engine/source/tools/curve/vinter.F
29!||====================================================================
30 SUBROUTINE mat122_nice(
31 1 NEL ,NUPARAM ,NUVAR ,UPARAM ,UVAR ,RHO0 ,
32 2 EPSXX ,EPSYY ,EPSZZ ,PLA ,DPLA ,
33 3 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
34 4 SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
35 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
36 6 OFF ,SIGY ,ET ,DMG ,SEQ ,EPSD ,
37 7 SOUNDSP ,NFUNC ,IFUNC ,NPF ,TF ,NVARTMP ,
38 8 VARTMP )
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C C O M M O N
45C-----------------------------------------------
46#include "tabsiz_c.inc"
47C-----------------------------------------------
48C I N P U T A r g u m e n t s
49C-----------------------------------------------
50 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,
51 . NFUNC,IFUNC(NFUNC),NPF(SNPC),NVARTMP
52 my_real, INTENT(IN) ::
53 . UPARAM(NUPARAM),TF(STF)
54 my_real,DIMENSION(NEL), INTENT(IN) ::
55 . RHO0,EPSXX,EPSYY,EPSZZ,
56 . DEPSXX,DEPSYY,DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,
57 . SIGOYY,SIGOZZ,SIGOXY,SIGOYZ,SIGOZX
58 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
59 . signxx,signyy,signzz,signxy,signyz,signzx,
60 . soundsp,sigy,et
61 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
62 . pla,dpla,epsd,off,seq
63 my_real ,DIMENSION(NEL,6), INTENT(INOUT) ::
64 . dmg
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67 INTEGER, DIMENSION(NEL,NVARTMP), INTENT(INOUT) ::
68 . VARTMP
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER
73 . i,ii,j,niter,iter,nindx,index(nel),
74 . ish,itr,ires,ibuck,icomp,ltype11,ltype12,
75 . ltyper0,ipos(nel),ilen(nel),iad(nel)
76 my_real
77 . e10,e20,e30,nu12,nu21,nu23,nu32,nu13,nu31,
78 . g120,g230,g310,e1c,gamma,sigy0,beta,m,a,efti0,
79 . eftu0,dftu,efci0,efcu0,dfcu,dsat1,y00,yc0,b,
80 . dmax,yr,ysp,dsat2,y0p0,ycp0,dsat2c,y0pc0,ycpc0,
81 . epsd11,d11,n11,d11u,n11u,epsd12,d22,n22,d12,
82 . n12,epsdr0,dr0,nr0
83 my_real
84 . ddep,dfdsig2,dlam,normyy,normzz,normxy,normyz,
85 . normzx,sig_dfdsig,h(nel),dpyy(nel),dpzz(nel),
86 . dpxy(nel),dpyz(nel),dpzx(nel),phi(nel),dpla_dlam(nel),
87 . dphi_dlam(nel),dydx(nel),s11(nel),s12(nel),s13(nel),
88 . s22(nel),s23(nel),s33(nel),c11,c12,c13,c22,c23,c33,
89 . detc,dft(nel),dfc(nel),e1(nel),e2(nel),epsf_eq,zd(nel),
90 . zdp(nel),y(nel),yp(nel),d(nel),dp(nel),df(nel),
91 . g12(nel),efti(nel),eftu(nel),efci(nel),efcu(nel),
92 . f11(nel),f22(nel),f12(nel),f11r(nel),fr0(nel),e3(nel),
93 . y0(nel),yc(nel),y0p(nel),ycp(nel),y0pc(nel),ycpc(nel),
94 . dpy(nel),dpz(nel),g23(nel),g31(nel),seq0(nel),phi0(nel),
95 . dsigyy(nel),dsigzz(nel),dsigxy(nel),dsigyz(nel),dsigzx(nel),
96 . dphi,var(nel),epspyy(nel),epspzz(nel)
97C-----------------------------------------------
98C USER VARIABLES INITIALIZATION
99C-----------------------------------------------
100C
101 ! Elastic parameters
102 e10 = uparam(1)
103 e20 = uparam(2)
104 e30 = uparam(3)
105 nu12 = uparam(4)
106 nu21 = uparam(5)
107 nu13 = uparam(6)
108 nu31 = uparam(7)
109 nu23 = uparam(8)
110 nu32 = uparam(9)
111 g120 = uparam(10)
112 g230 = uparam(11)
113 g310 = uparam(12)
114 e1c = uparam(13)
115 gamma = uparam(14)
116 ish = nint(uparam(15))
117 itr = nint(uparam(16))
118 ires = nint(uparam(17))
119 sigy0 = uparam(18)
120 beta = uparam(19)
121 m = uparam(20)
122 a = uparam(21)
123 efti0 = uparam(22)
124 eftu0 = uparam(23)
125 dftu = uparam(24)
126 efci0 = uparam(25)
127 efcu0 = uparam(26)
128 dfcu = uparam(27)
129 ibuck = nint(uparam(28))
130 dsat1 = uparam(29)
131 y00 = uparam(30)
132 yc0 = uparam(31)
133 b = uparam(32)
134 dmax = uparam(33)
135 yr = uparam(34)
136 ysp = uparam(35)
137 dsat2 = uparam(36)
138 y0p0 = uparam(37)
139 ycp0 = uparam(38)
140 dsat2c = uparam(39)
141 y0pc0 = uparam(40)
142 ycpc0 = uparam(41)
143 epsd11 = uparam(42)
144 d11 = uparam(43)
145 n11 = uparam(44)
146 d11u = uparam(45)
147 n11u = uparam(46)
148 epsd12 = uparam(47)
149 d22 = uparam(48)
150 n22 = uparam(49)
151 d12 = uparam(50)
152 n12 = uparam(51)
153 epsdr0 = uparam(52)
154 dr0 = uparam(53)
155 nr0 = uparam(54)
156 ltype11 = nint(uparam(55))
157 ltype12 = nint(uparam(56))
158 ltyper0 = nint(uparam(57))
159C
160 ! Recovering internal variables
161 DO i=1,nel
162 ! Checking deletion flag value
163 IF (off(i) < one) off(i) = four_over_5*off(i)
164 IF (off(i) < em01) off(i) = zero
165 ! Hourglass coefficient
166 et(i) = one
167 h(i) = zero
168 ! Old yield stress
169 seq0(i) = seq(i)
170 ! Previous yield stress
171 phi0(i) = uvar(i,13)
172 ! Plastic strain increment
173 dpla(i) = zero
174 dpyy(i) = zero
175 dpzz(i) = zero
176 dpxy(i) = zero
177 dpyz(i) = zero
178 dpzx(i) = zero
179 ! Damage variables and user variables
180 df(i) = dmg(i,2)
181 d(i) = dmg(i,3)
182 dp(i) = dmg(i,4)
183 dft(i) = dmg(i,5)
184 dfc(i) = dmg(i,6)
185 y(i) = uvar(i,1)
186 yp(i) = uvar(i,2)
187 dpy(i) = uvar(i,14)
188 dpz(i) = uvar(i,15)
189 epspyy(i) = uvar(i,16)
190 epspzz(i) = uvar(i,17)
191 ENDDO
192C
193 ! Compute strain rate dependency factor
194 DO i = 1,nel
195 ! Rate dependency in fiber direction 1 for Young modulus
196 IF (ltype11 == 1) THEN
197 f11(i) = d11*(abs(epsd(i))/epsd11)**n11
198 ELSEIF (ltype11 == 2) THEN
199 f11(i) = d11*(abs(epsd(i))/epsd11) + n11
200 ELSEIF (ltype11 == 3) THEN
201 f11(i) = d11*log(max(abs(epsd(i))/epsd11,one)) + log(n11)
202 ELSEIF (ltype11 == 4) THEN
203 f11(i) = d11*tanh(n11*(max(abs(epsd(i))-epsd11,zero)))
204 ENDIF
205 ! Rate dependency in matrix transverse direction 2 for Young modulus
206 IF (ltype12 == 1) THEN
207 f22(i) = d22*(abs(epsd(i))/epsd12)**n22
208 ELSEIF (ltype12 == 2) THEN
209 f22(i) = d22*(abs(epsd(i))/epsd12) + n22
210 ELSEIF (ltype12 == 3) THEN
211 f22(i) = d22*log(max(abs(epsd(i))/epsd12,one)) + log(n22)
212 ELSEIF (ltype12 == 4) THEN
213 f22(i) = d22*tanh(n22*(max(abs(epsd(i))-epsd12,zero)))
214 ENDIF
215 ! Rate dependency for shear plane 12 modulus
216 IF (ltype12 == 1) THEN
217 f12(i) = d12*(abs(epsd(i))/epsd12)**n12
218 ELSEIF (ltype12 == 2) THEN
219 f12(i) = d12*(abs(epsd(i))/epsd12) + n12
220 ELSEIF (ltype12 == 3) THEN
221 f12(i) = d12*log(max(abs(epsd(i))/epsd12,one)) + log(n12)
222 ELSEIF (ltype12 == 4) THEN
223 f12(i) = d12*tanh(n12*(max(abs(epsd(i))-epsd12,zero)))
224 ENDIF
225 ! Rate dependency in fiber direction 1 for Young modulus
226 IF (ltype11 == 1) THEN
227 f11r(i) = d11u*(abs(epsd(i))/epsd11)**n11u
228 ELSEIF (ltype11 == 2) THEN
229 f11r(i) = d11u*(abs(epsd(i))/epsd11) + n11u
230 ELSEIF (ltype11 == 3) THEN
231 f11r(i) = d11u*log(max(abs(epsd(i))/epsd11,one)) + log(n11u)
232 ELSEIF (ltype11 == 4) THEN
233 f11r(i) = d11u*tanh(n11u*(max(abs(epsd(i))-epsd11,zero)))
234 ENDIF
235 ! Rate dependency in fiber direction 1 for Young modulus
236 IF (ltyper0 == 1) THEN
237 fr0(i) = dr0*(abs(epsd(i))/epsdr0)**nr0
238 ELSEIF (ltyper0 == 2) THEN
239 fr0(i) = dr0*(abs(epsd(i))/epsdr0) + nr0
240 ELSEIF (ltyper0 == 3) THEN
241 fr0(i) = dr0*log(max(abs(epsd(i))/epsdr0,one)) + log(nr0)
242 ELSEIF (ltyper0 == 4) THEN
243 fr0(i) = dr0*tanh(nr0*(max(abs(epsd(i))-epsdr0,zero)))
244 ENDIF
245 ENDDO
246c
247 ! Elastic parameters, yield stress and strain rate dependency
248 DO i = 1,nel
249 ! Fiber (direction 1)
250 ! -> Tension
251 IF (epsxx(i) >= zero) THEN
252 e1(i) = e10
253 ! -> Compression
254 ELSE
255 e1(i) = e1c/(one + (gamma*e1c*abs(epsxx(i))))
256 ENDIF
257 e1(i) = e1(i)*(one + f11(i))
258 ! Matrix (direction 2)
259 e2(i) = e20*(one + f22(i))
260 ! Matrix (direction 3)
261 e3(i) = e30*(one + f22(i))
262 ! Shear moduli
263 g12(i) = g120*(one + f12(i))
264 g23(i) = g230*(one + f12(i))
265 g31(i) = g310*(one + f12(i))
266 ! Compliance matrix for 3D
267 c11 = one/e1(i)
268 c22 = one/e2(i)
269 c33 = one/e3(i)
270 c12 = -nu12/e1(i)
271 c13 = -nu31/e3(i)
272 c23 = -nu23/e2(i)
273 detc = c11*c22*c33-c11*c23*c23-c12*c12*c33+c12*c13*c23
274 . + c13*c12*c23-c13*c22*c13
275 ! Stiffness matrix for 3D
276 s11(i) = (c22*c33-c23*c23)/detc
277 s12(i) = -(c12*c33-c13*c23)/detc
278 s13(i) = (c12*c23-c13*c22)/detc
279 s22(i) = (c11*c33-c13*c13)/detc
280 s23(i) = -(c11*c23-c13*c12)/detc
281 s33(i) = (c11*c22-c12*c12)/detc
282 ! Yield stress
283 sigy(i) = sigy0*(one + fr0(i)) + beta*exp(m*log(pla(i)+em20))
284 ! Rate dependency on fiber failure
285 efti(i) = uvar(i,3)
286 IF (uvar(i,3) == zero) efti(i) = efti0*(one + f11r(i))
287 eftu(i) = uvar(i,4)
288 IF (uvar(i,4) == zero) eftu(i) = eftu0*(one + f11r(i))
289 efci(i) = uvar(i,5)
290 IF (uvar(i,5) == zero) efci(i) = efci0*(one + f11r(i))
291 efcu(i) = uvar(i,6)
292 IF (uvar(i,6) == zero) efcu(i) = efcu0*(one + f11r(i))
293 ! Rate dependency on matrix failure
294 y0(i) = uvar(i,7)
295 IF (uvar(i,7) == zero) y0(i) = y00*sqrt(one + f12(i))
296 yc(i) = uvar(i,8)
297 IF (uvar(i,8) == zero) yc(i) = yc0*sqrt(one + f12(i))
298 y0p(i) = uvar(i,9)
299 IF (uvar(i,9) == zero) y0p(i) = y0p0*sqrt(one + f22(i))
300 ycp(i) = uvar(i,10)
301 IF (uvar(i,10) == zero) ycp(i) = ycp0*sqrt(one + f22(i))
302 y0pc(i) = uvar(i,11)
303 IF (uvar(i,11) == zero) y0pc(i) = y0pc0*sqrt(one + f22(i))
304 ycpc(i) = uvar(i,12)
305 IF (uvar(i,12) == zero) ycpc(i) = ycpc0*sqrt(one + f22(i))
306 ENDDO
307c
308 !========================================================================
309 ! - COMPUTATION OF TRIAL VALUES
310 !========================================================================
311 DO i=1,nel
312c
313 ! computation of the trial stress tensor
314 signyy(i) = sigoyy(i)/max((one-dpy(i)),em20)
315 . + s12(i)*depsxx(i) + s22(i)*depsyy(i) + s23(i)*depszz(i)
316 signzz(i) = sigozz(i)/max((one-dpz(i)),em20)
317 . + s13(i)*depsxx(i) + s23(i)*depsyy(i) + s33(i)*depszz(i)
318 signxy(i) = sigoxy(i)/max((one- d(i)),em20) + g12(i)*depsxy(i)
319 signyz(i) = sigoyz(i)/max((one- d(i)),em20) + g23(i)*depsyz(i)
320 signzx(i) = sigozx(i)/max((one- d(i)),em20) + g31(i)*depszx(i)
321C
322 ! Equivalent stress
323 seq(i) = sqrt((signxy(i))**2 + (signyz(i))**2 + (signzx(i))**2 +
324 . a*(signyy(i))**2 + a*(signzz(i))**2)
325C
326 ENDDO
327c
328 !========================================================================
329 ! - COMPUTATION OF YIELD FONCTION
330 !========================================================================
331 phi(1:nel) = seq(1:nel) - sigy(1:nel)
332 ! Checking plastic behavior for all elements
333 nindx = 0
334 index(1:nel) = 0
335 DO i=1,nel
336 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
337 nindx = nindx+1
338 index(nindx) = i
339 ENDIF
340 ENDDO
341c
342 !====================================================================================
343 ! - PLASTIC CORRECTION WITH N.I.C.E EXPLICIT ALGORITHM (NEXT INCREMENT CORRECT ERROR)
344 !====================================================================================
345c
346#include "vectorize.inc"
347 ! Loop over yielding elements
348 DO ii=1,nindx
349 i = index(ii)
350c
351 ! Computation of the trial stress increment
352 dsigyy(i) = s12(i)*depsxx(i) + s22(i)*depsyy(i) + s23(i)*depszz(i)
353 dsigzz(i) = s13(i)*depsxx(i) + s23(i)*depsyy(i) + s33(i)*depszz(i)
354 dsigxy(i) = g12(i)*depsxy(i)
355 dsigyz(i) = g23(i)*depsyz(i)
356 dsigzx(i) = g31(i)*depszx(i)
357c
358 ! Note : in this part, the purpose is to compute for each iteration
359 ! a plastic multiplier allowing to update internal variables to satisfy
360 ! the consistency condition.
361 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
362 ! -> PHI : old value of yield function (known)
363 ! -> DPHI : yield function prediction (to compute)
364 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
365 ! into account of internal variables kinetic :
366 ! plasticity, temperature and damage (to compute)
367c
368 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
369 !-------------------------------------------------------------
370 normyy = a*(sigoyy(i)/max((one-dpy(i)),em20))/max(seq0(i),em20)
371 normzz = a*(sigozz(i)/max((one-dpz(i)),em20))/max(seq0(i),em20)
372 normxy = (sigoxy(i)/max((one - d(i)),em20))/max(seq0(i),em20)
373 normyz = (sigoyz(i)/max((one - d(i)),em20))/max(seq0(i),em20)
374 normzx = (sigozx(i)/max((one - d(i)),em20))/max(seq0(i),em20)
375c
376 ! Restoring previous value of the yield function
377 phi(i) = phi0(i)
378c
379 ! Computation of yield surface trial increment DPHI
380 dphi = normyy*dsigyy(i) + normzz*dsigzz(i) + normxy*dsigxy(i) +
381 . normyz*dsigyz(i) + normzx*dsigzx(i)
382c
383 ! 2 - Computation of DPHI_DLAMBDA
384 !---------------------------------------------------------
385c
386 ! a) Derivative with respect stress increments tensor DSIG
387 ! --------------------------------------------------------
388 dfdsig2 = normyy*(s22(i)*normyy + s23(i)*normzz)
389 . + normzz*(s23(i)*normyy + s33(i)*normzz)
390 . + normxy*normxy*g12(i)
391 . + normyz*normyz*g23(i)
392 . + normzx*normzx*g31(i)
393c
394 ! b) Derivatives with respect to plastic strain P
395 ! ------------------------------------------------
396c
397 ! i) Derivative of the yield stress with respect to plastic strain dSIGY / dPLA
398 ! ----------------------------------------------------------------------------
399 h(i) = beta*m*exp((m-1)*log(pla(i)+em20))
400 h(i) = min(h(i),max(two*g120,e20))
401c
402 ! ii) Derivative of dPLA with respect to DLAM
403 ! -------------------------------------------
404 sig_dfdsig = (sigoyy(i)/max((one-dpy(i)),em20)) * normyy
405 . + (sigozz(i)/max((one-dpz(i)),em20)) * normzz
406 . + (sigoxy(i)/max((one - d(i)),em20)) * normxy
407 . + (sigoyz(i)/max((one - d(i)),em20)) * normyz
408 . + (sigozx(i)/max((one - d(i)),em20)) * normzx
409 dpla_dlam(i) = sig_dfdsig/max(sigy(i),em20)
410c
411 ! 3 - Computation of plastic multiplier and variables update
412 !----------------------------------------------------------
413c
414 ! Derivative of PHI with respect to DLAM
415 dphi_dlam(i) = - dfdsig2 - h(i)*dpla_dlam(i)
416 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
417c
418 ! Computation of the plastic multiplier
419 dlam = -(phi(i)+dphi)/dphi_dlam(i)
420c
421 ! Plastic strains tensor increment
422 dpyy(i) = dlam*normyy
423 dpzz(i) = dlam*normzz
424 dpxy(i) = dlam*normxy
425 dpyz(i) = dlam*normyz
426 dpzx(i) = dlam*normzx
427c
428 ! Total plastic strains along Y and Z
429 epspyy(i) = epspyy(i) + dpyy(i)
430 epspzz(i) = epspzz(i) + dpzz(i)
431c
432 ! Elasto-plastic stresses update
433 signxx(i) = signxx(i) - (s12(i)*dpyy(i) + s13(i)*dpzz(i))
434 signyy(i) = signyy(i) - (s22(i)*dpyy(i) + s23(i)*dpzz(i))
435 signzz(i) = signzz(i) - (s23(i)*dpyy(i) + s33(i)*dpzz(i))
436 signxy(i) = signxy(i) - dpxy(i)*g12(i)
437 signyz(i) = signyz(i) - dpyz(i)*g23(i)
438 signzx(i) = signzx(i) - dpzx(i)*g31(i)
439c
440 ! Cumulated plastic strain and strain rate update
441 ddep = dlam*dpla_dlam(i)
442 dpla(i) = max(zero, dpla(i) + ddep)
443 pla(i) = max(zero, pla(i) + ddep)
444c
445 ! Update equivalent stress
446 seq(i) = sqrt((signxy(i))**2 + (signyz(i))**2 + (signzx(i))**2 +
447 . a*(signyy(i))**2 + a*(signzz(i))**2)
448c
449 ! Update the hardening yield stress
450 sigy(i) = sigy(i) + h(i)*dlam*dpla_dlam(i)
451c
452 ! Update yield function value
453 phi(i) = seq(i) - sigy(i)
454c
455 ENDDO
456 ! End of the loop over the yielding elements
457c
458 !===================================================================
459 ! - END OF PLASTIC CORRECTION WITH N.I.C.E EXPLICIT ALGORITHM
460 !===================================================================
461c
462 !===================================================================
463 ! - DAMAGE VARIABLES COMPUTATION AND UPDATE STRESS TENSOR
464 !===================================================================
465 DO i = 1,nel
466c
467 ! Fiber damage (direction 1)
468 ! ------------------------------------------
469 ! -> Fiber equivalent strain
470 epsf_eq = (one - nu23*nu32)*epsxx(i) +
471 . (nu23*nu31 + nu21)*(epsyy(i)-epspyy(i)) +
472 . (nu21*nu32 + nu31)*(epszz(i)-epspzz(i))
473 ! -> Tension
474 IF (epsf_eq >= zero) THEN
475 IF (epsf_eq < efti(i)) THEN
476 dft(i) = max(zero,dft(i))
477 ELSEIF ((epsf_eq >= efti(i)).AND.(epsf_eq < eftu(i))) THEN
478 dft(i) = max(dftu*((epsf_eq-efti(i))/(eftu(i)-efti(i))),dft(i))
479 ! Save damage threshold in case of strain rate dependency
480 IF (uvar(i,3) == zero) uvar(i,3) = efti(i)
481 IF (uvar(i,4) == zero) uvar(i,4) = eftu(i)
482 ELSEIF (epsf_eq >= eftu(i)) THEN
483 dft(i) = max(one - (one - dftu)*(eftu(i)/epsf_eq),dft(i))
484 ENDIF
485 dft(i) = max(dft(i),zero)
486 dft(i) = min(dft(i),one)
487 df(i) = dft(i)
488 ! -> Compression
489 ELSEIF ((epsf_eq < zero).AND.(ibuck > 1)) THEN
490 IF (abs(epsf_eq) < efci(i)) THEN
491 dfc(i) = max(zero,dfc(i))
492 ELSEIF ((abs(epsf_eq) >= efci(i)).AND.(abs(epsf_eq) < efcu(i))) THEN
493 dfc(i) = max(dfcu*((abs(epsf_eq)-efci(i))/(efcu(i)-efci(i))),dfc(i))
494 ! Save damage threshold in case of strain rate dependency
495 IF (uvar(i,5) == zero) uvar(i,5) = efci(i)
496 IF (uvar(i,6) == zero) uvar(i,6) = efcu(i)
497 ELSEIF (abs(epsf_eq) >= efcu(i)) THEN
498 dfc(i) = max(one - (one - dfcu)*(efcu(i)/abs(epsf_eq)),dfc(i))
499 ENDIF
500 dfc(i) = max(dfc(i),zero)
501 dfc(i) = min(dfc(i),one)
502 df(i) = dfc(i)
503 ENDIF
504c
505 ! Matrix damage (direction 2 and 3)
506 ! ------------------------------------------
507 ! -> Damage functions (derivatives of elastic energy)
508 zd(i) = half*((signxy(i)**2/g120) + (signyz(i)**2/g230) + (signzx(i)**2/g310))
509 zdp(i) = half*((((max(signyy(i),zero))**2)/e20) + (((max(signzz(i),zero))**2)/e30))
510 ! -> Damage evolution functions
511 y(i) = max(y(i) ,sqrt(zd(i)+b*zdp(i)))
512 yp(i) = max(yp(i),sqrt(zdp(i)))
513 ENDDO
514c
515 ! -> Shear damage
516 ! Linear
517 IF (ish == 1) THEN
518 DO i = 1,nel
519 IF (y(i)<y0(i)) THEN
520 d(i) = zero
521 ELSEIF ((d(i)<dmax).AND.(y(i)<ysp).AND.(y(i)<yr)) THEN
522 d(i) = max(y(i)-y0(i),zero)/max(yc(i),em20)
523 d(i) = min(d(i),dmax)
524 ! Save damage threshold in case of strain rate dependency
525 IF (uvar(i,7) == zero) uvar(i,7) = y0(i)
526 IF (uvar(i,8) == zero) uvar(i,8) = yc(i)
527 ELSE
528 d(i) = one - (one - dmax)*uvar(i,1)/max(y(i),em20)
529 ENDIF
530 d(i) = max(d(i),zero)
531 d(i) = min(d(i), one)
532 ENDDO
533 ! Exponential
534 ELSEIF (ish == 2) THEN
535 DO i = 1,nel
536 IF (y(i)>y0(i)) THEN
537 d(i) = dsat1*(one - exp((y0(i)-y(i))/max(yc(i),em20)))
538 ! Save damage threshold in case of strain rate dependency
539 IF (uvar(i,7) == zero) uvar(i,7) = y0(i)
540 IF (uvar(i,8) == zero) uvar(i,8) = yc(i)
541 ELSE
542 d(i) = zero
543 ENDIF
544 d(i) = max(d(i),zero)
545 d(i) = min(d(i), one)
546 ENDDO
547 ! Tabulated function
548 ELSEIF (ish == 3) THEN
549 ipos(1:nel) = vartmp(1:nel,1)
550 iad(1:nel) = npf(ifunc(1)) / 2 + 1
551 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
552 DO i = 1,nel
553 var(i) = y(i)/y0(i)
554 ENDDO
555 CALL vinter(tf,iad,ipos,ilen,nel,var,dydx,d)
556 vartmp(1:nel,1) = ipos(1:nel)
557 DO i = 1,nel
558 ! Save damage threshold in case of strain rate dependency
559 IF (uvar(i,7) == zero .AND. d(i) /= zero) uvar(i,7) = y0(i)
560 d(i) = max(d(i),zero)
561 d(i) = min(d(i), one)
562 ENDDO
563 ENDIF
564c
565 ! -> Transverse damage
566 ! Linear
567 IF (itr == 1) THEN
568 DO i = 1,nel
569 IF (yp(i)<y0p(i)) THEN
570 dp(i) = zero
571 ELSEIF ((dp(i)<dmax).AND.(yp(i)<ysp).AND.(yp(i)<yr)) THEN
572 dp(i) = max(yp(i)-y0p(i),zero)/max(ycp(i),em20)
573 dp(i) = min(dp(i),dmax)
574 ! Save damage threshold in case of strain rate dependency
575 IF (uvar(i,9) == zero) uvar(i,9) = y0p(i)
576 IF (uvar(i,10) == zero) uvar(i,10) = ycp(i)
577 ELSE
578 dp(i) = one - (one - dmax)*uvar(i,2)/max(yp(i),em20)
579 ENDIF
580 dp(i) = max(dp(i),zero)
581 dp(i) = min(dp(i), one)
582 ENDDO
583 ! Exponential
584 ELSEIF (itr == 2) THEN
585 DO i = 1,nel
586 IF (yp(i)>y0p(i)) THEN
587 dp(i) = dsat2*(one - exp((y0p(i)-yp(i))/max(ycp(i),em20)))
588 ! Save damage threshold in case of strain rate dependency
589 IF (uvar(i,9) == zero) uvar(i,9) = y0p(i)
590 IF (uvar(i,10) == zero) uvar(i,10) = ycp(i)
591 ELSE
592 dp(i) = zero
593 ENDIF
594 dp(i) = max(dp(i),zero)
595 dp(i) = min(dp(i), one)
596 ENDDO
597 ! Tabulated function
598 ELSEIF (itr == 3) THEN
599 ipos(1:nel) = vartmp(1:nel,2)
600 iad(1:nel) = npf(ifunc(2)) / 2 + 1
601 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
602 DO i = 1,nel
603 var(i) = yp(i)/y0p(i)
604 ENDDO
605 CALL vinter(tf,iad,ipos,ilen,nel,var,dydx,dp)
606 vartmp(1:nel,2) = ipos(1:nel)
607 DO i = 1,nel
608 ! Save damage threshold in case of strain rate dependency
609 IF (uvar(i,9) == zero .AND. dp(i) /= zero) uvar(i,9) = y0p(i)
610 dp(i) = max(dp(i),zero)
611 dp(i) = min(dp(i), one)
612 ENDDO
613 ENDIF
614c
615 DO i = 1,nel
616c
617 ! Sound-speed
618 soundsp(i) = sqrt(max(s11(i),s22(i),s33(i),
619 . two*g12(i),two*g23(i),two*g31(i))/rho0(i))
620c
621 ! Computation of effective DP for each transverse directions
622 IF (epsyy(i) >= zero) THEN
623 dpy(i) = dp(i)
624 ELSE
625 dpy(i) = zero
626 ENDIF
627 IF (epszz(i) >= zero) THEN
628 dpz(i) = dp(i)
629 ELSE
630 dpz(i) = zero
631 ENDIF
632c
633 ! compute damaged stiffness matrix
634 s11(i) = s11(i)*(one - df(i))
635 s12(i) = s12(i)*(one - df(i))*(one - dpy(i))
636 s13(i) = s13(i)*(one - df(i))*(one - dpz(i))
637 s22(i) = s22(i)*(one - dpy(i))
638 s23(i) = s23(i)*(one - dpy(i))*(one - dpz(i))
639 s33(i) = s33(i)*(one - dpz(i))
640c
641 ! Update stresses with damage softening effect
642 ! --------------------------------------------
643 ! -> If non-linear compression young modulus is used for fibers
644 IF (gamma > zero .AND. epsxx(i) < zero) THEN
645 signxx(i) = -(one/gamma)*log(one + gamma*e1c*abs(epsxx(i)))*(one - df(i))*(one + f11(i))
646 ! -> If linear elasticity is used for fibers
647 ELSE
648 signxx(i) = s11(i)*epsxx(i)
649 ENDIF
650 signxx(i) = signxx(i) +
651 . s12(i)*(epsyy(i) - epspyy(i)) +
652 . s13(i)*(epszz(i) - epspzz(i))
653 signyy(i) = s12(i)*epsxx(i) +
654 . s22(i)*(epsyy(i) - epspyy(i)) +
655 . s23(i)*(epszz(i) - epspzz(i))
656 signzz(i) = s13(i)*epsxx(i) +
657 . s23(i)*(epsyy(i) - epspyy(i)) +
658 . s33(i)*(epszz(i) - epspzz(i))
659 signxy(i) = signxy(i)*(one - d(i))
660 signyz(i) = signyz(i)*(one - d(i))
661 signzx(i) = signzx(i)*(one - d(i))
662c
663 ! -> store internal variables and damages
664 ! ------------------------------------------
665 dmg(i,1) = max(df(i),d(i),dp(i))
666 dmg(i,2) = df(i)
667 dmg(i,3) = d(i)
668 dmg(i,4) = dp(i)
669 dmg(i,5) = dft(i)
670 dmg(i,6) = dfc(i)
671 uvar(i,1) = y(i)
672 uvar(i,2) = yp(i)
673 uvar(i,14) = dpy(i)
674 uvar(i,15) = dpz(i)
675 uvar(i,16) = epspyy(i)
676 uvar(i,17) = epspzz(i)
677c
678 ENDDO
679 !===================================================================
680c
681 ! Sound-speed and thickness update
682 DO i=1,nel
683 ! Save old yield stress
684 IF (dpla(i) > zero) THEN
685 uvar(i,13) = phi(i)
686 ELSE
687 uvar(i,13) = zero
688 ENDIF
689 ! Coefficient for hourglass
690 IF (dpla(i)>zero) THEN
691 et(i) = h(i) / (h(i) + max(e1(i),e2(i),e3(i)))
692 ELSE
693 et(i) = one
694 ENDIF
695 ENDDO
696C
697 END
698C
the
end diagonal values have been computed in the(sparse) matrix id.SOL
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat122_nice
subroutine mat122_nice(nel, nuparam, nuvar, uparam, uvar, rho0, epsxx, epsyy, epszz, pla, dpla, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, off, sigy, et, dmg, seq, epsd, soundsp, nfunc, ifunc, npf, tf, nvartmp, vartmp)
Definition mat122_nice.F:39
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72