OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps124.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!|| sigeps124 ../engine/source/materials/mat/mat124/sigeps124.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
29!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.f
30!||====================================================================
31 SUBROUTINE sigeps124(
32 1 NEL ,NGL ,NUPARAM ,NUVAR ,TIMESTEP,TIME ,
33 2 UPARAM ,UVAR ,RHO ,PLA ,DPLA ,
34 3 SOUNDSP ,EPSD ,OFF ,LOFF ,IPG ,NPG ,
35 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
36 5 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
37 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
38 7 SIGY ,ET ,DMG ,LE )
39 !=======================================================================
40 ! Implicit types
41 !=======================================================================
42#include "implicit_f.inc"
43 !=======================================================================
44 ! Common
45 !=======================================================================
46#include "com08_c.inc"
47#include "units_c.inc"
48#include "comlock.inc"
49#include "scr05_c.inc"
50#include "scr17_c.inc"
51#include "mvsiz_p.inc"
52 !=======================================================================
53 ! Dummy arguments
54 !=======================================================================
55 INTEGER,INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,NPG,NGL(NEL)
56 my_real
57 . TIME,TIMESTEP,UPARAM(NUPARAM)
58 my_real,DIMENSION(NEL), INTENT(IN) ::
59 . RHO,LE,
60 . DEPSXX,DEPSYY,DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,
61 . EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX
62 my_real ,DIMENSION(NEL), INTENT(OUT) ::
63 . soundsp,signxx,signyy,signzz,signxy,signyz,signzx
64 my_real,DIMENSION(NEL) ::
65 . sigy,et
66 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
67 . pla,dpla,epsd,loff,off
68 my_real ,DIMENSION(NEL,3), INTENT(INOUT) ::
69 . dmg
70 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
71 . uvar
72 !=======================================================================
73 ! Local Variables
74 !=======================================================================
75 INTEGER I,J,II,ITER,NITER,NINDX,INDX(NEL),NINDX2,INDX2(NEL),
76 . IRATE,DTYPE,DFLAG,IREG,INDX3(NEL),NINDX3,IDEL
77 my_real
78 . young,bulk,g,nu,lam,g2,afiltr,ah,bh,ch,dh,hp,as,qh0,
79 . m0,ecc,df,bs,epsi,epst0,epstmax,deltas,betas,epsc0,
80 . epscmax,alphas,gammas
81 my_real
82 . dlam,ag,bg,dfp_dkap,dfp_dqh1,dfp_dqh2,dfp_drhob,dfp_dsigm,
83 . normxx,normyy,normzz,normxy,normyz,normzx,denom,dfpdsig2,
84 . dfp_dtheta,dgp_drhob,dgp_dsigm,dj3dsxx,dj3dsxy,dj3dsyy,
85 . dj3dsyz,dj3dszx,dj3dszz,dkap_dlam,dmg_dsigm,dqh1_dkap,
86 . dqh2_dkap,drcos_dtheta,dtheta_dj2,dtheta_dj3,eh,fh,
87 . normdp,normgp,normpxx,normpxy,normpyy,normpyz,normpzx,
88 . normpzz,rh,trdep,trdgpds,xsih,normfp,ldav,dfapex_dkap,dfp,
89 . rs,u,uprim,v,vprim,dfp_dth,dr_dcosth,dth_dj2,dth_dj3,
90 . wf,wf1,efc,betac,sigtens(mvsiz,6),evv(mvsiz,3),dirprv(mvsiz,3,3),
91 . sigpt(nel,3),sigpc(nel,3)
92 my_real, DIMENSION(NEL) ::
93 . sxx,syy,szz,sxy,syz,szx,j2,j3,theta,qh1,qh2,hardp,fp,trsig,
94 . dpxx,dpyy,dpxy,dpzz,dpyz,dpzx,rcos,kap,sigm,rhob,dfp_dlam,
95 . costh,sinth,cos3th,al,bl,cl,plaxx,playy,plazz,plaxy,playz,plazx,
96 . apex,fapex,eps_eq,kapdt,kapdt1,kapdt2,fst,omegat,eps_inel,h,fc,ft,
97 . alphart,alpharc,arate,eps0,omegac,kapdc,kapdc1,xsis,normsigp,
98 . kapdc2,fsc,alphac,sigtxx,sigtyy,sigtzz,sigtxy,sigtyz,sigtzx,
99 . sigcxx,sigcyy,sigczz,sigcxy,sigcyz,sigczx,elxx,elyy,elzz,elxy,
100 . elyz,elzx,sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx,ft1
101 !=======================================================================
102c
103 !=======================================================================
104 ! - INITIALISATION OF COMPUTATION ON TIME STEP
105 !=======================================================================
106 ! Recovering model parameters
107 ! -> Elastic parameters
108 young = uparam(1) ! Young modulus
109 nu = uparam(2) ! Poisson ratio
110 g = uparam(3) ! Shear modulus
111 g2 = uparam(4) ! 2*Shear modulus
112 lam = uparam(5) ! Lame coefficient
113 bulk = uparam(6) ! Bulk modulus
114 ! -> Plastic parameters
115 ft(1:nel) = uparam(7) ! Uniaxial tensile strength
116 fc(1:nel) = uparam(8) ! Uniaxial compressive strength
117 ecc = uparam(9) ! Eccentricity
118 m0 = uparam(10) ! Friction parameter
119 qh0 = uparam(11) ! Initial hardening
120 hp = uparam(12) ! Hardening modulus
121 ah = uparam(13) ! Hardening ductility parameter 1
122 bh = uparam(14) ! Hardening ductility parameter 2
123 ch = uparam(15) ! Hardening ductility parameter 3
124 dh = uparam(16) ! Hardening ductility parameter 4
125 ! -> Damage parameters
126 as = uparam(17) ! Damage ductility measure
127 bs = uparam(18) ! Damage ductility parameter
128 df = uparam(19) ! Dilation parameter
129 dflag = nint(uparam(20)) ! Damage flag
130 dtype = nint(uparam(21)) ! Tensile damage type
131 ireg = nint(uparam(22)) ! Regularization flag
132 wf = uparam(23) ! First displacement threshold
133 wf1 = uparam(24) ! Second displacement threshold
134 ft1(1:nel) = uparam(25) ! Second uniaxial tensile strength
135 efc = uparam(26) ! Compressive inelastic strain threshold
136 ! -> Strain rate effect parameters
137 irate = nint(uparam(27)) ! Strain rate effect flag
138 epst0 = uparam(28) ! Reference tensile strain rate
139 epstmax = uparam(29) ! Maximum tensile strain rate threshold
140 deltas = uparam(30) ! Tensile strain rate effect parameter 1
141 betas = uparam(31) ! Tensile strain rate effect parameter 2
142 epsc0 = uparam(32) ! Reference tensile strain rate
143 epscmax = uparam(33) ! Maximum compressive strain rate threshold
144 alphas = uparam(34) ! Compressive strain rate effect parameter 1
145 gammas = uparam(35) ! Compressive strain rate effect parameter 2
146 ! -> Element deletion flag
147 idel = nint(uparam(36)) ! Element deletion flag
148c
149 !--------------------------------------------------------------------
150 ! recovering internal variables
151 !--------------------------------------------------------------------
152 DO i=1,nel
153 ! Elastic strain tensor
154 elxx(i) = (epsxx(i)-depsxx(i)) - uvar(i,2)
155 elyy(i) = (epsyy(i)-depsyy(i)) - uvar(i,3)
156 elzz(i) = (epszz(i)-depszz(i)) - uvar(i,4)
157 elxy(i) = (epsxy(i)-depsxy(i)) - uvar(i,5)
158 elyz(i) = (epsyz(i)-depsyz(i)) - uvar(i,6)
159 elzx(i) = (epszx(i)-depszx(i)) - uvar(i,7)
160 ldav = (elxx(i) + elyy(i) + elzz(i)) * lam
161 ! Old undamaged stress tensor
162 sigoxx(i) = elxx(i)*g2 + ldav
163 sigoyy(i) = elyy(i)*g2 + ldav
164 sigozz(i) = elzz(i)*g2 + ldav
165 sigoxy(i) = elxy(i)*g
166 sigoyz(i) = elyz(i)*g
167 sigozx(i) = elzx(i)*g
168 ! User variables
169 kap(i) = uvar(i,1) ! Kappa hardening parameter
170 plaxx(i) = uvar(i,2) ! Plastic strain tensor component X
171 playy(i) = uvar(i,3) ! Plastic strain tensor component Y
172 plazz(i) = uvar(i,4) ! Plastic strain tensor component Z
173 plaxy(i) = uvar(i,5) ! 2*Plastic strain tensor component XY
174 playz(i) = uvar(i,6) ! 2*Plastic strain tensor component YZ
175 plazx(i) = uvar(i,7) ! 2*Plastic strain tensor component ZX
176 kapdt(i) = uvar(i,9) ! Tensile damage strain threshold
177 kapdt1(i) = uvar(i,10) ! First tensile damage variable
178 kapdt2(i) = uvar(i,11) ! Second tensile damage variable
179 kapdc(i) = uvar(i,12) ! Compressive damage strain threshold
180 kapdc1(i) = uvar(i,13) ! First compressive damage variable
181 kapdc2(i) = uvar(i,14) ! Second compressive damage variable
182 arate(i) = uvar(i,15) ! Rate dependency factor
183 ! Damage variables
184 omegat(i) = dmg(i,2) ! Tensile damage
185 omegac(i) = dmg(i,3) ! Compressive damage
186 ! Standard inputs
187 dpla(i) = zero ! Initialization of the plastic strain increment
188 et(i) = one ! Initialization of hourglass coefficient
189 hardp(i) = zero ! Initialization of hardening rate
190 dpxx(i) = zero ! Initialization of the XX plastic strain increment
191 dpyy(i) = zero ! Initialization of the YY plastic strain increment
192 dpzz(i) = zero ! Initialization of the ZZ plastic strain increment
193 dpxy(i) = zero ! Initialization of the XY plastic strain increment
194 dpyz(i) = zero ! Initialization of the YZ plastic strain increment
195 dpzx(i) = zero ! Initialization of the ZX plastic strain increment
196 ! Regularization factor
197 IF (ireg > 1) THEN
198 IF (uvar(i,16) == zero) uvar(i,16) = le(i)
199 h(i) = uvar(i,16)
200 ELSE
201 h(i) = one
202 ENDIF
203 ENDDO
204c
205 !========================================================================
206 ! - COMPUTATION OF TRIAL VALUES + YIELD FONCTION
207 !========================================================================
208 ! Computation of the trial stress tensor
209 DO i=1,nel
210 IF (loff(i) == one) THEN
211 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
212 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
213 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
214 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
215 signxy(i) = sigoxy(i) + depsxy(i)*g
216 signyz(i) = sigoyz(i) + depsyz(i)*g
217 signzx(i) = sigozx(i) + depszx(i)*g
218 ELSE
219 signxx(i) = zero
220 signyy(i) = zero
221 signzz(i) = zero
222 signxy(i) = zero
223 signyz(i) = zero
224 signzx(i) = zero
225 ENDIF
226 ENDDO
227c
228 ! Strain rate effects
229 IF (irate > 1) THEN
230 ! Computation of principal stresses and vectors
231 DO i = 1,nel
232 sigtens(i,1) = signxx(i)
233 sigtens(i,2) = signyy(i)
234 sigtens(i,3) = signzz(i)
235 sigtens(i,4) = signxy(i)
236 sigtens(i,5) = signyz(i)
237 sigtens(i,6) = signzx(i)
238 ENDDO
239 evv(1:mvsiz,1:3) = zero
240 dirprv(1:mvsiz,1:3,1:3) = zero
241 sigpt(1:nel,1:3) = zero
242 sigpc(1:nel,1:3) = zero
243 IF (iresp == 1) THEN
244 CALL valpvecdp_v(sigtens,evv,dirprv,nel)
245 ELSE
246 CALL valpvec_v(sigtens,evv,dirprv,nel)
247 ENDIF
248c
249 ! Computation of ALPHAC compression factor
250 DO i = 1,nel
251 normsigp(i) = evv(i,1)**2 + evv(i,2)**2 + evv(i,3)**2
252 ENDDO
253 alphac(1:nel) = zero
254 DO j = 1,3
255 DO i = 1,nel
256 IF (evv(i,j) > zero) THEN
257 sigpt(i,j) = evv(i,j)
258 ELSE
259 sigpc(i,j) = evv(i,j)
260 ENDIF
261 alphac(i) = alphac(i) + sigpc(i,j)*(sigpt(i,j) + sigpc(i,j))/max(normsigp(i),em20)
262 ENDDO
263 ENDDO
264c
265 ! computation of strain rate factor
266 DO i = 1,nel
267 ! Tension strain rate factor
268 IF (epsd(i) <= epst0) THEN
269 alphart(i) = one
270 ELSEIF (epsd(i) > epst0 .AND. epsd(i) <= epstmax) THEN
271 alphart(i) = (epsd(i)/epst0)**(deltas)
272 ELSE
273 alphart(i) = betas*(epsd(i)/epst0)**third
274 ENDIF
275c
276 ! Compression strain rate factor
277 IF (epsd(i) <= epsc0) THEN
278 alpharc(i) = one
279 ELSEIF (epsd(i) > epsc0 .AND. epsd(i) <= epscmax) THEN
280 alpharc(i) = (epsd(i)/epsc0)**(1.026d0*alphas)
281 ELSE
282 alpharc(i) = gammas*(epsd(i)/epsc0)**third
283 ENDIF
284c
285 ! Rate factor
286 IF (arate(i) == zero) arate(i) = (one - alphac(i))*alphart(i) + alphac(i)*alpharc(i)
287c
288 ! Compute strain rate dependent parameters
289 ft(i) = ft(i)*arate(i)
290 ft1(i) = ft1(i)*arate(i)
291 fc(i) = fc(i)*arate(i)
292 ENDDO
293 ENDIF
294c
295 ! Computation of the yield functions
296 DO i = 1,nel
297 ! Computation of the trace of stress tensor and hydrostatic stress
298 trsig(i) = signxx(i) + signyy(i) + signzz(i)
299 sigm(i) = trsig(i)*third
300 ! Computation of the deviatoric trial stress tensor
301 sxx(i) = signxx(i) - sigm(i)
302 syy(i) = signyy(i) - sigm(i)
303 szz(i) = signzz(i) - sigm(i)
304 sxy(i) = signxy(i)
305 syz(i) = signyz(i)
306 szx(i) = signzx(i)
307 ! Second deviatoric invariant
308 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2)
309 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
310 j2(i) = max(j2(i),em20)
311 ! Radius, Third deviatoric invariant and Lode angle
312 rhob(i) = sqrt(two*j2(i))
313 j3(i) = sxx(i)*syy(i)*szz(i) + two*sxy(i)*szx(i)*syz(i) -
314 . sxx(i)*(syz(i)**2) - szz(i)*(sxy(i)**2) - syy(i)*(szx(i)**2)
315 cos3th(i) = (three*sqr3/two)*j3(i)/((j2(i)**(three/two)))
316 IF (cos3th(i) > one) THEN
317 cos3th(i) = one
318 ELSE IF (cos3th(i) < -one) THEN
319 cos3th(i) = -one
320 ENDIF
321 theta(i) = third*acos(cos3th(i))
322 ! Save cosinus and sinus values to keep good performances
323 costh(i) = cos(theta(i))
324 sinth(i) = sin(theta(i))
325c
326 ! Compute the R function value for Lode Angle
327 rcos(i) = (four*(one - ecc**2)*(costh(i)**2) + (two*ecc - one)**2)/
328 . max(((two*(one - ecc**2)*costh(i)) +
329 . (two*ecc - one)*sqrt(four*(one - ecc**2)*(costh(i)**2) +
330 . five*(ecc**2) - four*ecc)),em20)
331c
332 ! Compute hardening variables QH1 and QH2
333 IF (kap(i) < zero) THEN
334 qh1(i) = qh0
335 qh2(i) = one
336 ELSEIF (kap(i) >= zero .AND. kap(i) < one) THEN
337 qh1(i) = qh0 +
338 . (one - qh0)*((kap(i)**3) - three*(kap(i)**2) + three*kap(i)) -
339 . hp*((kap(i)**3) - three*(kap(i)**2) + two*kap(i))
340 qh2(i) = one
341 ELSE
342 qh1(i) = one
343 qh2(i) = one + (kap(i)-one)*hp
344 ENDIF
345c
346 ! Apex stress
347 apex(i) = qh2(i)*fc(i)/m0
348c
349 ! Compute the regular yield function value
350 cl(i) = (rhob(i)*rcos(i)/(sqr6*fc(i))) + sigm(i)/fc(i)
351 bl(i) = (sigm(i)/fc(i)) + rhob(i)/(sqr6*fc(i))
352 al(i) = (one - qh1(i))*(bl(i)**2) + sqrt(three/two)*(rhob(i)/fc(i))
353 fp(i) = al(i)**2 + m0*(qh1(i)**2)*qh2(i)*cl(i) - (qh1(i)**2)*(qh2(i)**2)
354c
355 ! Compute the apex yield function value
356 fapex(i) = sigm(i) - apex(i)
357 ENDDO
358c
359 ! Checking plastic behavior for all elements
360 nindx = 0
361 indx(1:nel) = 0
362 nindx2 = 0
363 indx2(1:nel) = 0
364 DO i=1,nel
365 ! Apex return mapping
366 IF ((fapex(i) > zero) .AND. (loff(i) == one)) THEN
367 nindx = nindx+1
368 indx(nindx) = i
369 ! Regular return mapping
370 ELSEIF ((fp(i) > zero) .AND. (loff(i) == one)) THEN
371 nindx2 = nindx2+1
372 indx2(nindx2) = i
373 ENDIF
374 ENDDO
375c
376 !====================================================================
377 ! - APEX RETURN MAPPING WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
378 !====================================================================
379 IF (nindx > 0) THEN
380c
381 ! Number of Newton iterations
382 niter = 6
383c
384 ! Loop over the iterations
385 DO iter = 1, niter
386#include "vectorize.inc"
387 ! Loop over yielding elements
388 DO ii=1,nindx
389 i = indx(ii)
390 ! -> Derivative of QH2 with respect to KAP
391 IF (kap(i) < one) THEN
392 dqh2_dkap = zero
393 ELSE
394 dqh2_dkap = hp
395 ENDIF
396 hardp(i) = dqh2_dkap
397 ! Compute the ductility measure
398 rh = -(sigm(i)/fc(i)) - third
399 IF (rh >= zero) THEN
400 xsih = ah-(ah-bh)*exp(-rh/ch)
401 ELSE
402 eh = bh - dh
403 fh = (bh-dh)*ch/(ah-bh)
404 xsih = eh*exp(rh/fh) + dh
405 ENDIF
406 ! Derivative of FAPEX with respect to kappa
407 dfapex_dkap = dqh2_dkap*(fc(i)/m0) + three*bulk*xsih
408 ! Update Kappa parameter
409 kap(i) = kap(i) + fapex(i)/dfapex_dkap
410 ! Update hydrostatic stress
411 sigm(i) = sigm(i) - three*bulk*xsih*fapex(i)/dfapex_dkap
412 ! Update APEX stress
413 IF (kap(i) < one) THEN
414 qh2(i) = one
415 ELSE
416 qh2(i) = one + (kap(i)-one)*hp
417 ENDIF
418 apex(i) = qh2(i)*fc(i)/m0
419 ! Update the APEX yield function
420 fapex(i) = sigm(i) - apex(i)
421 ENDDO
422 ENDDO
423c
424 ! Update the plastic strain and the stress tensors
425#include "vectorize.inc"
426 DO ii=1,nindx
427 i = indx(ii)
428c
429 ! Plastic strain
430 plaxx(i) = plaxx(i) + (one/young)*(
431 & (signxx(i)-sigm(i))
432 & - nu*(signyy(i)-sigm(i))
433 & - nu*(signzz(i)-sigm(i)))
434 playy(i) = playy(i) + (one/young)*(
435 & - nu*(signxx(i)-sigm(i))
436 & + (signyy(i)-sigm(i))
437 & - nu*(signzz(i)-sigm(i)))
438 plazz(i) = plazz(i) + (one/young)*(
439 & - nu*(signxx(i)-sigm(i))
440 & - nu*(signyy(i)-sigm(i))
441 & + (signzz(i)-sigm(i)))
442 plaxy(i) = plaxy(i) + (one/g)*(signxy(i))
443 playz(i) = playz(i) + (one/g)*(signyz(i))
444 plazx(i) = plazx(i) + (one/g)*(signzx(i))
445c
446 ! Stress tensor
447 signxx(i) = sigm(i)
448 signyy(i) = sigm(i)
449 signzz(i) = sigm(i)
450 signxy(i) = zero
451 signyz(i) = zero
452 signzx(i) = zero
453 j2(i) = em20
454 rhob(i) = sqrt(two*j2(i))
455 cl(i) = sigm(i)/fc(i)
456 ENDDO
457c
458 ENDIF
459c
460 !====================================================================
461 ! - REGULAR RETURN MAPPING WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
462 !====================================================================
463 IF (nindx2 > 0) THEN
464c
465 ! Number of Newton iterations
466 niter = 6
467c
468 ! Loop over the iterations
469 DO iter = 1, niter
470#include "vectorize.inc"
471 ! Loop over yielding elements
472 DO ii=1,nindx2
473 i = indx2(ii)
474c
475 ! Note : in this part, the purpose is to compute for each iteration
476 ! a plastic multiplier allowing to update internal variables to satisfy
477 ! the consistency condition using the cutting plane semi-implicit
478 ! method
479 ! Its expression at each iteration is : DLAMBDA = - FP/DFP_DLAMBDA
480 ! -> FP : current value of yield function (known)
481 ! -> DFP_DLAMBDA : derivative of FP with respect to DLAMBDA by taking
482 ! into account of internal variables kinetic :
483 ! hardening parameter ...
484c
485 ! 1 a) - Computation of DFP_DSIG the normal to the yield surface
486 !-------------------------------------------------------------
487 ! -> Computation of derivative of FP with respect to RHOB
488 dfp_drhob = (two*al(i)/fc(i))*(two*((one-qh1(i))/sqr6)*bl(i) + sqrt(three/two)) +
489 . m0*(qh1(i)**2)*qh2(i)*rcos(i)/(sqr6*fc(i))
490 ! -> Computation of derivative of FP with respect to SIGM
491 dfp_dsigm = (four*al(i)*bl(i)*(one-qh1(i))/fc(i)) + m0*(qh1(i)**2)*qh2(i)/fc(i)
492 ! -> Computation of derivative of FP with respect to COSTH
493 u = four*(one - ecc**2)*(costh(i)**2) + (two*ecc - one)**2
494 uprim = eight*(one - ecc**2)*costh(i)
495 v = two*(one - ecc**2)*costh(i) +
496 . (two*ecc - one)*sqrt(four*(one-ecc**2)*(costh(i)**2)
497 . + five*(ecc**2) - four*ecc)
498 vprim = two*(one-ecc**2) + (two*ecc - one)*((eight*(one-ecc**2)*costh(i))/
499 . (two*sqrt(four*(one-ecc**2)*(costh(i)**2)
500 . + five*(ecc**2) - four*ecc)))
501 dr_dcosth = (uprim*v - u*vprim)/(v**2)
502 dfp_dth = m0*(qh1(i)**2)*(qh2(i))*(rhob(i)/(sqr6*fc(i)))*dr_dcosth*(-sinth(i))
503 dth_dj2 = three*sqr3*j3(i)/(four*(j2(i)**2)*sqrt(j2(i))*max(sqrt(one - (cos3th(i))**2),em08))
504 dth_dj3 = -sqr3/(two*j2(i)*sqrt(j2(i))*max(sqrt(one - (cos3th(i))**2),em08))
505c
506 ! Derivative of third deviatoric invariant with respect to stresses
507 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)
508 . - third*(sxx(i)*szz(i)-szx(i)**2)
509 . - third*(sxx(i)*syy(i)-sxy(i)**2)
510 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)
511 . + two_third*(sxx(i)*szz(i)-szx(i)**2)
512 . - third*(sxx(i)*syy(i)-sxy(i)**2)
513 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)
514 . - third*(sxx(i)*szz(i)-szx(i)**2)
515 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)
516 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))
517 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))
518 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))
519c
520 ! -> Assembling derivatives
521 normxx = dfp_drhob*(sxx(i)/rhob(i)) + dfp_dsigm*third +
522 . dfp_dth*(dth_dj2*sxx(i) + dth_dj3*dj3dsxx)
523 normyy = dfp_drhob*(syy(i)/rhob(i)) + dfp_dsigm*third +
524 . dfp_dth*(dth_dj2*syy(i) + dth_dj3*dj3dsyy)
525 normzz = dfp_drhob*(szz(i)/rhob(i)) + dfp_dsigm*third +
526 . dfp_dth*(dth_dj2*szz(i) + dth_dj3*dj3dszz)
527 normxy = two*dfp_drhob*(sxy(i)/rhob(i)) +
528 . dfp_dth*(dth_dj2*two*sxy(i) + dth_dj3*dj3dsxy)
529 normyz = two*dfp_drhob*(syz(i)/rhob(i)) +
530 . dfp_dth*(dth_dj2*two*syz(i) + dth_dj3*dj3dsyz)
531 normzx = two*dfp_drhob*(szx(i)/rhob(i)) +
532 . dfp_dth*(dth_dj2*two*szx(i) + dth_dj3*dj3dszx)
533c
534 ! 1 b) - Computation of DGP_DSIG the normal to the plastic potential
535 !-------------------------------------------------------------
536 ! -> Computation of derivative of GP with respect to RHOB
537 dgp_drhob = (two*al(i)/fc(i))*(two*((one-qh1(i))/sqr6)*bl(i) + sqrt(three/two)) +
538 . (qh1(i)**2)*m0/(sqr6*fc(i))
539 ! -> Computation of AG and BG parameters
540 ag = (three*ft(i)*qh2(i)/fc(i)) + m0*half
541 bg = ((qh2(i)*third)*(one + (ft(i)/fc(i))))/max((
542 . log(ag) - log(two*df - one) - log(three*qh2(i) + m0*half) +
543 . log(df + one)),em20)
544 ! -> computation of derivative of mg with respect to sigm
545 dmg_dsigm = ag*exp((sigm(i) - qh2(i)*ft(i)*third)/(bg*fc(i)))
546 ! -> Computation of derivative of GP with respect to SIGM
547 dgp_dsigm = (four*al(i)*bl(i)*(one-qh1(i))/fc(i)) + ((qh1(i)**2)/fc(i))*dmg_dsigm
548 ! -> Assembling derivatives
549 normpxx = dgp_drhob*(sxx(i)/rhob(i)) + dgp_dsigm*third
550 normpyy = dgp_drhob*(syy(i)/rhob(i)) + dgp_dsigm*third
551 normpzz = dgp_drhob*(szz(i)/rhob(i)) + dgp_dsigm*third
552 normpxy = two*dgp_drhob*(sxy(i)/rhob(i))
553 normpyz = two*dgp_drhob*(syz(i)/rhob(i))
554 normpzx = two*dgp_drhob*(szx(i)/rhob(i))
555c
556 ! 2 - Computation of DPHI_DLAMBDA
557 !---------------------------------------------------------
558c
559 ! a) Derivative with respect stress increments tensor DSIG
560 ! --------------------------------------------------------
561 trdgpds = normpxx + normpyy + normpzz
562 dfpdsig2 = normxx * (normpxx*g2 + lam*trdgpds)
563 . + normyy * (normpyy*g2 + lam*trdgpds)
564 . + normzz * (normpzz*g2 + lam*trdgpds)
565 . + normxy * normpxy*g
566 . + normyz * normpyz*g
567 . + normzx * normpzx*g
568c
569 ! b) Derivatives with respect to hardening Kap
570 ! ------------------------------------------------
571 ! -> Derivative of FP with respect to QH1
572 dfp_dqh1 = -two*al(i)*(bl(i)**2) + two*qh1(i)*qh2(i)*(m0*cl(i) - qh2(i))
573 ! -> Derivative of FP with respect to QH2
574 dfp_dqh2 = (qh1(i)**2)*(m0*cl(i) - two*qh2(i))
575 ! -> Derivative of QH1 and QH2 with respect to KAP
576 IF (kap(i) < zero) THEN
577 dqh1_dkap = (one - qh0)*three - hp*two
578 dqh2_dkap = zero
579 ELSEIF (kap(i) >= zero .AND. kap(i) < one) THEN
580 dqh1_dkap = (one - qh0)*(three*(kap(i)**2) - six*kap(i) + three) -
581 . hp*(three*(kap(i)**2) - six*kap(i) + two)
582 dqh2_dkap = zero
583 ELSE
584 dqh1_dkap = zero
585 dqh2_dkap = hp
586 ENDIF
587 ! -> Derivative of FP with respect to KAP
588 dfp_dkap = dfp_dqh1*dqh1_dkap + dfp_dqh2*dqh2_dkap
589 dfp_dkap = min(dfp_dkap,zero)
590 ! -> Computation of RH and XSIH parameters + Euclidean norm of derivative
591 rh = -(sigm(i)/fc(i)) - third
592 IF (rh >= zero) THEN
593 xsih = ah-(ah-bh)*exp(-rh/ch)
594 ELSE
595 eh = bh-dh
596 fh = (bh-dh)*ch/(ah-bh)
597 xsih = eh*exp(rh/fh) + dh
598 ENDIF
599 normgp = sqrt(third*(dgp_dsigm)**2 + dgp_drhob**2)
600 ! -> Derivative of KAP with respect to DLAM
601 dkap_dlam = (normgp/xsih)*(two*costh(i))**2
602c
603 ! 3 - Computation of plastic multiplier and variables update
604 !----------------------------------------------------------
605c
606 ! Derivative of PHI with respect to DLAM
607 dfp_dlam(i) = - dfpdsig2 + dfp_dkap
608 dfp_dlam(i) = sign(max(abs(dfp_dlam(i)),em20),dfp_dlam(i))
609c
610 ! Computation of the plastic multiplier
611 dlam = -fp(i)/dfp_dlam(i)
612c
613 ! Plastic strains tensor update
614 dpxx(i) = dlam * normpxx
615 dpyy(i) = dlam * normpyy
616 dpzz(i) = dlam * normpzz
617 dpxy(i) = dlam * normpxy
618 dpyz(i) = dlam * normpyz
619 dpzx(i) = dlam * normpzx
620 trdep = dpxx(i) + dpyy(i) + dpzz(i)
621 plaxx(i) = plaxx(i) + dpxx(i)
622 playy(i) = playy(i) + dpyy(i)
623 plazz(i) = plazz(i) + dpzz(i)
624 plaxy(i) = plaxy(i) + dpxy(i)
625 playz(i) = playz(i) + dpyz(i)
626 plazx(i) = plazx(i) + dpzx(i)
627c
628 ! Update the kap hardening parameter KAP
629 kap(i) = kap(i) + dkap_dlam*dlam
630 kap(i) = max(kap(i),zero)
631c
632 ! Update the kap hardening parameters
633! IF (KAP(I) < ZERO) THEN
634! QH1(I) = QH0
635! QH2(I) = ONE
636! ELSEIF (KAP(I) >= ZERO .AND. KAP(I) < ONE) THEN
637! QH1(I) = QH0 +
638! . (ONE - QH0)*((KAP(I)**3) - THREE*(KAP(I)**2) + THREE*KAP(I)) -
639! . HP*((KAP(I)**3) - THREE*(KAP(I)**2) + TWO*KAP(I))
640! QH2(I) = ONE
641! ELSE
642! QH1(I) = ONE
643! qh2(i) = one + (kap(i)-one)*hp
644! ENDIF
645c
646 ! Elasto-plastic stresses update
647 ldav = trdep * lam
648 signxx(i) = signxx(i) - (dpxx(i)*g2 + ldav)
649 signyy(i) = signyy(i) - (dpyy(i)*g2 + ldav)
650 signzz(i) = signzz(i) - (dpzz(i)*g2 + ldav)
651 signxy(i) = signxy(i) - dpxy(i)*g
652 signyz(i) = signyz(i) - dpyz(i)*g
653 signzx(i) = signzx(i) - dpzx(i)*g
654c
655 ! Computation of the trace and hydrostatic stress the new stress tensor
656 trsig(i) = signxx(i) + signyy(i) + signzz(i)
657 sigm(i) = trsig(i)*third
658 ! Computation of the new deviatoric stress tensor
659 sxx(i) = signxx(i) - sigm(i)
660 syy(i) = signyy(i) - sigm(i)
661 szz(i) = signzz(i) - sigm(i)
662 sxy(i) = signxy(i)
663 syz(i) = signyz(i)
664 szx(i) = signzx(i)
665 ! Second deviatoric invariant
666 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2)
667 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
668 j2(i) = max(j2(i),em20)
669 ! Radius, Third deviatoric invariant and Lode angle
670 rhob(i) = sqrt(two*j2(i))
671 j3(i) = sxx(i)*syy(i)*szz(i) + two*sxy(i)*szx(i)*syz(i) -
672 . sxx(i)*(syz(i)**2) - szz(i)*(sxy(i)**2) - syy(i)*(szx(i)**2)
673 cos3th(i) = (three*sqr3/two)*j3(i)/((j2(i)**(three/two)))
674 IF (cos3th(i) > one) THEN
675 cos3th(i) = one
676 ELSE IF (cos3th(i) < -one) THEN
677 cos3th(i) = -one
678 ENDIF
679 theta(i) = third*acos(cos3th(i))
680 ! Save cosinus and sinus values to keep good performances
681 costh(i) = cos(theta(i))
682 sinth(i) = sin(theta(i))
683c
684 ! Compute the R function value for Lode Angle
685 rcos(i) = (four*(one - ecc**2)*(costh(i)**2) + (two*ecc - one)**2)/
686 . max(((two*(one - ecc**2)*costh(i)) +
687 . (two*ecc - one)*sqrt(four*(one - ecc**2)*(costh(i)**2) +
688 . five*(ecc**2) - four*ecc)),em20)
689c
690 ! Update the yield function value
691 cl(i) = (rhob(i)*rcos(i)/(sqr6*fc(i))) + sigm(i)/fc(i)
692 bl(i) = (sigm(i)/fc(i)) + rhob(i)/(sqr6*fc(i))
693 al(i) = (one - qh1(i))*(bl(i)**2) + sqrt(three/two)*(rhob(i)/fc(i))
694 fp(i) = al(i)**2 + m0*(qh1(i)**2)*qh2(i)*cl(i) - (qh1(i)**2)*(qh2(i)**2)
695c
696 ENDDO
697 ! End of the loop over yielding elements
698 ENDDO
699 !End of the loop over the Newton iterations
700 ENDIF
701 !=========================================================================
702 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
703 !=========================================================================
704c
705 !====================================================================
706 ! - UPDATE EQUIVALENT PLASTIC STRAIN
707 !====================================================================
708 DO i=1,nel
709 ! Update cumulated plastic strain for damage and output
710 dpla(i) = sqrt((plaxx(i)-uvar(i,2))**2 + (playy(i)-uvar(i,3))**2 +
711 . (plazz(i)-uvar(i,4))**2 + half*((plaxy(i)-uvar(i,5))**2 +
712 . (playz(i)-uvar(i,6))**2 + (plazx(i)-uvar(i,7))**2))
713 pla(i) = pla(i) + dpla(i)
714 ! Hardening rate
715 IF (dpla(i) > zero) THEN
716 hardp(i) = sqrt((signxx(i)-sigoxx(i))**2 +
717 . (signyy(i)-sigoyy(i))**2 +
718 . (signzz(i)-sigozz(i))**2 +
719 . two*((signxy(i)-sigoxy(i))**2 +
720 . (signyz(i)-sigoyz(i))**2 +
721 . (signzx(i)-sigozx(i))**2))
722 hardp(i) = hardp(i)/dpla(i)
723 ELSE
724 hardp(i) = zero
725 ENDIF
726 ENDDO
727c
728 !====================================================================
729 ! - COMPUTE DAMAGE VARIABLES EVOLUTION
730 !====================================================================
731 nindx3 = 0
732 indx3(1:nel) = 0
733 IF (dflag < 4) THEN
734c
735 ! Compute eigen value and principal stress tensor
736 DO i = 1,nel
737 sigtens(i,1) = signxx(i)
738 sigtens(i,2) = signyy(i)
739 sigtens(i,3) = signzz(i)
740 sigtens(i,4) = signxy(i)
741 sigtens(i,5) = signyz(i)
742 sigtens(i,6) = signzx(i)
743 ENDDO
744 evv(1:mvsiz,1:3) = zero
745 dirprv(1:mvsiz,1:3,1:3) = zero
746 sigpt(1:nel,1:3) = zero
747 sigpc(1:nel,1:3) = zero
748 IF (iresp == 1) THEN
749 CALL valpvecdp_v(sigtens,evv,dirprv,nel)
750 ELSE
751 CALL valpvec_v(sigtens,evv,dirprv,nel)
752 ENDIF
753c
754 ! Loop over elements
755 DO i = 1,nel
756c
757 ! Equivalent strain damage threshold
758 eps0(i) = ft(i)/young
759c
760 ! Compute equivalent strain
761 eps_eq(i) = eps0(i)*m0*half*cl(i)
762 eps_eq(i) = eps_eq(i) + sqrt((eps_eq(i))**2 +
763 . three*half*(eps0(i)*rhob(i)/fc(i))**2)
764c
765 ! Compute ductility measure
766 IF (sigm(i) <= zero) THEN
767 rs = -sqr6*sigm(i)/max(rhob(i),em20)
768 ELSE
769 rs = zero
770 ENDIF
771 xsis(i) = one + (as - one)*(rs**bs)
772c
773 ! Computation of norm of principal stress tensor
774 normsigp(i) = evv(i,1)**2 + evv(i,2)**2 + evv(i,3)**2
775 ENDDO
776c
777 ! Computation of ALPHAC compression factor
778 alphac(1:nel) = zero
779 DO j = 1,3
780 DO i = 1,nel
781 IF (evv(i,j) > zero) THEN
782 sigpt(i,j) = evv(i,j)
783 ELSE
784 sigpc(i,j) = evv(i,j)
785 ENDIF
786 alphac(i) = alphac(i) + sigpc(i,j)*(sigpt(i,j) + sigpc(i,j))/max(normsigp(i),em20)
787 ENDDO
788 ENDDO
789c
790 ! Computation of damage variables
791 DO i = 1,nel
792c
793 ! Damage limit is reached
794 IF (eps_eq(i)>=(eps0(i)-em08)) THEN
795 ! Compute tension damage
796 kapdt2(i) = kapdt2(i) + (max(eps_eq(i) - kapdt(i),zero))/xsis(i)
797 kapdt1(i) = kapdt1(i) + dpla(i)/xsis(i)
798 kapdt(i) = max(eps_eq(i),kapdt(i))
799 eps_inel(i) = kapdt1(i) + omegat(i)*kapdt2(i)
800 ! -> Linear type
801 IF (dtype == 1) THEN
802 omegat(i) = (young*kapdt(i)*wf - ft(i)*wf +
803 . ft(i)*kapdt1(i)*h(i))/(young*kapdt(i)*wf - ft(i)*h(i)*kapdt2(i))
804 ! -> Bilinear type
805 ELSEIF (dtype == 2) THEN
806 IF (h(i)*eps_inel(i) >= zero .AND. h(i)*eps_inel(i) < wf1) THEN
807 omegat(i) = (young*kapdt(i)*wf1 - ft(i)*wf1 - (ft1(i) -
808 . ft(i))*kapdt1(i)*h(i))/(young*kapdt(i)*wf1 +
809 . (ft1(i) - ft(i))*h(i)*kapdt2(i))
810 ELSEIF (h(i)*eps_inel(i) >= wf1 .AND. h(i)*eps_inel(i) < wf) THEN
811 omegat(i) = (young*kapdt(i)*(wf - wf1) -
812 . ft1(i)*(wf - wf1) + ft1(i)*h(i)*kapdt1(i) - ft1(i)*wf1)/
813 . (young*kapdt(i)*(wf - wf1) - ft1(i)*h(i)*kapdt2(i))
814 ELSE
815 omegat(i) = zep999
816 ENDIF
817 ! -> Exponential type
818 ELSEIF (dtype == 3) THEN
819 omegat(i) = one - exp(-(h(i)*eps_inel(i)/wf))
820 ENDIF
821 omegat(i) = max(omegat(i),dmg(i,2))
822 omegat(i) = min(max(omegat(i),zero),zep999)
823 ! Compute compression damage
824 kapdc2(i) = kapdc2(i) + alphac(i)*(max(eps_eq(i) - kapdc(i),zero))/xsis(i)
825 betac = ft(i)*qh2(i)*sqrt(two_third)/
826 . (rhob(i)*sqrt(one + two*(df**2)))
827 kapdc1(i) = kapdc1(i) + alphac(i)*betac*dpla(i)/xsis(i)
828 kapdc(i) = max(alphac(i)*eps_eq(i),kapdc(i))
829 eps_inel(i) = kapdc1(i) + omegac(i)*kapdc2(i)
830 omegac(i) = one - exp(-(eps_inel(i)/efc))
831 omegac(i) = max(omegac(i),dmg(i,3))
832 omegac(i) = min(max(omegac(i),zero),zep999)
833 ENDIF
834c
835 ! Apply damage on stress computation
836 ! -> Standard model with two damage variables
837 IF (dflag == 1) THEN
838 dmg(i,1) = min(min(omegat(i),omegac(i)),zep999)
839 ! Tensile stress tensor
840 sigtxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigpt(i,1)
841 . + dirprv(i,1,2)*dirprv(i,1,2)*sigpt(i,2)
842 . + dirprv(i,1,3)*dirprv(i,1,3)*sigpt(i,3)
843 sigtyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigpt(i,2)
844 . + dirprv(i,2,3)*dirprv(i,2,3)*sigpt(i,3)
845 . + dirprv(i,2,1)*dirprv(i,2,1)*sigpt(i,1)
846 sigtzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigpt(i,3)
847 . + dirprv(i,3,1)*dirprv(i,3,1)*sigpt(i,1)
848 . + dirprv(i,3,2)*dirprv(i,3,2)*sigpt(i,2)
849 sigtxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigpt(i,1)
850 . + dirprv(i,1,2)*dirprv(i,2,2)*sigpt(i,2)
851 . + dirprv(i,1,3)*dirprv(i,2,3)*sigpt(i,3)
852 sigtyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigpt(i,2)
853 . + dirprv(i,2,3)*dirprv(i,3,3)*sigpt(i,3)
854 . + dirprv(i,2,1)*dirprv(i,3,1)*sigpt(i,1)
855 sigtzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigpt(i,3)
856 . + dirprv(i,3,1)*dirprv(i,1,1)*sigpt(i,1)
857 . + dirprv(i,3,2)*dirprv(i,1,2)*sigpt(i,2)
858 ! Compressive stress tensor
859 sigcxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigpc(i,1)
860 . + dirprv(i,1,2)*dirprv(i,1,2)*sigpc(i,2)
861 . + dirprv(i,1,3)*dirprv(i,1,3)*sigpc(i,3)
862 sigcyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigpc(i,2)
863 . + dirprv(i,2,3)*dirprv(i,2,3)*sigpc(i,3)
864 . + dirprv(i,2,1)*dirprv(i,2,1)*sigpc(i,1)
865 sigczz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigpc(i,3)
866 . + dirprv(i,3,1)*dirprv(i,3,1)*sigpc(i,1)
867 . + dirprv(i,3,2)*dirprv(i,3,2)*sigpc(i,2)
868 sigcxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigpc(i,1)
869 . + dirprv(i,1,2)*dirprv(i,2,2)*sigpc(i,2)
870 . + dirprv(i,1,3)*dirprv(i,2,3)*sigpc(i,3)
871 sigcyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigpc(i,2)
872 . + dirprv(i,2,3)*dirprv(i,3,3)*sigpc(i,3)
873 . + dirprv(i,2,1)*dirprv(i,3,1)*sigpc(i,1)
874 sigczx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigpc(i,3)
875 . + dirprv(i,3,1)*dirprv(i,1,1)*sigpc(i,1)
876 . + dirprv(i,3,2)*dirprv(i,1,2)*sigpc(i,2)
877 ! New stress tensor
878 signxx(i) = (one-omegat(i))*sigtxx(i) + (one-omegac(i))*sigcxx(i)
879 signyy(i) = (one-omegat(i))*sigtyy(i) + (one-omegac(i))*sigcyy(i)
880 signzz(i) = (one-omegat(i))*sigtzz(i) + (one-omegac(i))*sigczz(i)
881 signxy(i) = (one-omegat(i))*sigtxy(i) + (one-omegac(i))*sigcxy(i)
882 signyz(i) = (one-omegat(i))*sigtyz(i) + (one-omegac(i))*sigcyz(i)
883 signzx(i) = (one-omegat(i))*sigtzx(i) + (one-omegac(i))*sigczx(i)
884 ! -> Isotropic damage model with one damage variable
885 ELSEIF (dflag == 2) THEN
886 dmg(i,1) = min(omegat(i),zep999)
887 signxx(i) = (one-dmg(i,1))*signxx(i)
888 signyy(i) = (one-dmg(i,1))*signyy(i)
889 signzz(i) = (one-dmg(i,1))*signzz(i)
890 signxy(i) = (one-dmg(i,1))*signxy(i)
891 signyz(i) = (one-dmg(i,1))*signyz(i)
892 signzx(i) = (one-dmg(i,1))*signzx(i)
893 ! -> Multiplicative damage model
894 ELSEIF (dflag == 3) THEN
895 dmg(i,1) = min(one - (one-omegat(i))*(one-omegac(i)),zep999)
896 signxx(i) = (one-dmg(i,1))*signxx(i)
897 signyy(i) = (one-dmg(i,1))*signyy(i)
898 signzz(i) = (one-dmg(i,1))*signzz(i)
899 signxy(i) = (one-dmg(i,1))*signxy(i)
900 signyz(i) = (one-dmg(i,1))*signyz(i)
901 signzx(i) = (one-dmg(i,1))*signzx(i)
902 ENDIF
903c
904 ! Element deletion check
905 IF ((dmg(i,1) == zep999).AND.(loff(i) == one).AND.(off(i) > zero)) THEN
906 signxx(i) = zero
907 signyy(i) = zero
908 signzz(i) = zero
909 signxy(i) = zero
910 signyz(i) = zero
911 signzx(i) = zero
912 dmg(i,1) = one
913 loff(i) = four_over_5
914 IF (idel > 1) THEN
915 nindx3 = nindx3 + 1
916 indx3(nindx3) = i
917 off(i) = zero
918 ENDIF
919 ENDIF
920c
921 ENDDO
922 ENDIF
923c
924 ! Storing new values
925 DO i=1,nel
926 ! Update user variable
927 uvar(i,1) = kap(i)
928 uvar(i,2) = plaxx(i)
929 uvar(i,3) = playy(i)
930 uvar(i,4) = plazz(i)
931 uvar(i,5) = plaxy(i)
932 uvar(i,6) = playz(i)
933 uvar(i,7) = plazx(i)
934 uvar(i,8) = eps_eq(i)
935 uvar(i,9) = kapdt(i)
936 uvar(i,10) = kapdt1(i)
937 uvar(i,11) = kapdt2(i)
938 uvar(i,12) = kapdc(i)
939 uvar(i,13) = kapdc1(i)
940 uvar(i,14) = kapdc2(i)
941 ! Damage variables
942 dmg(i,2) = omegat(i)
943 dmg(i,3) = omegac(i)
944 ! Coefficient for hourglass
945 IF (dpla(i) > zero) THEN
946 et(i) = hardp(i) / (hardp(i) + young)
947 ELSE
948 et(i) = one
949 ENDIF
950 ! Computation of the sound speed
951 soundsp(i) = sqrt((bulk + four_over_3*g)/rho(i))
952 ENDDO
953c
954 IF (irate > 1) THEN
955#include "vectorize.inc"
956 DO ii = 1,nindx
957 i = indx(ii)
958 IF (uvar(i,15) == zero) uvar(i,15) = arate(i)
959 ENDDO
960#include "vectorize.inc"
961 DO ii = 1,nindx2
962 i = indx2(ii)
963 IF (uvar(i,15) == zero) uvar(i,15) = arate(i)
964 ENDDO
965 ENDIF
966
967 ! Printout failure
968 IF (nindx3>0) THEN
969 DO i=1,nindx3
970#include "lockon.inc"
971 WRITE(iout, 2000) ngl(indx3(i))
972 WRITE(istdo,2100) ngl(indx3(i)),tt
973#include "lockoff.inc"
974 ENDDO
975 ENDIF
976c
977 2000 FORMAT(1x,'FAILURE (CDPM2) OF SOLID ELEMENT ',i10)
978 2100 FORMAT(1x,'FAILURE (CDPM2) OF SOLID ELEMENT ',i10,1x,'AT TIME :',1pe12.4)
979c
980 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps124
subroutine sigeps124(nel, ngl, nuparam, nuvar, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, loff, ipg, npg, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, dmg, le)
Definition sigeps124.F:39
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2231
sigeps33
subroutine sigeps33(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off)
Definition sigeps33.F:46
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1902