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