OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat122_nice.F File Reference
#include "implicit_f.inc"
#include "tabsiz_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ mat122_nice()

subroutine mat122_nice ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
dimension(nuparam), intent(in) uparam,
intent(inout) uvar,
intent(in) rho0,
intent(in) epsxx,
intent(in) epsyy,
intent(in) epszz,
intent(inout) pla,
intent(inout) dpla,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(inout) signxx,
intent(inout) signyy,
intent(inout) signzz,
intent(inout) signxy,
intent(inout) signyz,
intent(inout) signzx,
intent(inout) off,
intent(inout) sigy,
intent(inout) et,
intent(inout) dmg,
intent(inout) seq,
intent(inout) epsd,
intent(inout) soundsp,
integer, intent(in) nfunc,
integer, dimension(nfunc), intent(in) ifunc,
integer, dimension(snpc), intent(in) npf,
dimension(stf), intent(in) tf,
integer, intent(in) nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp )

Definition at line 30 of file mat122_nice.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72