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

subroutine mat122c_newton ( 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_newton.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 . var(nel),dpt(nel),dptc(nel),epspyy(nel)
96C-----------------------------------------------
97C USER VARIABLES INITIALIZATION
98C-----------------------------------------------
99C
100 ! Elastic parameters
101 e10 = uparam(1)
102 e20 = uparam(2)
103 e30 = uparam(3)
104 nu12 = uparam(4)
105 nu21 = uparam(5)
106 nu13 = uparam(6)
107 nu31 = uparam(7)
108 nu23 = uparam(8)
109 nu32 = uparam(9)
110 g120 = uparam(10)
111 g23 = uparam(11)
112 g31 = uparam(12)
113 e1c = uparam(13)
114 gamma = uparam(14)
115 ish = nint(uparam(15))
116 itr = nint(uparam(16))
117 ires = nint(uparam(17))
118 sigy0 = uparam(18)
119 beta = uparam(19)
120 m = uparam(20)
121 a = uparam(21)
122 efti0 = uparam(22)
123 eftu0 = uparam(23)
124 dftu = uparam(24)
125 efci0 = uparam(25)
126 efcu0 = uparam(26)
127 dfcu = uparam(27)
128 ibuck = nint(uparam(28))
129 dsat1 = uparam(29)
130 y00 = uparam(30)
131 yc0 = uparam(31)
132 b = uparam(32)
133 dmax = uparam(33)
134 yr = uparam(34)
135 ysp = uparam(35)
136 dsat2 = uparam(36)
137 y0p0 = uparam(37)
138 ycp0 = uparam(38)
139 dsat2c = uparam(39)
140 y0pc0 = uparam(40)
141 ycpc0 = uparam(41)
142 epsd11 = uparam(42)
143 d11 = uparam(43)
144 n11 = uparam(44)
145 d11u = uparam(45)
146 n11u = uparam(46)
147 epsd12 = uparam(47)
148 d22 = uparam(48)
149 n22 = uparam(49)
150 d12 = uparam(50)
151 n12 = uparam(51)
152 epsdr0 = uparam(52)
153 dr0 = uparam(53)
154 nr0 = uparam(54)
155 ltype11 = nint(uparam(55))
156 ltype12 = nint(uparam(56))
157 ltyper0 = nint(uparam(57))
158C
159 ! Compression matrix damage flag
160 icomp = 0
161 IF (y0pc0 > zero) icomp = 1
162C
163 ! Recovering internal variables
164 DO i=1,nel
165 ! Checking deletion flag value
166 IF (off(i) < one) off(i) = four_over_5*off(i)
167 IF (off(i) < em01) off(i) = zero
168 ! Hourglass coefficient
169 etse(i) = one
170 h(i) = zero
171 ! Plastic strain increment
172 dpla(i) = zero
173 dpyy(i) = zero
174 dpxy(i) = zero
175 dpzz(i) = zero
176 ! Damage variables and user variables
177 df(i) = dmg(i,2)
178 d(i) = dmg(i,3)
179 dp(i) = dmg(i,4)
180 dft(i) = dmg(i,5)
181 dfc(i) = dmg(i,6)
182 y(i) = uvar(i,1)
183 yp(i) = uvar(i,2)
184 epspyy(i) = uvar(i,16)
185 ENDDO
186C
187 ! Compute strain rate dependency factor
188 DO i = 1,nel
189 ! Rate dependency in fiber direction 1 for Young modulus
190 IF (ltype11 == 1) THEN
191 f11(i) = d11*(abs(epspl(i))/epsd11)**n11
192 ELSEIF (ltype11 == 2) THEN
193 f11(i) = d11*(abs(epspl(i))/epsd11) + n11
194 ELSEIF (ltype11 == 3) THEN
195 f11(i) = d11*log(max(abs(epspl(i))/epsd11,one)) + log(n11)
196 ELSEIF (ltype11 == 4) THEN
197 f11(i) = d11*tanh(n11*(max(abs(epspl(i))-epsd11,zero)))
198 ENDIF
199 ! Rate dependency in matrix transverse direction 2 for Young modulus
200 IF (ltype12 == 1) THEN
201 f22(i) = d22*(abs(epspl(i))/epsd12)**n22
202 ELSEIF (ltype12 == 2) THEN
203 f22(i) = d22*(abs(epspl(i))/epsd12) + n22
204 ELSEIF (ltype12 == 3) THEN
205 f22(i) = d22*log(max(abs(epspl(i))/epsd12,one)) + log(n22)
206 ELSEIF (ltype12 == 4) THEN
207 f22(i) = d22*tanh(n22*(max(abs(epspl(i))-epsd12,zero)))
208 ENDIF
209 ! Rate dependency for shear plane 12 modulus
210 IF (ltype12 == 1) THEN
211 f12(i) = d12*(abs(epspl(i))/epsd12)**n12
212 ELSEIF (ltype12 == 2) THEN
213 f12(i) = d12*(abs(epspl(i))/epsd12) + n12
214 ELSEIF (ltype12 == 3) THEN
215 f12(i) = d12*log(max(abs(epspl(i))/epsd12,one)) + log(n12)
216 ELSEIF (ltype12 == 4) THEN
217 f12(i) = d12*tanh(n12*(max(abs(epspl(i))-epsd12,zero)))
218 ENDIF
219 ! Rate dependency in fiber direction 1 for Young modulus
220 IF (ltype11 == 1) THEN
221 f11r(i) = d11u*(abs(epspl(i))/epsd11)**n11u
222 ELSEIF (ltype11 == 2) THEN
223 f11r(i) = d11u*(abs(epspl(i))/epsd11) + n11u
224 ELSEIF (ltype11 == 3) THEN
225 f11r(i) = d11u*log(max(abs(epspl(i))/epsd11,one)) + log(n11u)
226 ELSEIF (ltype11 == 4) THEN
227 f11r(i) = d11u*tanh(n11u*(max(abs(epspl(i))-epsd11,zero)))
228 ENDIF
229 ! Rate dependency in fiber direction 1 for Young modulus
230 IF (ltyper0 == 1) THEN
231 fr0(i) = dr0*(abs(epspl(i))/epsdr0)**nr0
232 ELSEIF (ltyper0 == 2) THEN
233 fr0(i) = dr0*(abs(epspl(i))/epsdr0) + nr0
234 ELSEIF (ltyper0 == 3) THEN
235 fr0(i) = dr0*log(max(abs(epspl(i))/epsdr0,one)) + log(nr0)
236 ELSEIF (ltyper0 == 4) THEN
237 fr0(i) = dr0*tanh(nr0*(max(abs(epspl(i))-epsdr0,zero)))
238 ENDIF
239 ENDDO
240c
241 ! Elastic parameters, yield stress and strain rate dependency
242 DO i = 1,nel
243 ! Fiber (direction 1)
244 ! -> Tension
245 IF (epsxx(i) >= zero) THEN
246 e1(i) = e10
247 ! -> Compression
248 ELSE
249 e1(i) = e1c/(one + (gamma*e1c*abs(epsxx(i))))
250 ENDIF
251 e1(i) = e1(i)*(one + f11(i))
252 ! Matrix (direction 2)
253 e2(i) = e20*(one + f22(i))
254 ! In-plane shear 12
255 g12(i) = g120*(one + f12(i))
256 ! Plane stress stiffness matrix terms
257 a11(i) = e1(i)/(one - nu12*nu21)
258 a12(i) = nu21*a11(i)
259 a22(i) = e2(i)/(one - nu12*nu21)
260 ! Yield stress
261 sigy(i) = sigy0*(one + fr0(i)) + beta*exp(m*log(pla(i)+em20))
262 ! Rate dependency on fiber failure
263 efti(i) = uvar(i,3)
264 IF (uvar(i,3) == zero) efti(i) = efti0*(one + f11r(i))
265 eftu(i) = uvar(i,4)
266 IF (uvar(i,4) == zero) eftu(i) = eftu0*(one + f11r(i))
267 efci(i) = uvar(i,5)
268 IF (uvar(i,5) == zero) efci(i) = efci0*(one + f11r(i))
269 efcu(i) = uvar(i,6)
270 IF (uvar(i,6) == zero) efcu(i) = efcu0*(one + f11r(i))
271 ! Rate dependency on matrix failure
272 y0(i) = uvar(i,7)
273 IF (uvar(i,7) == zero) y0(i) = y00*sqrt(one + f12(i))
274 yc(i) = uvar(i,8)
275 IF (uvar(i,8) == zero) yc(i) = yc0*sqrt(one + f12(i))
276 y0p(i) = uvar(i,9)
277 IF (uvar(i,9) == zero) y0p(i) = y0p0*sqrt(one + f22(i))
278 ycp(i) = uvar(i,10)
279 IF (uvar(i,10) == zero) ycp(i) = ycp0*sqrt(one + f22(i))
280 y0pc(i) = uvar(i,11)
281 IF (uvar(i,11) == zero) y0pc(i) = y0pc0*sqrt(one + f22(i))
282 ycpc(i) = uvar(i,12)
283 IF (uvar(i,12) == zero) ycpc(i) = ycpc0*sqrt(one + f22(i))
284 ENDDO
285c
286 !========================================================================
287 ! - COMPUTATION OF TRIAL VALUES
288 !========================================================================
289 DO i=1,nel
290c
291 ! Computation of the trial stress tensor
292 signyy(i) = sigoyy(i)/max((one-dp(i)),em20) + a12(i)*depsxx(i) + a22(i)*depsyy(i)
293 signxy(i) = sigoxy(i)/max((one- d(i)),em20) + g12(i)*depsxy(i)
294 signyz(i) = sigoyz(i)/max(min(one-d(i),one-dp(i)),em20) + g23*depsyz(i)*shf(i)
295 signzx(i) = sigozx(i)/max(min(one-d(i),one-dp(i)),em20) + g31*depszx(i)*shf(i)
296C
297 ! Equivalent stress
298 seq(i) = sqrt((signxy(i))**2 + a*(signyy(i))**2)
299C
300 ENDDO
301c
302 !========================================================================
303 ! - COMPUTATION OF YIELD FONCTION
304 !========================================================================
305 phi(1:nel) = seq(1:nel) - sigy(1:nel)
306 ! Checking plastic behavior for all elements
307 nindx = 0
308 index(1:nel) = 0
309 DO i=1,nel
310 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
311 nindx = nindx+1
312 index(nindx) = i
313 ENDIF
314 ENDDO
315c
316 !====================================================================
317 ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM (NEWTON ITERATION)
318 !====================================================================
319c
320 ! Number of maximum Newton iterations
321 niter = 3
322c
323 ! Loop over the iterations
324 DO iter = 1, niter
325#include "vectorize.inc"
326 ! Loop over yielding elements
327 DO ii=1,nindx
328 i = index(ii)
329c
330 ! Note : in this part, the purpose is to compute for each iteration
331 ! a plastic multiplier allowing to update internal variables to satisfy
332 ! the consistency condition using the cutting plane semi-implicit method
333 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
334 ! -> PHI : current value of yield function (known)
335 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
336 ! into account of internal variables kinetic (to compute)
337c
338 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
339 !-------------------------------------------------------------
340 normyy = a*signyy(i)/max(seq(i),em20)
341 normxy = signxy(i)/max(seq(i),em20)
342c
343 ! 2 - Computation of DPHI_DLAMBDA
344 !---------------------------------------------------------
345c
346 ! a) Derivative with respect stress increments tensor DSIG
347 ! --------------------------------------------------------
348 dfdsig2 = normyy*normyy*a22(i) + normxy*normxy*g12(i)
349c
350 ! b) Derivatives with respect to plastic strain P
351 ! ------------------------------------------------
352c
353 ! i) Derivative of the yield stress with respect to plastic strain dSIGY / dPLA
354 ! ----------------------------------------------------------------------------
355 h(i) = beta*m*exp((m-1)*log(pla(i)+em20))
356 h(i) = min(h(i),max(two*g120,e20))
357c
358 ! ii) Derivative of dPLA with respect to DLAM
359 ! -------------------------------------------
360 sig_dfdsig = signyy(i)*normyy + signxy(i)*normxy
361 dpla_dlam(i) = sig_dfdsig/max(sigy(i),em20)
362c
363 ! 3 - Computation of plastic multiplier and variables update
364 !----------------------------------------------------------
365c
366 ! Derivative of PHI with respect to DLAM
367 dphi_dlam(i) = - dfdsig2 - h(i)*dpla_dlam(i)
368 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
369c
370 ! Computation of the plastic multiplier
371 dlam = -phi(i)/dphi_dlam(i)
372c
373 ! Plastic strains tensor increment
374 dpyy(i) = dlam*normyy
375 dpxy(i) = dlam*normxy
376c
377 ! Total plastic strain along matrix direction
378 epspyy(i) = epspyy(i) + dpyy(i)
379c
380 ! Elasto-plastic stresses update
381 signyy(i) = signyy(i) - a22(i)*dpyy(i)
382 signxy(i) = signxy(i) - dpxy(i)*g12(i)
383c
384 ! Cumulated plastic strain and strain rate update
385 ddep = dlam*dpla_dlam(i)
386 dpla(i) = max(zero, dpla(i) + ddep)
387 pla(i) = pla(i) + ddep
388c
389 ! Update equivalent stress
390 seq(i) = sqrt((signxy(i))**2 + a*(signyy(i))**2)
391c
392 ! Transverse strain update
393 dpzz(i) = dpzz(i) - dpyy(i)
394c
395 ! Update the hardening yield stress
396 sigy(i) = sigy(i) + h(i)*dlam*dpla_dlam(i)
397c
398 ! Update yield function value
399 phi(i) = seq(i) - sigy(i)
400c
401 ENDDO
402 ! End of the loop over the yielding elements
403c
404 ENDDO
405 ! End of the loop over the iterations
406 !===================================================================
407 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
408 !===================================================================
409c
410 !===================================================================
411 ! - DAMAGE VARIABLES COMPUTATION AND UPDATE STRESS TENSOR
412 !===================================================================
413 DO i = 1,nel
414c
415 ! Fiber damage (direction 1)
416 ! ------------------------------------------
417 ! -> Fiber equivalent strain
418 epsf_eq = epsxx(i) + nu21*(epsyy(i)-epspyy(i))
419 ! -> Tension
420 IF (epsf_eq >= zero) THEN
421 IF (epsf_eq < efti(i)) THEN
422 dft(i) = max(zero,dft(i))
423 ELSEIF ((epsf_eq >= efti(i)).AND.(epsf_eq < eftu(i))) THEN
424 dft(i) = max(dftu*((epsf_eq-efti(i))/(eftu(i)-efti(i))),dft(i))
425 ! Save damage threshold in case of strain rate dependency
426 IF (uvar(i,3) == zero) uvar(i,3) = efti(i)
427 IF (uvar(i,4) == zero) uvar(i,4) = eftu(i)
428 ELSEIF (epsf_eq >= eftu(i)) THEN
429 dft(i) = max(one - (one - dftu)*(eftu(i)/epsf_eq),dft(i))
430 ENDIF
431 dft(i) = max(dft(i),zero)
432 dft(i) = min(dft(i),one)
433 df(i) = dft(i)
434 ! -> Compression
435 ELSEIF ((epsf_eq < zero).AND.(ibuck > 1)) THEN
436 IF (abs(epsf_eq) < efci(i)) THEN
437 dfc(i) = max(zero,dfc(i))
438 ELSEIF ((abs(epsf_eq) >= efci(i)).AND.(abs(epsf_eq) < efcu(i))) THEN
439 dfc(i) = max(dfcu*((abs(epsf_eq)-efci(i))/(efcu(i)-efci(i))),dfc(i))
440 ! Save damage threshold in case of strain rate dependency
441 IF (uvar(i,5) == zero) uvar(i,5) = efci(i)
442 IF (uvar(i,6) == zero) uvar(i,6) = efcu(i)
443 ELSEIF (abs(epsf_eq) >= efcu(i)) THEN
444 dfc(i) = max(one - (one - dfcu)*(efcu(i)/abs(epsf_eq)),dfc(i))
445 ENDIF
446 dfc(i) = max(dfc(i),zero)
447 dfc(i) = min(dfc(i),one)
448 df(i) = dfc(i)
449 ENDIF
450c
451 ! Matrix damage (direction 2 and 3)
452 ! ------------------------------------------
453 ! -> Damage functions (derivatives of elastic energy)
454 zd(i) = half*((signxy(i)**2/g120) + (signzx(i)**2/g31))
455 ! -> If compression damage is activated
456 IF (icomp > 0) THEN
457 zdp(i) = half*((signyy(i)**2)/e20)
458 ! -> Tension only
459 ELSE
460 zdp(i) = half*(((max(signyy(i),zero))**2)/e20)
461 ENDIF
462 ! -> Damage evolution functions
463 y(i) = max(y(i) ,sqrt(zd(i)+b*zdp(i)))
464 yp(i) = max(yp(i),sqrt(zdp(i)))
465 ENDDO
466c
467 ! -> Shear damage
468 ! Linear
469 IF (ish == 1) THEN
470 DO i = 1,nel
471 IF (y(i)<y0(i)) THEN
472 d(i) = zero
473 ELSEIF ((d(i)<dmax).AND.(y(i)<ysp).AND.(y(i)<yr)) THEN
474 d(i) = max(y(i)-y0(i),zero)/max(yc(i),em20)
475 d(i) = min(d(i),dmax)
476 ! Save damage threshold in case of strain rate dependency
477 IF (uvar(i,7) == zero) uvar(i,7) = y0(i)
478 IF (uvar(i,8) == zero) uvar(i,8) = yc(i)
479 ELSE
480 d(i) = one - (one - dmax)*uvar(i,1)/max(y(i),em20)
481 ENDIF
482 d(i) = max(d(i),zero)
483 d(i) = min(d(i), one)
484 ENDDO
485 ! Exponential
486 ELSEIF (ish == 2) THEN
487 DO i = 1,nel
488 IF (y(i)>y0(i)) THEN
489 d(i) = dsat1*(one - exp((y0(i)-y(i))/max(yc(i),em20)))
490 ! Save damage threshold in case of strain rate dependency
491 IF (uvar(i,7) == zero) uvar(i,7) = y0(i)
492 IF (uvar(i,8) == zero) uvar(i,8) = yc(i)
493 ELSE
494 d(i) = zero
495 ENDIF
496 d(i) = max(d(i),zero)
497 d(i) = min(d(i), one)
498 ENDDO
499 ! Tabulated function
500 ELSEIF (ish == 3) THEN
501 ipos(1:nel) = vartmp(1:nel,1)
502 iad(1:nel) = npf(ifunc(1)) / 2 + 1
503 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
504 DO i = 1,nel
505 var(i) = y(i)/y0(i)
506 ENDDO
507 CALL vinter(tf,iad,ipos,ilen,nel,var,dydx,d)
508 vartmp(1:nel,1) = ipos(1:nel)
509 DO i = 1,nel
510 ! Save damage threshold in case of strain rate dependency
511 IF (uvar(i,7) == zero .AND. d(i) /= zero) uvar(i,7) = y0(i)
512 d(i) = max(d(i),zero)
513 d(i) = min(d(i), one)
514 ENDDO
515 ENDIF
516c
517 ! -> Transverse damage
518 ! Linear
519 IF (itr == 1) THEN
520 DO i = 1,nel
521 ! -> Compression damage if activated
522 IF ((icomp > 0).AND.(epsyy(i)<zero)) THEN
523 IF (yp(i)<y0pc(i)) THEN
524 dp(i) = zero
525 ELSEIF ((dp(i)<dmax).AND.(yp(i)<ysp).AND.(yp(i)<yr)) THEN
526 dp(i) = max(yp(i)-y0pc(i),zero)/max(ycpc(i),em20)
527 dp(i) = min(dp(i),dmax)
528 ! Save damage threshold in case of strain rate dependency
529 IF (uvar(i,11) == zero) uvar(i,11) = y0pc(i)
530 IF (uvar(i,12) == zero) uvar(i,12) = ycpc(i)
531 ELSE
532 dp(i) = one - (one - dmax)*uvar(i,2)/max(yp(i),em20)
533 ENDIF
534 ! -> Tension damage
535 ELSE
536 IF (yp(i)<y0p(i)) THEN
537 dp(i) = zero
538 ELSEIF ((dp(i)<dmax).AND.(yp(i)<ysp).AND.(yp(i)<yr)) THEN
539 dp(i) = max(yp(i)-y0p(i),zero)/max(ycp(i),em20)
540 dp(i) = min(dp(i),dmax)
541 ! Save damage threshold in case of strain rate dependency
542 IF (uvar(i,9) == zero) uvar(i,9) = y0p(i)
543 IF (uvar(i,10) == zero) uvar(i,10) = ycp(i)
544 ELSE
545 dp(i) = one - (one - dmax)*uvar(i,2)/max(yp(i),em20)
546 ENDIF
547 ENDIF
548 dp(i) = max(dp(i),zero)
549 dp(i) = min(dp(i), one)
550 ENDDO
551 ! Exponential
552 ELSEIF (itr == 2) THEN
553 DO i = 1,nel
554 ! -> Compression damage if activated
555 IF ((icomp > 0).AND.(epsyy(i)<zero)) THEN
556 IF (yp(i)>y0pc(i)) THEN
557 dp(i) = dsat2c*(one - exp((y0pc(i)-yp(i))/max(ycpc(i),em20)))
558 ! Save damage threshold in case of strain rate dependency
559 IF (uvar(i,11) == zero) uvar(i,11) = y0pc(i)
560 IF (uvar(i,12) == zero) uvar(i,12) = ycpc(i)
561 ELSE
562 dp(i) = zero
563 ENDIF
564 ! -> Tension damage
565 ELSE
566 IF (yp(i)>y0p(i)) THEN
567 dp(i) = dsat2*(one - exp((y0p(i)-yp(i))/max(ycp(i),em20)))
568 ! Save damage threshold in case of strain rate dependency
569 IF (uvar(i,9) == zero) uvar(i,9) = y0p(i)
570 IF (uvar(i,10) == zero) uvar(i,10) = ycp(i)
571 ELSE
572 dp(i) = zero
573 ENDIF
574 ENDIF
575 dp(i) = max(dp(i),zero)
576 dp(i) = min(dp(i), one)
577 ENDDO
578 ! Tabulated function
579 ELSEIF (itr == 3) THEN
580 ! -> Compression tabulated damage
581 IF (icomp > 0) THEN
582 ipos(1:nel) = vartmp(1:nel,3)
583 iad(1:nel) = npf(ifunc(3)) / 2 + 1
584 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
585 DO i = 1,nel
586 var(i) = yp(i)/y0pc(i)
587 ENDDO
588 CALL vinter(tf,iad,ipos,ilen,nel,var,dydx,dptc)
589 vartmp(1:nel,3) = ipos(1:nel)
590 ENDIF
591 ! -> Tension tabulated damage
592 ipos(1:nel) = vartmp(1:nel,2)
593 iad(1:nel) = npf(ifunc(2)) / 2 + 1
594 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
595 DO i = 1,nel
596 var(i) = yp(i)/y0p(i)
597 ENDDO
598 CALL vinter(tf,iad,ipos,ilen,nel,var,dydx,dpt)
599 vartmp(1:nel,2) = ipos(1:nel)
600 DO i = 1,nel
601 ! -> Compression damage if activated
602 IF ((icomp > 0).AND.(epsyy(i)<zero)) THEN
603 dp(i) = dptc(i)
604 ! Save damage threshold in case of strain rate dependency
605 IF ((uvar(i,11) == zero) .AND. (dp(i) /= zero)) uvar(i,11) = y0pc(i)
606 ! -> Tension damage
607 ELSE
608 dp(i) = dpt(i)
609 ! Save damage threshold in case of strain rate dependency
610 IF ((uvar(i,9) == zero) .AND. (dp(i) /= zero)) uvar(i,9) = y0p(i)
611 ENDIF
612 dp(i) = max(dp(i),zero)
613 dp(i) = min(dp(i), one)
614 ENDDO
615 ENDIF
616c
617 DO i = 1,nel
618c
619 ! Sound-speed
620 soundsp(i) = sqrt(max(a11(i),a22(i))/rho(i))
621c
622 ! Compute damaged plane stress stiffness matrix
623 a11(i) = a11(i)*(one - df(i))
624 a12(i) = nu21*a11(i)*(one - dp(i))
625 a22(i) = a22(i)*(one - dp(i))
626c
627 ! Update stresses with damage softening effect
628 ! --------------------------------------------
629 ! -> If non-linear compression young modulus is used for fibers
630 IF (gamma > zero .AND. epsxx(i) < zero) THEN
631 signxx(i) = -(one/gamma)*log(one + gamma*e1c*abs(epsxx(i)))*(one - df(i))*(one + f11(i))
632 ! -> If linear elasticity is used for fibers
633 ELSE
634 signxx(i) = a11(i)*epsxx(i)
635 ENDIF
636 signxx(i) = signxx(i) + a12(i)*(epsyy(i) - epspyy(i))
637 signyy(i) = a12(i)*epsxx(i) + a22(i)*(epsyy(i) - epspyy(i))
638 signxy(i) = signxy(i)*(one - d(i))
639 signyz(i) = signyz(i)*min(one - d(i),one - dp(i))
640 signzx(i) = signzx(i)*min(one - d(i),one - dp(i))
641c
642 ! -> Store internal variables
643 ! ------------------------------------------
644 dmg(i,1) = max(df(i),d(i),dp(i))
645 dmg(i,2) = df(i)
646 dmg(i,3) = d(i)
647 dmg(i,4) = dp(i)
648 dmg(i,5) = dft(i)
649 dmg(i,6) = dfc(i)
650 uvar(i,1) = y(i)
651 uvar(i,2) = yp(i)
652 uvar(i,16) = epspyy(i)
653c
654 ENDDO
655 !===================================================================
656c
657 ! Sound-speed and thickness update
658 DO i=1,nel
659 ! Coefficient for hourglass
660 IF (dpla(i)>zero) THEN
661 etse(i) = h(i) / (h(i) + max(e1(i),e2(i)))
662 ELSE
663 etse(i) = one
664 ENDIF
665 ! Computation of the thickness variation
666 deelzz(i) = -(nu13/e1(i))*(signxx(i)-sigoxx(i))-(nu23/e2(i))*(signyy(i)-sigoyy(i))
667 dezz(i) = deelzz(i) + dpzz(i)
668 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
669 ENDDO
670C
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