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