OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps78c.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!|| sigeps78c ../engine/source/materials/mat/mat078/sigeps78c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| vinter ../engine/source/tools/curve/vinter.F
29!||====================================================================
30 SUBROUTINE sigeps78c(
31 1 NEL ,NUPARAM ,NUVAR ,NVARTMP ,TIME ,
32 2 NFUNC ,IFUNC ,NPF ,TF ,UPARAM ,
33 3 THKLY ,THK ,GS ,ETSE ,SIGY ,
34 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
35 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
36 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
37 7 SOUNDSP ,UVAR ,SIGA ,SIGB ,SIGC ,
38 8 RHO ,OFF ,PLA ,DEP ,VARTMP ,
39 9 INLOC ,DPLANL ,ETIMP ,SEQ ,LOFF )
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 "com01_c.inc"
48C-----------------------------------------------
49C I N P U T A r g u m e n t s
50C-----------------------------------------------
51 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,
52 . NVARTMP,INLOC
53 my_real ,INTENT(IN) ::
54 . TIME
55 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) ::
56 . UPARAM
57 my_real ,DIMENSION(NEL) ,INTENT(IN) ::
58 . RHO,THKLY,OFF,GS,DPLANL,
59 . DEPSXX,DEPSYY,DEPSXY,DEPSYZ,DEPSZX,
60 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
61 INTEGER ,INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
62 my_real, DIMENSION(NEL), INTENT(IN) :: loff
63C-----------------------------------------------
64C I N P U T O U T P U T A r g u m e n t s
65C-----------------------------------------------
66 my_real ,DIMENSION(NEL,3) ,INTENT(INOUT) ::
67 . siga,sigb
68 my_real ,DIMENSION(NEL,4) ,INTENT(INOUT) ::
69 . sigc
70 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) ::
71 . uvar
72 my_real ,DIMENSION(NEL) ,INTENT(INOUT) ::
73 . pla,thk
74 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
75 . soundsp,dep,etse,sigy,seq,
76 . signxx,signyy,signxy,signyz,signzx,etimp
77C-----------------------------------------------
78C VARIABLES FOR FUNCTION INTERPOLATION
79C-----------------------------------------------
80 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
81 my_real
82 . tf(*)
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER :: I,J,K,OPTE,OPTR,NITER,NINDX,IPLAS,ITER
87 INTEGER ,DIMENSION(NEL) :: INDEX,IPOS,ILEN,IAD2
88 my_real ::
89 . nu,bsat,einf,coe,fct,gsigma,dgsigma,rdot,
90 . eini,yield,byu,myu,hyu,rsat,p1,p2,p3,p4,n3,cst,
91 . cstt,df_dk1,df_dk2,dk1_dsigxx,dk1_dsigyy,dk2_dsigxx,
92 . dk2_dsigyy,dk2_dsigxy,normxx,normyy,normxy,dfdsig_c_2,
93 . dfdsig_a,dfdsig_alpha,dfdsig_beta,dphi_dlam,dpla_dlam,
94 . drnih,dmu,c1_kh,c2_kh,norm_sp,t1,dvxx_dlam,dvyy_dlam,dvxy_dlam,
95 . daxx_dlam,dayy_dlam,daxy_dlam,du_dlam
96 REAL(KIND=8) :: dlam,dep_dp(nel)
97 my_real ,DIMENSION(NEL) ::
98 . young,shear,ayu,scalee,h,a1,a2,axx,ayy,axy,k1,k2,cyu,asta,camb,
99 . seff,r,rnih,depszz,deelzz,devboxx,devboyy,devbozz,devboxy,dydxe,
100 . deplxx,deplyy,deplxy,deplzz,devbxx,devbyy,devbzz,devbxy,phi,
101 . dbxx,dbyy,dbzz,dbxy,bqxx,bqyy,bqzz,bqxy,yld,max_asta,canbnxx,canbnyy,canbnxy,
102 . cannxx,cannyy,cannxy,ddep,u,vxx,vyy,vxy,ka1,ka2
103 my_real ,DIMENSION(NEL,3) :: sigbo
104C=======================================================================
105C SIGA (alpha) = CENTER OF YIELD SURFACE ALPHA (Backstress)
106C SIGB (beta) = CENTER OF BOUNDING SURFACE
107C SIGC (q) = CENTER OF NON-IH SURFACE G_SIGMA
108c
109C UVAR(I,1) = R : ISOTROPIC HARDENING OF BOUNDING SURFACE
110C UVAR(I,2) = r : RADIUS OF G_SIGMA (RNIH)
111C UVAR(I,3) = a = B + R - YIELD (AYU)
112C UVAR(I,4) = Yld Stress
113C=======================================================================
114 norm_sp = -huge(norm_sp)
115 normxx = -huge(normxx)
116 normyy = -huge(normyy)
117 normxy = -huge(normxy)
118 !=======================================================================
119 ! - INITIALISATION OF COMPUTATION ON TIME STEP
120 !=======================================================================
121 ! Recovering model parameters
122 eini = uparam(1) ! Initial Young Modulus
123 nu = uparam(2) ! Poisson's ratio
124 yield = uparam(3) ! Yield stress
125 byu = uparam(4) ! Center of the bounding surface
126 c2_kh = uparam(5) ! Parameter for kinematic hardening rule of yield surface (Lower case if C1_KH is defined)
127 hyu = uparam(6) ! Material parameter for controlling work hardening stagnation
128 bsat = uparam(7) ! Initial size of the bounding surface
129 myu = uparam(8) ! Parameter for isotropic and kinematic hardening of the bounding surface
130 rsat = uparam(9) ! Saturated value of the isotropic hardening stress
131 einf = uparam(10) ! Asymptotic value of Young's modulus
132 coe = uparam(11) ! Parameter controlling the dependency of Young's modulus on the effective plastic strain
133 opte = uparam(12) ! Modified Young's modulus
134 optr = uparam(13) ! Modified isotropic hardening rule flag
135 p1 = uparam(14) ! First yield criterion parameter
136 p2 = uparam(15) ! First yield criterion parameter
137 p3 = uparam(16) ! First yield criterion parameter
138 p4 = uparam(17) ! First yield criterion parameter
139 n3 = uparam(18) ! Barlat 1989 exponent
140 cst = uparam(19) ! Constant used in the modified formulation of the isotropic hardening of bounding surface
141 cstt = uparam(20) !
142 iplas = nint(uparam(21)) ! Flag for choosing yield criterion (Hill 1948 - Barlat 1989)
143 c1_kh = uparam(22) ! Upper case for controlling work hardening stagnation
144c
145 IF (iplas == 1) THEN
146 norm_sp = one
147 ELSEIF(iplas == 2) THEN
148 norm_sp = ep03 / eini
149 ENDIF
150 ! Initial value
151 IF (isigi == 0 .and. time == zero) THEN
152 DO i=1,nel
153 ! Initial AYU value
154 uvar(i,3) = bsat - yield
155 ! Initial yield stress
156 uvar(i,4) = yield
157 ENDDO
158 ENDIF
159C
160 ! Update of Young modulus
161 ! -> Constant young modulus
162 IF (coe == zero .and. ifunc(1) == 0) THEN
163 young(1:nel) = eini
164 ! -> Tabulated young modulus
165 ELSE IF (opte == 1) THEN
166 ipos(1:nel) = vartmp(1:nel,1)
167 iad2(1:nel) = npf(ifunc(1)) / 2 + 1
168 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad2(1:nel) - ipos(1:nel)
169 CALL vinter(tf,iad2,ipos,ilen,nel,pla,dydxe,scalee)
170 young(1:nel) = eini * scalee(1:nel)
171 vartmp(1:nel,1) = ipos(1:nel)
172 ! -> Non-linear young modulus
173 ELSE
174 DO i=1,nel
175 IF (pla(i) > zero)THEN
176 young(i) = eini-(eini-einf)*(one-exp(-coe*pla(i)))
177 ELSE
178 young(i) = eini
179 ENDIF
180 ENDDO
181 ENDIF
182C
183 ! Recovering elastic parameters and hardening variables
184 DO i=1,nel
185 a1(i) = young(i)/(one-nu*nu) ! Plane stress elastic matrix diagonal component
186 a2(i) = nu*young(i)/(one-nu*nu) ! Plane stress elastic matrix non-diagonal component
187 shear(i) = young(i)*half/(one+nu) ! Shear modulus
188 r(i) = uvar(i,1) ! Isotropic hardening
189 rnih(i) = uvar(i,2) ! Non-IH hardening
190 ayu(i) = uvar(i,3) ! AYU variable
191 max_asta(i) = uvar(i,6) ! Maximal value of NORM(ALPHA_STAR)
192 dep_dp(i) = zero
193 ddep(i) = zero
194 deplxx(i) = zero
195 deplyy(i) = zero
196 deplxy(i) = zero
197 deplzz(i) = zero ! Plastic strain increment for thickness variation
198 etimp(i) = young(i)/eini ! For implicit
199 sigbo(i,1) = sigb(i,1)
200 sigbo(i,2) = sigb(i,2)
201 sigbo(i,3) = sigb(i,3)
202 etse(i) = one ! For hourglassing
203 u(i) = one
204 ENDDO
205c
206 !========================================================================
207 ! - COMPUTATION OF TRIAL VALUES
208 !========================================================================
209 DO i=1,nel
210 ! Trial stress tensor
211 signxx(i) = sigoxx(i) + a1(i)*depsxx(i)+a2(i)*depsyy(i)
212 signyy(i) = sigoyy(i) + a2(i)*depsxx(i)+a1(i)*depsyy(i)
213 signxy(i) = sigoxy(i) + shear(i) * depsxy(i)
214 signyz(i) = sigoyz(i) + gs(i) * depsyz(i)
215 signzx(i) = sigozx(i) + gs(i) * depszx(i)
216 ! Update with the back stress (Backstress ALPHA = SIGA + SIGB)
217 axx(i) = signxx(i) - siga(i,1) - sigb(i,1)
218 ayy(i) = signyy(i) - siga(i,2) - sigb(i,2)
219 axy(i) = signxy(i) - siga(i,3) - sigb(i,3)
220 ! Choice of yield criterion
221 ! -> Hill 1948
222 IF (iplas == 1) THEN
223 seff(i) = sqrt(axx(i)**2 - two*p2*axx(i)*ayy(i) + p1*ayy(i)**2 + p3*axy(i)**2)
224 !effective back stress alpha*
225 asta(i) = sqrt(siga(i,1)**2 - two*p2*siga(i,1)*siga(i,2) + p1*siga(i,2)**2 + p3*siga(i,3)**2)
226 ! -> Barlat 1989
227 ELSEIF (iplas == 2) THEN
228 k1(i) = (axx(i) + p3*ayy(i))/two
229 k2(i) = (sqrt(((axx(i) - p3*ayy(i))/two)**two + (p4*axy(i))**two))
230 seff(i) = p1*(abs(k1(i)+k2(i))*norm_sp)**n3 + p1*(abs(k1(i)-k2(i))*norm_sp)**n3 + p2*(abs(two*k2(i))*norm_sp)**n3
231 seff(i) = (half*max(seff(i),em20))**(one/n3)
232 !effective back stress alpha*
233 ka1(i) = (siga(i,1) + p3*siga(i,2))/two
234 ka2(i) = sqrt(((siga(i,1) - p3*siga(i,2))/two)**two + (p4*siga(i,3))**two)
235 asta(i) = p1*abs(ka1(i)+ka2(i))**n3 + p1*abs(ka1(i)-ka2(i))**n3 + p2*abs(two*ka2(i))**n3
236 asta(i) = (half*max(em20,asta(i)))**(one/n3)
237 ENDIF
238
239 asta(i) = max(asta(i),em20)
240 max_asta(i) = max(max_asta(i),asta(i))
241 IF (max_asta(i) < bsat - yield) THEN
242 cyu(i) = c1_kh
243 ELSE
244 cyu(i) = c2_kh
245 ENDIF
246 camb(i) = (cyu(i)*ayu(i) + myu*byu)/yield
247 ! component used in sig - alpha:
248 t1 = cyu(i)*sqrt(ayu(i)/asta(i))
249 canbnxx(i) = t1 * siga(i,1) + myu*sigb(i,1)
250 canbnyy(i) = t1 * siga(i,2) + myu*sigb(i,2)
251 canbnxy(i) = t1 * siga(i,3) + myu*sigb(i,3)
252
253 cannxx(i) = t1 * siga(i,1)
254 cannyy(i) = t1 * siga(i,2)
255 cannxy(i) = t1 * siga(i,3)
256 !------------------------------------------
257 ! yield function value
258 phi(i) = seff(i) /norm_sp - yield
259 !------------------------------------------
260 ENDDO
261c
262 !========================================================================
263 ! - CHECKING PLASTIC BEHAVIOR
264 !========================================================================
265 nindx = 0
266 DO i=1,nel
267 IF (phi(i) >= zero .AND. off(i) == one) THEN
268 nindx = nindx+1
269 index(nindx) = i
270 ENDIF
271 ENDDO
272#include "vectorize.inc"
273 DO j=1,nindx
274 i = index(j)
275 vxx(i) = axx(i)
276 vyy(i) = ayy(i)
277 vxy(i) = axy(i)
278 ENDDO
279
280c
281 !====================================================================
282 ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
283 !====================================================================
284c
285 ! Number of Newton iterations
286 niter = 3
287c
288 DO iter = 1,niter
289#include "vectorize.inc"
290 ! Loop over yielding elements
291 DO j=1,nindx
292 ! Number of the element with plastic behaviour
293 i = index(j)
294c
295 ! Note : in this part, the purpose is to compute for each iteration
296 ! a plastic multiplier allowing to update internal variables to satisfy
297 ! the consistency condition using the cutting plane method
298 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
299 ! -> PHI : current value of yield function (known)
300 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
301 ! into account of internal variables kinetic :
302 ! hardening parameters
303C
304 !---------------------------------------------------------
305 ! 1 - Computation of the normal to the yield surface
306 !---------------------------------------------------------
307 ! -> Hill 1948
308 IF (iplas == 1) THEN
309 normxx = (axx(i)-p2*ayy(i))/(max(seff(i),em20))
310 normyy = (p1*ayy(i)-p2*axx(i))/(max(seff(i),em20))
311 normxy = (p3*axy(i))/(max(seff(i),em20))
312 ! -> Barlat 1989
313 ELSEIF (iplas == 2) THEN
314 ! Computation of the derivatives
315 df_dk1 = ((seff(i)/norm_sp)**(1-n3))*(p1/two)* (
316 . + sign(one,k1(i)+k2(i))*(abs(k1(i)+k2(i))**(n3-1))
317 . + sign(one,k1(i)-k2(i))*(abs(k1(i)-k2(i))**(n3-1)))
318 df_dk2 = ((seff(i)/norm_sp)**(1-n3))*((p1/two)*(
319 . + sign(one,k1(i)+k2(i))*(abs(k1(i)+k2(i))**(n3-1))
320 . - sign(one,k1(i)-k2(i))*(abs(k1(i)-k2(i))**(n3-1)))
321 . + p2*(abs(two*k2(i))**(n3-1)))
322 dk1_dsigxx = half
323 dk1_dsigyy = p3 /two
324 dk2_dsigxx = (axx(i)-p3*ayy(i)) /(max(four*k2(i),em20))
325 dk2_dsigyy = -p3*(axx(i)-p3*ayy(i)) /(max(four*k2(i),em20))
326 dk2_dsigxy = (p4**two)*axy(i) /max(k2(i),em20)
327 ! Assembling the normal
328 normxx = df_dk1*dk1_dsigxx+ df_dk2*dk2_dsigxx
329 normyy = df_dk1*dk1_dsigyy+ df_dk2*dk2_dsigyy
330 normxy = df_dk2*dk2_dsigxy
331 ENDIF
332c
333
334 ! 2 - Computation of DPHI_DLAMBDA
335 !---------------------------------------------------------
336
337
338 ! Computation of the tensor product Norm * (SigN - Alpha)
339 dfdsig_a = normxx * axx(i)
340 . + normyy * ayy(i)
341 . + normxy * axy(i)
342
343 dpla_dlam = dfdsig_a/yield
344
345 du_dlam = -u(i)**2 * camb(i) * dpla_dlam
346 dvxx_dlam = -(a1(i)*normxx + a2(i)*normyy) + canbnxx(i) * dpla_dlam
347 dvyy_dlam = -(a1(i)*normyy + a2(i)*normxx) + canbnyy(i) * dpla_dlam
348 dvxy_dlam = -shear(i) *normxy + canbnxy(i) * dpla_dlam
349
350 daxx_dlam = du_dlam * vxx(i) + dvxx_dlam * u(i)
351 dayy_dlam = du_dlam * vyy(i) + dvyy_dlam * u(i)
352 daxy_dlam = du_dlam * vxy(i) + dvxy_dlam * u(i)
353
354 ! Assembling the derivative dphi/dlam = dphi/dA * dA/dlam
355
356 dphi_dlam = normxx * daxx_dlam + normyy * dayy_dlam + normxy * daxy_dlam
357 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
358
359 ! 3 - Computation of plastic multiplier and variables update
360 !----------------------------------------------------------
361c
362 ! Plastic multiplier
363 dlam = - phi(i)/dphi_dlam
364c
365 ! Plastic strains tensor update
366 deplxx(i) = dlam*normxx
367 deplyy(i) = dlam*normyy
368 deplxy(i) = dlam*normxy
369c
370 ! Cumulated plastic strain update
371 pla(i) = pla(i) + dlam*dpla_dlam
372c
373 ddep(i) = dlam * dpla_dlam
374 ! Total plastic strain increment on the time step
375 dep_dp(i) = dep_dp(i) + dlam*dpla_dlam
376 dep_dp(i) = max(dep_dp(i),zero)
377c
378 ! Cauchy stress tensor update
379 signxx(i) = signxx(i) - a1(i)*deplxx(i) - a2(i)*deplyy(i)
380 signyy(i) = signyy(i) - a2(i)*deplxx(i) - a1(i)*deplyy(i)
381 signxy(i) = signxy(i) - shear(i)*deplxy(i)
382
383 !update a_ij = SigN - Alpha
384 u(i) = one /(one + camb(i) * ddep(i) )
385 vxx(i) = signxx(i) - siga(i,1) - sigb(i,1) + canbnxx(i) * ddep(i)
386 vyy(i) = signyy(i) - siga(i,2) - sigb(i,2) + canbnyy(i) * ddep(i)
387 vxy(i) = signxy(i) - siga(i,3) - sigb(i,3) + canbnxy(i) * ddep(i)
388 axx(i) = u(i) * vxx(i)
389 ayy(i) = u(i) * vyy(i)
390 axy(i) = u(i) * vxy(i)
391
392c
393 ! Alpha* back stress tensor update
394 siga(i,1) = siga(i,1) + ((cyu(i)*ayu(i)*axx(i)/yield) - cannxx(i) ) * ddep(i)
395 siga(i,2) = siga(i,2) + ((cyu(i)*ayu(i)*ayy(i)/yield) - cannyy(i) ) * ddep(i)
396 siga(i,3) = siga(i,3) + ((cyu(i)*ayu(i)*axy(i)/yield) - cannxy(i) ) * ddep(i)
397c
398 ! Beta back stress tensor update
399 sigb(i,1) = sigb(i,1) + ((myu *byu *axx(i)/yield) - myu*sigbo(i,1))* ddep(i)
400 sigb(i,2) = sigb(i,2) + ((myu *byu *ayy(i)/yield) - myu*sigbo(i,2))* ddep(i)
401 sigb(i,3) = sigb(i,3) + ((myu *byu *axy(i)/yield) - myu*sigbo(i,3))* ddep(i)
402
403 ! New equivalent stress
404 ! -> Hill 1948
405 IF (iplas == 1) THEN
406 seff(i) = sqrt(axx(i)**2 - two*p2*axx(i)*ayy(i) + p1*ayy(i)**2 + p3*axy(i)**2)
407 ! -> Barlat 1989
408 ELSEIF (iplas == 2) THEN
409 k1(i) = (axx(i) + p3*ayy(i))/two
410 k2(i) = sqrt(((axx(i) - p3*ayy(i))/two)**two + (p4*axy(i))**two)
411 seff(i) = p1*(abs(k1(i)+k2(i))*norm_sp)**n3 + p1*(abs(k1(i)-k2(i))*norm_sp)**n3 + p2*(abs(two*k2(i))*norm_sp)**n3
412 seff(i) = (half*seff(i))**(one/n3)
413 ENDIF
414c
415 ! New yield function value
416 phi(i) = seff(i)/norm_sp - yield
417c
418 ! Thickness plastic strain update
419 IF (inloc == 0) deplzz(i) = deplzz(i) - (deplxx(i)+deplyy(i))
420 ENDDO
421 ENDDO
422 !===================================================================
423 ! - END OF PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
424 !===================================================================
425#include "vectorize.inc"
426 ! Loop over yielding elements
427 DO j=1,nindx
428 ! Number of the element with plastic behaviour
429 i = index(j)
430
431 IF (iplas == 1) THEN
432 asta(i) = sqrt(siga(i,1)**2 - two*p2*siga(i,1)*siga(i,2) + p1*siga(i,2)**2 + p3*siga(i,3)**2)
433 ! -> Barlat 1989
434 ELSEIF (iplas == 2) THEN
435 ka1(i) = (siga(i,1) + p3*siga(i,2))/two
436 ka2(i) = sqrt(((siga(i,1) - p3*siga(i,2))/two)**two + (p4*siga(i,3))**two)
437 asta(i) = p1*abs(ka1(i)+ka2(i))**n3 + p1*abs(ka1(i)-ka2(i))**n3 + p2*abs(two*ka2(i))**n3
438 asta(i) = (half*max(em20,asta(i)))**(one/n3)
439 asta(i) = asta(i)
440 ENDIF
441 asta(i) = max(asta(i),em20)
442 max_asta(i) = max(max_asta(i),asta(i))
443 IF (max_asta(i) < bsat - yield) THEN
444 cyu(i) = c1_kh
445 ELSE
446 cyu(i) = c2_kh
447 ENDIF
448
449 ! Old beta deviator
450 devboxx(i) = two_third*sigbo(i,1) - third*sigbo(i,2)
451 devboyy(i) = two_third*sigbo(i,2) - third*sigbo(i,1)
452 devbozz(i) = -third*sigbo(i,1) - third*sigbo(i,2)
453 devboxy(i) = sigbo(i,3)
454c
455 ! New beta deviator
456 devbxx(i) = two_third*sigb(i,1) - third*sigb(i,2)
457 devbyy(i) = two_third*sigb(i,2) - third*sigb(i,1)
458 devbzz(i) = -third*sigb(i,1) - third*sigb(i,2)
459 devbxy(i) = sigb(i,3)
460c
461 !-------------------------------------------------------------------------------------
462 ! Updating the hardening
463 !-------------------------------------------------------------------------------------
464 ! -> Computation of the equivalent stress for the bounding surface Gsigma(DevBeta - q_n)
465 bqxx(i) = devbxx(i) - sigc(i,1)
466 bqyy(i) = devbyy(i) - sigc(i,2)
467 bqzz(i) = devbzz(i) - sigc(i,3)
468 bqxy(i) = devbxy(i) - sigc(i,4)
469 ! -> Computation of Beta increment
470 dbxx(i) = devbxx(i) - devboxx(i)
471 dbyy(i) = devbyy(i) - devboyy(i)
472 dbzz(i) = devbzz(i) - devbozz(i)
473 dbxy(i) = devbxy(i) - devboxy(i)
474 ! Computation of GSIGMA J2-type yield function
475 gsigma = three_half*(bqxx(i)*bqxx(i) + bqyy(i)*bqyy(i) + bqzz(i)*bqzz(i) + two*bqxy(i)*bqxy(i))
476 . - rnih(i)*rnih(i)
477 ! Computation of (Beta - q)*DBeta
478 dgsigma = bqxx(i)*dbxx(i) + bqyy(i)*dbyy(i) + bqzz(i)*dbzz(i) + two*bqxy(i)*dbxy(i)
479 ! If the GSIGMA yield function is greater that 0 -> SIGC must be updated
480 IF ((gsigma >= zero).AND.(dgsigma > zero)) THEN
481 ! Isotropic hardening
482 IF (optr == one) THEN
483 rdot = rsat*((pla(i)+cst)**cstt-cst**cstt) - r(i)
484 r(i) = rsat*((pla(i)+cst)**cstt-cst**cstt)
485 ELSE
486 rdot = myu*(rsat - r(i))*dep_dp(i)
487 r(i) = r(i) + rdot
488 ENDIF
489 ayu(i) = bsat + r(i) - yield
490 ! Update of q back stress tensor in case of work-hardening stagnation
491 IF (hyu > zero) THEN
492 ! Computation of DMU
493 IF (rnih(i) == zero) THEN
494 dmu = (gsigma/(three*hyu*dgsigma)) - one
495 ELSE
496 dmu = (one-hyu)*three_half*dgsigma/(rnih(i)**2)
497 ENDIF
498 ! q back stress update
499 sigc(i,1) = devbxx(i) - bqxx(i)/(one+dmu)
500 sigc(i,2) = devbyy(i) - bqyy(i)/(one+dmu)
501 sigc(i,3) = devbzz(i) - bqzz(i)/(one+dmu)
502 sigc(i,4) = devbxy(i) - bqxy(i)/(one+dmu)
503 ! -> Update of the equivalent stress for the bounding surface Gsigma(DevBeta - q_n)
504 bqxx(i) = devbxx(i) - sigc(i,1)
505 bqyy(i) = devbyy(i) - sigc(i,2)
506 bqzz(i) = devbzz(i) - sigc(i,3)
507 bqxy(i) = devbxy(i) - sigc(i,4)
508 ! Update of (Beta - q)*DBeta
509 dgsigma = bqxx(i)*dbxx(i) + bqyy(i)*dbyy(i) + bqzz(i)*dbzz(i) + two*bqxy(i)*dbxy(i)
510 ! Computation of RNIH evolution
511 IF (rdot > zero) THEN
512 IF (rnih(i) == zero) THEN
513 ! Radius initial value
514 rnih(i) = three*hyu*dgsigma
515 ELSE
516 ! Radius increment
517 drnih = hyu*three_half*dgsigma/rnih(i)
518 ! Radius update
519 rnih(i) = rnih(i) + drnih
520 ENDIF
521 ENDIF
522 ENDIF
523 ENDIF
524
525 ENDDO
526C
527 !============================================================
528 ! - STORING NEW VALUES
529 !============================================================
530#include "vectorize.inc"
531 DO j = 1,nindx
532 ! Number of the yielding element
533 i = index(j)
534 ! Equivalent stress
535 ! -> Hill 1948
536 IF (iplas == 1) THEN
537 yld(i) = sqrt(signxx(i)**2 - two*p2*signxx(i)*signyy(i) + p1*signyy(i)**2 + p3*signxy(i)**2)
538 ! -> Barlat 1989
539 ELSEIF (iplas == 2) THEN
540 k1(i) = (signxx(i) + p3*signyy(i))/two
541 k2(i) = sqrt(((signxx(i) - p3*signyy(i))/two)**two + (p4*signxy(i))**two)
542 yld(i) = p1*abs(k1(i)+k2(i))**n3 + p1*abs(k1(i)-k2(i))**n3 + p2*abs(two*k2(i))**n3
543 yld(i) = (half*max(yld(i),em20))**(one/n3)
544 yld(i) = yld(i)
545 ENDIF
546 ! Hardening rate
547 IF (dep_dp(i) > zero) THEN
548 h(i) = abs(yld(i)-uvar(i,4))/dep_dp(i)
549 !H(I) = MAX(EM20,H(I))
550 etse(i) = h(i) / (h(i)+young(i))
551 ENDIF
552 ! Storing new values
553 uvar(i,1) = r(i)
554 uvar(i,2) = rnih(i)
555 uvar(i,3) = ayu(i)
556 uvar(i,4) = yld(i)
557 uvar(i,6) = max_asta(i)
558 END DO
559 ! Equivalent stress output, thickness update and soundspeed
560 DO i=1,nel
561 ! New equivalent stress
562 axx(i) = signxx(i) - siga(i,1) - sigb(i,1)
563 ayy(i) = signyy(i) - siga(i,2) - sigb(i,2)
564 axy(i) = signxy(i) - siga(i,3) - sigb(i,3)
565 ! -> Hill 1948
566 IF (iplas == 1) THEN
567 seff(i) = sqrt(axx(i)**2 - two*p2*axx(i)*ayy(i) + p1*ayy(i)**2 + p3*axy(i)**2)
568 ! -> Barlat 1989
569 ELSEIF (iplas == 2) THEN
570 k1(i) = (axx(i) + p3*ayy(i))/two
571 k2(i) = sqrt(((axx(i) - p3*ayy(i))/two)**two + (p4*axy(i))**two)
572 seff(i) = p1*abs(k1(i)+k2(i))**n3 + p1*abs(k1(i)-k2(i))**n3 + p2*abs(two*k2(i))**n3
573 seff(i) = (half*seff(i))**(one/n3)
574 ENDIF
575 seq(i) = seff(i)
576 deelzz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i)
577
578
579
580 !NON LOCAL
581 IF ((inloc > 0).AND.(loff(i) == one)) THEN
582 ! -> Hill 1948
583 IF (iplas == 1) THEN
584 seff(i) = sqrt(axx(i)**2 - two*p2*axx(i)*ayy(i) + p1*ayy(i)**2 + p3*axy(i)**2)
585 normxx = (axx(i)-p2*ayy(i))/(max(seff(i),em20))
586 normyy = (p1*ayy(i)-p2*axx(i))/(max(seff(i),em20))
587 normxy = (p3*axy(i))/(max(seff(i),em20))
588 ! -> Barlat 1989
589 ELSEIF (iplas == 2) THEN
590 k1(i) = (axx(i) + p3*ayy(i))/two
591 k2(i) = (sqrt(((axx(i) - p3*ayy(i))/two)**two + (p4*axy(i))**two))
592 seff(i) = p1*abs(k1(i)+k2(i))**n3 + p1*abs(k1(i)-k2(i))**n3 + p2*abs(two*k2(i))**n3
593 seff(i) = (half*max(seff(i),em20))**(one/n3)
594 df_dk1 = (seff(i)**(1-n3))*(p1/two)*(
595 . + sign(one,k1(i)+k2(i))*(abs(k1(i)+k2(i))**(n3-1))
596 . + sign(one,k1(i)-k2(i))*(abs(k1(i)-k2(i))**(n3-1)))
597 df_dk2 = (seff(i)**(1-n3))*((p1/two)*(
598 . + sign(one,k1(i)+k2(i))*(abs(k1(i)+k2(i))**(n3-1))
599 . - sign(one,k1(i)-k2(i))*(abs(k1(i)-k2(i))**(n3-1)))
600 . + p2*(abs(two*k2(i))**(n3-1)))
601 dk1_dsigxx = half
602 dk1_dsigyy = p3/two
603 dk2_dsigxx = (axx(i)-p3*ayy(i))/(max(four*k2(i),em20))
604 dk2_dsigyy = -p3*(axx(i)-p3*ayy(i))/(max(four*k2(i),em20))
605 dk2_dsigxy = (p4**two)*axy(i)/max(k2(i),em20)
606 normxx = df_dk1*dk1_dsigxx + df_dk2*dk2_dsigxx
607 normyy = df_dk1*dk1_dsigyy + df_dk2*dk2_dsigyy
608 normxy = df_dk2*dk2_dsigxy
609 ENDIF
610 dfdsig_a = normxx * axx(i)
611 . + normyy * ayy(i)
612 . + normxy * axy(i)
613 IF (dfdsig_a /= zero) THEN
614 deplzz(i) = - max(dplanl(i),zero)*(yield/dfdsig_a)*(normxx + normyy)
615 ELSE
616 deplzz(i) = zero
617 ENDIF
618 ENDIF
619
620 depszz(i) = deelzz(i) + deplzz(i)
621 thk(i) = thk(i) + depszz(i)*thkly(i)*off(i)
622 soundsp(i) = sqrt(a1(i)/rho(i))
623 sigy(i) = uvar(i,4)
624 dep(i) = dep_dp(i)
625 ENDDO
626c-----------
627 RETURN
628 END
max
#define max(a, b)
Definition macros.h:21
sigeps78c
subroutine sigeps78c(nel, nuparam, nuvar, nvartmp, time, nfunc, ifunc, npf, tf, uparam, thkly, thk, gs, etse, sigy, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, soundsp, uvar, siga, sigb, sigc, rho, off, pla, dep, vartmp, inloc, dplanl, etimp, seq, loff)
Definition sigeps78c.F:40
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72