OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
no_asso_plas76.F File Reference
#include "implicit_f.inc"
#include "tabsiz_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine no_asso_plas76 (nel, nuparam, nuvar, nfunc, ifunc, nvartmp, npf, tf, time, timestep, uparam, vartmp, rho0, pla, dpla, et, numtabl, table, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, uvar, off, epsd, yld)

Function/Subroutine Documentation

◆ no_asso_plas76()

subroutine no_asso_plas76 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, intent(in) nvartmp,
integer, dimension(snpc) npf,
tf,
intent(in) time,
intent(in) timestep,
dimension(nuparam), intent(in) uparam,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
intent(in) rho0,
intent(inout) pla,
intent(out) dpla,
intent(out) et,
integer numtabl,
type(table_4d_), dimension(numtabl), intent(in) table,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) soundsp,
intent(inout) uvar,
intent(inout) off,
intent(inout) epsd,
intent(inout) yld )

Definition at line 34 of file no_asso_plas76.F.

42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE table4d_mod
46 USE table_mat_vinterp_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C G l o b a l P a r a m e t e r s
53C-----------------------------------------------
54#include "tabsiz_c.inc"
55C-----------------------------------------------
56C I N P U T A r g u m e n t s
57C-----------------------------------------------
58 INTEGER NEL,NUPARAM,NUVAR,NFUNC,NUMTABL
59 INTEGER ,INTENT(IN) :: NVARTMP
60 my_real,INTENT(IN) :: time,timestep
61 my_real,DIMENSION(NUPARAM),INTENT(IN) :: uparam(nuparam)
62 my_real,DIMENSION(NEL),INTENT(IN) ::
63 . rho0,depsxx,depsyy,depszz,
64 . depsxy,depsyz,depszx,sigoxx,sigoyy,sigozz,
65 . sigoxy,sigoyz,sigozx
66C-----------------------------------------------
67C O U T P U T A r g u m e n t s
68C-----------------------------------------------
69 my_real,DIMENSION(NEL),INTENT(OUT) ::
70 . signxx,signyy,signzz,signxy,signyz,signzx,
71 . soundsp,dpla,et
72C-----------------------------------------------
73C I N P U T O U T P U T A r g u m e n t s
74C-----------------------------------------------
75 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
76 my_real,DIMENSION(NEL),INTENT(INOUT) ::
77 . off,yld,pla,epsd
78 my_real,DIMENSION(NEL,NUVAR),INTENT(INOUT) ::
79 . uvar
80 TYPE(TABLE_4D_),DIMENSION(NUMTABL),INTENT(IN) :: TABLE
81C-----------------------------------------------
82C VARIABLES FOR FUNCTION INTERPOLATION
83C-----------------------------------------------
84 INTEGER :: NPF(SNPC), IFUNC(NFUNC)
85 my_real :: tf(stf)
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER I,J,II,ICONV,NINDX,ICAS,ITER,IQUAD
90 INTEGER,PARAMETER :: FUNC_TRAC = 1
91 INTEGER,PARAMETER :: FUNC_COMP = 2
92 INTEGER,PARAMETER :: FUNC_SHEAR = 3
93 INTEGER,DIMENSION(NEL) :: INDX,IAD,ILEN
94 my_real :: g,aa1,aa2,c1,nupc,xfac
95 my_real ,DIMENSION(NEL) :: plat,plac,plas,epspt,epspc,epsps,sigt,
96 . sigs,dsigs_dp,sigc,dsigc_dp,sdxx,sdyy,sdzz,sdxy,sdyz,sdzx,dsigt_dp,
97 . nup,p,svm,a0,a1,a2,phi,dplac,dplat,dplas,psi,alpha,dydx
98 my_real
99 . normxx,normyy,normzz,normxy,normyz,normzx,norm,aa,cb,
100 . normxx_n,normyy_n,normzz_n,normxy_n,normyz_n,normzx_n,
101 . dsigt_dlam,dsigc_dlam,dsigs_dlam,da1_dsigs,da1_dsigt,
102 . da1_dsigc,da2_dsigs,da2_dsigt,da2_dsigc,da0_dlam,da1_dlam,
103 . da2_dlam,dfdsig2,dphi_dlam,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,
104 . dpdt_c,dpdt_s,dpdt_t,dpdt,dlam,cc,epdt_min,da0_dsigs,
105 . epdt_max,epdc_max,epdc_min,epds_min,epds_max,asrate
106 my_real, DIMENSION(NEL,2) :: xvec
107 my_real, DIMENSION(NUMTABL) :: tfac
108 my_real, DIMENSION(NFUNC) :: yfac
109 LOGICAL :: CONV(NEL)
110 my_real, PARAMETER :: sfac = 1.05d0 ! Security factor of ICONV
111C
112 !=======================================================================
113 ! - INITIALISATION OF COMPUTATION ON TIME STEP
114 !=======================================================================
115 ! Recovering model parameters
116 ! -> Elastic parameters
117 g = uparam(4) ! Shear modulus
118 aa1 = uparam(6) ! First component of elastic matrix
119 aa2 = uparam(7) ! Second component of elastic matrix
120 c1 = uparam(8) ! Bulk modulus
121 ! -> Plastic parameters
122 nupc = uparam(9) ! Plastic Poisson ratio
123 ! -> Flags
124 iquad = nint(uparam(14)) ! Flag for quadratic or non-quatratic yield surface
125 asrate = min(one,uparam(16)*timestep)
126 iconv = nint(uparam(15)) ! Flag to ensure convexity
127 icas = nint(uparam(17)) ! Flag for tabulated case
128 !icas ifunt | ifunc | ifuncs
129 ! -1 1 | 1 | 1
130 ! 0 1 | 0 | 0
131 ! 1 1 | 1 | 0
132 ! 2 1 | 0 | 1
133 xfac = uparam(18) ! Strain-rate scale factor
134 epdt_min = uparam(19)
135 epdt_max = uparam(20)
136 epdc_min = uparam(21)
137 epdc_max = uparam(22)
138 epds_min = uparam(23)
139 epds_max = uparam(24)
140 tfac(1) = uparam(25)
141 tfac(2) = uparam(26)
142 tfac(3) = uparam(27)
143 yfac(1) = uparam(28)
144 yfac(2) = uparam(29)
145c
146 ! Initialize plastic Poisson ratio
147 IF (time == zero) THEN
148 nup(1:nel) = nupc
149 uvar(1:nel,7) = nupc
150 ELSE
151 nup(1:nel) = uvar(1:nel,7)
152 END IF
153c
154 ! Recovering internal variables
155 DO i=1,nel
156 plat(i) = uvar(i,1) ! Uniaxial tension plastic strain
157 plac(i) = uvar(i,2) ! Uniaxial compression plastic strain
158 plas(i) = uvar(i,3) ! Shear plastic strain
159 epspt(i) = uvar(i,4) ! Uniaxial tension plastic strain-rate
160 epspc(i) = uvar(i,5) ! Uniaxial compression plastic strain-rate
161 epsps(i) = uvar(i,6) ! Shear plastic strain-rate
162 epspt(i) = min(epdt_max, max(epspt(i),epdt_min))
163 epspc(i) = min(epdc_max, max(epspc(i),epdc_min))
164 epsps(i) = min(epds_max, max(epsps(i),epds_min))
165 alpha(i) = (nine/two)*((one-two*nup(i))/(one+nup(i)))
166 dpla(i) = zero
167 dplac(i) = zero
168 dplat(i) = zero
169 dplas(i) = zero
170 ENDDO
171C
172 ! Computation of yield stresses
173 xvec(1:nel,1) = plat(1:nel)
174 xvec(1:nel,2) = epspt(1:nel)*xfac
175 CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dsigt_dp)
176 sigt = sigt*tfac(1)
177 dsigt_dp = dsigt_dp*tfac(1)
178 IF (table(func_comp)%NOTABLE > 0) THEN
179 xvec(1:nel,1) = plac(1:nel)
180 xvec(1:nel,2) = epspc(1:nel)*xfac
181 CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dsigc_dp)
182 sigc = sigc*tfac(2)
183 dsigc_dp = dsigc_dp*tfac(2)
184 ENDIF
185 IF (table(func_shear)%NOTABLE > 0) THEN
186 xvec(1:nel,1) = plas(1:nel)
187 xvec(1:nel,2) = epsps(1:nel)*xfac
188 CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dsigs_dp)
189 sigs = sigs*tfac(3)
190 dsigs_dp = dsigs_dp*tfac(3)
191 ENDIF
192 ! Select case for tabulated yield stresses
193 IF (icas == 0) THEN
194 sigc(1:nel) = sigt(1:nel)
195 sigs(1:nel) = sigt(1:nel)/sqr3
196 ELSEIF (icas == 1) THEN
197 IF (iquad == 1) THEN
198 DO i = 1,nel
199 sigs(i) = sqrt(sigc(i)*sigt(i)/three)
200 ENDDO
201 ELSEIF (iquad == 0) THEN
202 DO i = 1,nel
203 sigs(i) = one /(sigt(i) + sigc(i))/sqr3
204 sigs(i) = two*sigt(i)*sigc(i)*sigs(i)
205 ENDDO
206 ENDIF
207 ENDIF
208C
209 !========================================================================
210 ! - COMPUTATION OF TRIAL VALUES
211 !========================================================================
212 DO i=1,nel
213C
214 ! Computation of the trial stress tensor
215 signxx(i) = sigoxx(i) + (aa1*depsxx(i)
216 . + aa2*(depsyy(i) + depszz(i)))
217 signyy(i) = sigoyy(i) + (aa1*depsyy(i)
218 . + aa2*(depsxx(i) + depszz(i)))
219 signzz(i) = sigozz(i) + (aa1*depszz(i)
220 . + aa2*(depsxx(i) + depsyy(i)))
221 signxy(i) = sigoxy(i) + g*depsxy(i)
222 signyz(i) = sigoyz(i) + g*depsyz(i)
223 signzx(i) = sigozx(i) + g*depszx(i)
224C
225 ! Computation of the pressure of the trial stress tensor
226 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
227 ! Computation of the Von Mises equivalent stress of the trial stress tensor
228 sdxx(i) = signxx(i) + p(i)
229 sdyy(i) = signyy(i) + p(i)
230 sdzz(i) = signzz(i) + p(i)
231 sdxy(i) = signxy(i)
232 sdyz(i) = signyz(i)
233 sdzx(i) = signzx(i)
234 svm(i) = half*(sdxx(i)**2 + sdyy(i)**2 + sdzz(i)**2)
235 . + (sdxy(i)**2 + sdzx(i)**2 + sdyz(i)**2)
236 svm(i) = sqrt(three*svm(i))
237 ENDDO
238c
239 !========================================================================
240 ! - COMPUTATION OF YIELD FONCTION
241 !========================================================================
242 ! Compute yield criterion parameters
243 IF (iconv == 1) THEN
244 ! Ensured convexity
245 DO i = 1,nel
246 conv(i) = .false.
247 aa = one /(sigt(i) + sigc(i))/sqr3
248 IF ((iquad == 1) .AND. (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three))) THEN
249 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
250 conv(i) = .true.
251 ELSEIF ((iquad == 0) .AND. (sigs(i) < sfac*two*sigt(i)*sigc(i)*aa)) THEN
252 sigs(i) = sfac*two*sigt(i)*sigc(i)*aa
253 conv(i) = .true.
254 ENDIF
255 ENDDO
256 ENDIF
257 ! Compute yield criterion parameters A0,A1,A2
258 IF (iquad == 1) THEN
259 DO i = 1,nel
260 aa = one/sigc(i)/sigt(i)
261 a0(i) = three*(sigs(i)**2)
262 a1(i) = nine*(sigs(i)**2)*(sigc(i) - sigt(i))*aa
263 a2(i) = nine*(sigc(i)*sigt(i) - three*(sigs(i)**2))*aa
264 ENDDO
265 ELSE
266 DO i = 1,nel
267 a0(i) = sigs(i)*sqr3
268 a1(i) = three*(((sigt(i)-sigc(i))/(sigt(i)+sigc(i))) -
269 . a0(i)*((sigt(i)-sigc(i))/(sigt(i)*sigc(i))))
270 a2(i) = eighteen*((one/(sigt(i)+sigc(i)))-a0(i)/(two*sigt(i)*sigc(i)))
271 ENDDO
272 ENDIF
273c
274 ! -> Checking plastic behavior for all elements
275 nindx = 0
276 DO i=1,nel
277 IF (iquad == 1) THEN
278 phi(i) = (svm(i)**2) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
279 ELSE
280 phi(i) = svm(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
281 ENDIF
282 IF (phi(i) >= zero .AND. off(i) == one) THEN
283 nindx = nindx + 1
284 indx(nindx) = i
285 ENDIF
286 ENDDO
287c
288 !====================================================================
289 ! - PLASTIC RETURN MAPPING WITH CUTTING PLANE METHOD
290 !====================================================================
291 IF (nindx > 0) THEN
292c
293 ! Loop over the iterations
294 DO iter = 1,3
295c
296 ! Loop over yielding elements
297 DO ii = 1,nindx
298 i = indx(ii)
299c
300 ! Note : in this part, the purpose is to compute for each iteration
301 ! a plastic multiplier allowing to update internal variables to satisfy
302 ! the consistency condition using the cutting plane algorithm
303 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
304 ! -> PHI : current value of yield function (known)
305 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
306 ! into account of internal variables kinetic :
307 ! plasticity ...
308c
309 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
310 !-------------------------------------------------------------
311C
312 ! Non associated flow surface
313 psi(i) = sqrt((svm(i)**2)+alpha(i)*p(i)**2)
314 psi(i) = max(psi(i),em20)
315C
316 ! Computation of normal to non-associated yield surface
317 normxx_n = ((three*sdxx(i)/(two*psi(i))) - third*(alpha(i)*p(i)/psi(i)))
318 normyy_n = ((three*sdyy(i)/(two*psi(i))) - third*(alpha(i)*p(i)/psi(i)))
319 normzz_n = ((three*sdzz(i)/(two*psi(i))) - third*(alpha(i)*p(i)/psi(i)))
320 normxy_n = (three*sdxy(i)/(two*psi(i)))
321 normyz_n = (three*sdyz(i)/(two*psi(i)))
322 normzx_n = (three*sdzx(i)/(two*psi(i)))
323C
324 ! Computation of the normal to the yield surface
325 IF (iquad == 1) THEN
326 normxx = three*sdxx(i) + third*(a1(i) + two*a2(i)*p(i))
327 normyy = three*sdyy(i) + third*(a1(i) + two*a2(i)*p(i))
328 normzz = three*sdzz(i) + third*(a1(i) + two*a2(i)*p(i))
329 normxy = three*sdxy(i)
330 normyz = three*sdyz(i)
331 normzx = three*sdzx(i)
332 ELSE
333 normxx = three_half*(sdxx(i)/svm(i)) + third*(a1(i) + two*a2(i)*p(i))
334 normyy = three_half*(sdyy(i)/svm(i)) + third*(a1(i) + two*a2(i)*p(i))
335 normzz = three_half*(sdzz(i)/svm(i)) + third*(a1(i) + two*a2(i)*p(i))
336 normxy = three_half*(sdxy(i)/svm(i))
337 normyz = three_half*(sdyz(i)/svm(i))
338 normzx = three_half*(sdzx(i)/svm(i))
339 ENDIF
340C
341 ! 2 - Computation of DPHI_DLAMBDA
342 !---------------------------------------------------------
343c
344 ! a) Derivative with respect stress increments tensor DSIG
345 ! --------------------------------------------------------
346 dfdsig2 = normxx*(aa1*normxx_n + aa2*(normyy_n + normzz_n)) +
347 . normyy*(aa1*normyy_n + aa2*(normxx_n + normzz_n)) +
348 . normzz*(aa1*normzz_n + aa2*(normxx_n + normyy_n)) +
349 . two*normxy*two*g*normxy_n +
350 . two*normyz*two*g*normyz_n +
351 . two*normzx*two*g*normzx_n
352c
353 ! b) Derivative of yield criterion parameters
354 ! --------------------------------------------
355C
356 ! Derivative of yield surfaces with respect to yield criterion parameter A0,A1,A2
357 dsigt_dlam = dsigt_dp(i)*((svm(i)/psi(i))*three_half/(one + nup(i)))
358 IF (table(func_comp)%NOTABLE > 0) dsigc_dlam = dsigc_dp(i)*((svm(i)/psi(i))*three_half/(one + nup(i)))
359 IF (table(func_shear)%NOTABLE > 0) dsigs_dlam = dsigs_dp(i)*(sqr3/two)*(svm(i)/psi(i))
360 IF (icas == 0) THEN
361 dsigc_dlam = dsigt_dlam
362 dsigs_dlam = (one/sqr3)*dsigt_dlam
363 ELSEIF (icas == 1) THEN
364 IF (iquad == 1) THEN
365 dsigs_dlam = (one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
366 . (dsigc_dlam*sigt(i) + sigc(i)*dsigt_dlam)
367 ELSEIF (iquad == 0) THEN
368 aa = one /(sigt(i) + sigc(i))/sqr3
369 dsigs_dlam = two*(dsigt_dlam*sigc(i) + dsigc_dlam*sigt(i))*aa
370 . - two*sqr3*sigc(i)*sigt(i)*(dsigt_dlam + dsigc_dlam)*aa*aa
371 ENDIF
372 ENDIF
373 IF (iconv == 1) THEN
374 IF (conv(i)) THEN
375 IF (iquad == 1) THEN
376 dsigs_dlam = sfac*(one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
377 . (dsigc_dlam*sigt(i) + sigc(i)*dsigt_dlam)
378 ELSEIF (iquad == 0) THEN
379 aa = one /(sigt(i) + sigc(i))/sqr3
380 dsigs_dlam = sfac*(two*(dsigt_dlam*sigc(i) + dsigc_dlam*sigt(i))*aa
381 . - two*sqr3*sigc(i)*sigt(i)*(dsigt_dlam + dsigc_dlam)*aa*aa)
382 ENDIF
383 ENDIF
384 ENDIF
385C
386 IF (iquad == 1) THEN
387 ! -> A0 derivatives
388 da0_dsigs = six*sigs(i)
389 ! -> A1 derivatives
390 cc = sigs(i)/sigc(i)/sigt(i)
391 ! -> With respect to SIGS
392 da1_dsigs = eighteen*(sigc(i) - sigt(i))*cc
393 ! -> With respect to SIGC
394 da1_dsigc = nine*(sigs(i)/sigc(i))**2
395 ! -> With respect to SIGT
396 da1_dsigt = -nine*(sigs(i)/sigt(i))**2
397 ! -> A2 derivatives
398 ! -> With respect to SIGS
399 da2_dsigs = -cinquante4*cc
400 ! -> With respect to SIGC
401 da2_dsigc = twenty7*cc*sigs(i)/sigc(i)
402 ! -> With respect to SIGT
403 da2_dsigt = twenty7*cc*sigs(i)/sigt(i)
404 ELSE
405 ! -> A0 derivatives
406 da0_dsigs = sqr3
407 ! -> A1 derivatives
408 ! -> With respect to SIGS
409 da1_dsigs = -three*sqr3*(sigt(i)-sigc(i))/(sigt(i)*sigc(i))
410 ! -> With respect to SIGC
411 da1_dsigc = three*((sigs(i)*sqr3/(sigc(i)**2))-
412 . two*sigt(i)/((sigt(i) + sigc(i))**2))
413 ! -> With respect to SIGT
414 da1_dsigt = three*(two*sigc(i)/((sigt(i) + sigc(i))**2) -
415 . (sigs(i)*sqr3/(sigt(i)**2)))
416 ! -> a2 derivatives
417 ! -> With respect to SIGS
418 da2_dsigs = -nine*sqr3/(sigt(i)*sigc(i))
419 ! -> With respect to SIGC
420 da2_dsigc = eighteen*((sigs(i)*sqr3/(two*sigt(i)*(sigc(i)**2))) -
421 . one/((sigt(i)+sigc(i))**2))
422 ! -> With respect to SIGT
423 da2_dsigt = eighteen*((sigs(i)*sqr3/(two*sigc(i)*(sigt(i)**2))) -
424 . one/((sigt(i)+sigc(i))**2))
425 ENDIF
426c
427 ! -> A parameters derivatives with respect to lambda
428 ! -> A0 with respect to lambda
429 da0_dlam = da0_dsigs*dsigs_dlam
430 ! -> A1 with respect to lambda
431 da1_dlam = da1_dsigs*dsigs_dlam + da1_dsigt*dsigt_dlam + da1_dsigc*dsigc_dlam
432 ! -> A2 with respect to lambda
433 da2_dlam = da2_dsigs*dsigs_dlam + da2_dsigt*dsigt_dlam + da2_dsigc*dsigc_dlam
434c
435c
436 ! 3 - Computation of plastic multiplier and variables update
437 !----------------------------------------------------------
438c
439 ! Derivative of PHI with respect to DLAM
440 dphi_dlam = - dfdsig2 - da0_dlam - p(i)*da1_dlam - (p(i)**2)*da2_dlam
441 dphi_dlam = sign(max(abs(dphi_dlam),em20),dphi_dlam)
442 dlam = -phi(i)/dphi_dlam
443c
444 ! Plastic strains tensor update
445 dpxx = dlam*normxx_n
446 dpyy = dlam*normyy_n
447 dpzz = dlam*normzz_n
448 dpxy = dlam*normxy_n
449 dpyz = dlam*normyz_n
450 dpzx = dlam*normzx_n
451c
452 ! Elasto-plastic stresses update
453 signxx(i) = signxx(i) - (aa1*dpxx + aa2*(dpyy + dpzz))
454 signyy(i) = signyy(i) - (aa1*dpyy + aa2*(dpxx + dpzz))
455 signzz(i) = signzz(i) - (aa1*dpzz + aa2*(dpxx + dpyy))
456 signxy(i) = signxy(i) - two*g*dpxy
457 signyz(i) = signyz(i) - two*g*dpyz
458 signzx(i) = signzx(i) - two*g*dpzx
459c
460 ! Cumulated plastic strains update
461 plat(i) = plat(i) + dlam*(svm(i)/psi(i))*three_half/(one + nup(i))
462 plac(i) = plat(i)
463 plas(i) = plas(i) + (sqr3/two)*(svm(i)/psi(i))*dlam
464 pla(i) = pla(i) + (svm(i)/psi(i))*dlam
465 ! Plastic strain increments update
466 dplat(i) = dplat(i) + dlam*(svm(i)/psi(i))*three_half/(one + nup(i))
467 dplac(i) = dplat(i)
468 dplas(i) = dplas(i) + (sqr3/two)*(svm(i)/psi(i))*dlam
469 dpla(i) = dpla(i) + (svm(i)/psi(i))*dlam
470C
471 ! Update pressure
472 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
473 ! Update Von Mises stress
474 sdxx(i) = signxx(i) + p(i)
475 sdyy(i) = signyy(i) + p(i)
476 sdzz(i) = signzz(i) + p(i)
477 sdxy(i) = signxy(i)
478 sdyz(i) = signyz(i)
479 sdzx(i) = signzx(i)
480 svm(i) = half*(sdxx(i)**2 + sdyy(i)**2 + sdzz(i)**2)
481 . + (sdxy(i)**2 + sdzx(i)**2 + sdyz(i)**2)
482 svm(i) = sqrt(three*svm(i))
483 ENDDO
484C
485 ! Update yield stresses
486 xvec(1:nel,1) = plat(1:nel)
487 xvec(1:nel,2) = epspt(1:nel)*xfac
488 CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dsigt_dp)
489 sigt = sigt*tfac(1)
490 dsigt_dp = dsigt_dp*tfac(1)
491 IF (table(func_comp)%NOTABLE > 0) THEN
492 xvec(1:nel,1) = plac(1:nel)
493 xvec(1:nel,2) = epspc(1:nel)*xfac
494 CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dsigc_dp)
495 sigc(1:nel) = sigc*tfac(2)
496 dsigc_dp(1:nel) = dsigc_dp*tfac(2)
497 ENDIF
498 IF (table(func_shear)%NOTABLE > 0) THEN
499 xvec(1:nel,1) = plas(1:nel)
500 xvec(1:nel,2) = epsps(1:nel)*xfac
501 CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dsigs_dp)
502 sigs(1:nel) = sigs*tfac(3)
503 dsigs_dp(1:nel) = dsigs_dp*tfac(3)
504 ENDIF
505 ! Select case for tabulated yield stresses
506 IF (icas == 0) THEN
507 DO ii = 1,nindx
508 i = indx(ii)
509 sigc(i) = sigt(i)
510 sigs(i) = sigt(i)/sqr3
511 ENDDO
512 ELSEIF (icas == 1) THEN
513 IF (iquad == 1) THEN
514 DO ii = 1,nindx
515 i = indx(ii)
516 sigs(i) = sqrt(sigc(i)*sigt(i)/three)
517 ENDDO
518 ELSEIF (iquad == 0) THEN
519 DO ii = 1,nindx
520 i = indx(ii)
521 sigs(i) = one /(sigt(i) + sigc(i))/sqr3
522 sigs(i) = two*sigt(i)*sigc(i)*sigs(i)
523 ENDDO
524 ENDIF
525 ENDIF
526c
527 ! Update yield function value
528 IF (iconv == 1) THEN
529 DO ii = 1,nindx
530 i = indx(ii)
531 aa = one /(sigt(i) + sigc(i))/sqr3
532 conv(i) = .false.
533 IF ((iquad == 1) .AND. (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three))) THEN
534 sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
535 conv(i) = .true.
536 ELSEIF ((iquad == 0) .AND. (sigs(i) < sfac*two*sigt(i)*sigc(i)*aa)) THEN
537 sigs(i) = sfac*two*sigt(i)*sigc(i)*aa
538 conv(i) = .true.
539 ENDIF
540 ENDDO
541 ENDIF
542 DO ii = 1,nindx
543 i = indx(ii)
544 IF (iquad == 1) THEN
545 aa = one/sigc(i)/sigt(i)
546 a0(i) = three*(sigs(i)**2)
547 a1(i) = nine*(sigs(i)**2)*(sigc(i) - sigt(i))*aa
548 a2(i) = nine*(sigc(i)*sigt(i) - three*(sigs(i)**2))*aa
549 ELSE
550 a0(i) = sigs(i)*sqr3
551 a1(i) = three*(((sigt(i)-sigc(i))/(sigt(i)+sigc(i))) -
552 . a0(i)*((sigt(i)-sigc(i))/(sigt(i)*sigc(i))))
553 a2(i) = eighteen*((one/(sigt(i)+sigc(i)))-a0(i)/(two*sigt(i)*sigc(i)))
554 ENDIF
555 IF (iquad == 1) THEN
556 phi(i) = (svm(i)**2) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
557 ELSE
558 phi(i) = svm(i) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
559 ENDIF
560 ENDDO
561 ENDDO
562 ENDIF
563c
564 !====================================================================
565 ! - UPDATE PLASTIC POISSON RATIO
566 !====================================================================
567 IF (ifunc(1) > 0) THEN
568 iad(1:nel) = npf(ifunc(1)) / 2 + 1
569 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - vartmp(1:nel,8)
570!
571 CALL vinter(tf,iad,vartmp(1:nel,8),ilen,nel,pla,dydx,nup)
572!
573 uvar(1:nel,7) = yfac(1) * nup(1:nel)
574 uvar(1:nel,7) = max(zero, min(nup(1:nel), half))
575 END IF
576c
577 !====================================================================
578 ! - STORING NEW VALUES
579 !====================================================================
580 DO i=1,nel
581 ! Update user variable
582 uvar(i,1) = plat(i)
583 uvar(i,2) = plac(i)
584 uvar(i,3) = plas(i)
585 dpdt_t = dplat(i)/max(timestep,em20)
586 uvar(i,4) = asrate*dpdt_t + (one-asrate)*epspt(i)
587 dpdt_c = dplac(i)/max(timestep,em20)
588 uvar(i,5) = asrate*dpdt_c + (one-asrate)*epspc(i)
589 dpdt_s = dplas(i)/max(timestep,em20)
590 uvar(i,6) = asrate*dpdt_s + (one-asrate)*epsps(i)
591 dpdt = dpla(i)/max(timestep,em20)
592 epsd(i) = asrate*dpdt + (one-asrate)*epsd(i)
593 ! Computation of soundspeed
594 soundsp(i) = sqrt((c1+four*g/three)/rho0(i))
595 ! Computation of the yield stress
596 yld(i) = a0(i) + a1(i)*p(i) + a2(i)*p(i)*p(i)
597 ! Computation of the hourglass coefficient
598 et(i) = half
599 ENDDO
600c
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
func_comp
subroutine func_comp(table, ntable, nptmax, npt_trac, npt_shear, npt_comp, x_comp, y_comp, nup)
Definition law76_upd.F:344
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
table4d_mod
Definition table4d_mod.F:52
table_mat_vinterp_mod
Definition table_mat_vinterp.F:32
table_mat_vinterp_mod::table_mat_vinterp
subroutine table_mat_vinterp(table, dimx, nel, ipos, xx, yy, dydx)
Definition table_mat_vinterp.F:45
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72