OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
no_asso_plas76.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!|| no_asso_plas76 ../engine/source/materials/mat/mat076/no_asso_plas76.F
25!||--- called by ------------------------------------------------------
26!|| sigeps76 ../engine/source/materials/mat/mat076/sigeps76.F
27!||--- calls -----------------------------------------------------
28!|| table_mat_vinterp ../engine/source/materials/tools/table_mat_vinterp.F
29!|| vinter ../engine/source/tools/curve/vinter.F
30!||--- uses -----------------------------------------------------
31!|| table4d_mod ../common_source/modules/table4d_mod.F
32!|| table_mat_vinterp_mod ../engine/source/materials/tools/table_mat_vinterp.F
33!||====================================================================
34 SUBROUTINE no_asso_plas76(
35 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,NVARTMP ,
36 2 NPF ,TF ,TIME ,TIMESTEP,UPARAM ,VARTMP ,
37 3 RHO0 ,PLA ,DPLA ,ET ,NUMTABL ,TABLE ,
38 3 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
39 4 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 6 SOUNDSP ,UVAR ,OFF ,EPSD ,YLD )
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
601 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
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
no_asso_plas76
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)
Definition no_asso_plas76.F:42
table4d_mod::table_4d_
Definition table4d_mod.F:60
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72