OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat122c_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 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)

Function/Subroutine Documentation

◆ mat122c_nice()

subroutine mat122c_nice ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
dimension(nuparam), intent(in) uparam,
intent(inout) uvar,
intent(in) epsxx,
intent(in) epsyy,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(inout) thk,
intent(in) thkly,
intent(inout) off,
intent(out) sigy,
intent(out) etse,
intent(inout) dmg,
intent(inout) seq,
intent(in) epspl,
intent(in) shf,
intent(out) 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 mat122c_nice.F.

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
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