OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat122c_newton.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_newton ../engine/source/materials/mat/mat122/mat122c_newton.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_newton(
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 . 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
671 END
672C
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat122c_newton
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)
Definition mat122c_newton.F:40
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72