OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2cplr.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!|| m2cplr ../engine/source/materials/mat/mat002/m2cplr.F
25!||--- called by ------------------------------------------------------
26!|| sigeps02c ../engine/source/materials/mat/mat002/sigeps02c.F
27!||====================================================================
28 SUBROUTINE m2cplr(NEL ,EZZ ,OFF ,PLA ,
29 2 IPLA ,TEMP ,C3 ,C4 ,M_EXP ,
30 3 IFORM ,ETSE ,GS ,EPSP ,VP ,
31 4 ISRATE ,YLD ,G ,A1 ,A2 ,
32 5 NU ,CA0 ,CB0 ,CN ,YMAX0 ,
33 6 EPCHK ,YOUNG ,CC ,EPDR ,ICC ,
34 7 DPLA ,TSTAR ,FISOKIN ,GAMA_IMP ,SIGNOR ,
35 8 HARDM ,DEPSXX ,DEPSYY ,DEPSXY ,
36 9 DEPSYZ ,DEPSZX ,SIGNXX ,SIGNYY ,SIGNXY ,
37 A SIGNYZ ,SIGNZX ,SIGBAKXX ,SIGBAKYY ,SIGBAKXY,
38 B SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX )
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C G l o b a l P a r a m e t e r s
45C-----------------------------------------------
46#include "mvsiz_p.inc"
47C-----------------------------------------------
48C C o m m o n B l o c k s
49C-----------------------------------------------
50#include "impl1_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER ,INTENT(IN) :: NEL
55 INTEGER ,INTENT(IN) :: IPLA
56 INTEGER ,INTENT(IN) :: IFORM
57 INTEGER ,INTENT(IN) :: ISRATE
58 INTEGER ,INTENT(IN) :: ICC
59 INTEGER ,INTENT(IN) :: VP
60 my_real ,INTENT(IN) :: C3,C4
61 my_real ,INTENT(IN) :: M_EXP
62 my_real ,INTENT(IN) :: G
63 my_real ,INTENT(IN) :: NU
64 my_real ,INTENT(IN) :: a1
65 my_real ,INTENT(IN) :: a2
66 my_real ,INTENT(IN) :: ca0
67 my_real ,INTENT(IN) :: cb0
68 my_real ,INTENT(IN) :: cn
69 my_real ,INTENT(IN) :: ymax0
70 my_real ,INTENT(IN) :: young
71 my_real ,INTENT(IN) :: cc
72 my_real ,INTENT(IN) :: epdr
73 my_real ,INTENT(IN) :: fisokin
74 my_real, DIMENSION(NEL) ,INTENT(IN) :: depsxx,depsyy,depsxy,depsyz,depszx
75 my_real, DIMENSION(NEL) ,INTENT(IN) :: sigoxx,sigoyy,sigoxy,sigoyz,sigozx
76 my_real, DIMENSION(NEL) ,INTENT(IN) :: gs
77 my_real, DIMENSION(NEL) ,INTENT(IN) :: tstar
78 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: signxx,signyy,signxy,signyz,signzx
79 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: sigbakxx,sigbakyy,sigbakxy,hardm
80 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: ezz,off,pla,etse,epchk,yld,epsp
81 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: temp,dpla,gama_imp
82 my_real, DIMENSION(NEL,5),INTENT(INOUT) :: signor
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER :: I,J,N,NINDX,NMAX
87 INTEGER :: INDEX(MVSIZ)
88 my_real :: f,df,pla_i,p2,q2,r,s1,s2,s3,yld_i,small,beta,aa,aaa,hkin,alpha
89 my_real :: a(mvsiz), b(mvsiz) ,dpla_i(mvsiz), dpla_j(mvsiz), dr(mvsiz),
90 . h(mvsiz), nu1(mvsiz), nu2(mvsiz) ,p(mvsiz) ,q(mvsiz),
91 . svm(mvsiz),ca(mvsiz), cb(mvsiz) ,ymax(mvsiz) ,
92 . svm2(mvsiz),yld2(mvsiz),hi(mvsiz),hk(mvsiz),logep(mvsiz),
93 . nu11(mvsiz),nu21(mvsiz),anu1(mvsiz),bnu2(mvsiz),h2(mvsiz)
94 DATA nmax/3/
95!===============================================================================
96 nindx = 0
97 small = em7
98 DO i = 1,nel
99 etse(i) = one
100 ca(i) = ca0
101 cb(i) = cb0
102 ymax(i) = ymax0
103 h(i) = zero
104!
105 signxx(i)=sigoxx(i)
106 signyy(i)=sigoyy(i)
107 signxy(i)=sigoxy(i)
108 signyz(i)=sigoyz(i)
109 signzx(i)=sigozx(i)
110 ENDDO
111C------------------------------------------
112C KINEMATIC HARDENING
113C------------------------------------------
114 IF (fisokin > zero) THEN
115 DO i = 1,nel
116 signxx(i)=signxx(i)-sigbakxx(i)
117 signyy(i)=signyy(i)-sigbakyy(i)
118 signxy(i)=signxy(i)-sigbakxy(i)
119 ENDDO
120 END IF ! FISOKIN > 0
121C---------------------------
122C ELASTIC STRESSES
123C---------------------------
124 DO i = 1,nel
125 signxx(i)=signxx(i)+a1*depsxx(i)+a2*depsyy(i)
126 signyy(i)=signyy(i)+a2*depsxx(i)+a1*depsyy(i)
127 signxy(i)=signxy(i)+g*depsxy(i)
128 signyz(i)=signyz(i)+gs(i)*depsyz(i)
129 signzx(i)=signzx(i)+gs(i)*depszx(i)
130 ENDDO
131 DO i = 1,nel
132 logep(i) = zero
133 ENDDO
134C-------------
135C STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
136C-------------
137 IF (cc /= zero) THEN
138 IF (iform == 0) THEN
139 DO i = 1,nel
140 IF (israte == 0.AND.vp == 2) epsp(i) =
141 . max( abs(depsxx(i)), abs(depsyy(i)), half*abs(depsxy(i)))
142 epsp(i) = max(epsp(i),epdr)
143 logep(i) = log(epsp(i)/epdr)
144 ENDDO
145#include "vectorize.inc"
146 DO i = 1,nel
147 IF (tstar(i) == zero) THEN
148 q(i) = (one + cc * logep(i))
149 ELSE
150 q(i) = (one + cc * logep(i))*(one-exp(m_exp*log(tstar(i))))
151 ENDIF
152 q(i) = max(q(i),em20)
153 ca(i) = ca(i) * q(i)
154 cb(i) = cb(i) * q(i)
155 IF (icc == 1) ymax(i) = ymax(i) * q(i)
156 ENDDO
157 ELSEIF (iform == 1) THEN ! Zerilli
158 DO i = 1,nel
159 IF (israte == 0 .AND. vp == 2) epsp(i) = max( abs(depsxx(i)),
160 . abs(depsyy(i)), half*abs(depsxy(i)))
161 epsp(i) = max(epsp(i),em20)
162 logep(i) = log(epsp(i)/epdr)
163 ENDDO
164 DO i = 1,nel
165 q(i) = logep(i)
166 q(i) = cc*exp((-c3 + c4*q(i))*temp(i))
167 IF (icc == 1) ymax(i) = ymax(i) + q(i)
168 ca(i) = ca(i) + q(i)
169 ENDDO
170 ENDIF ! IFORM
171 ELSE ! CC = ZERO , no strain rate effect
172 IF (iform == 0) THEN ! J-C
173 DO i = 1,nel
174 IF (tstar(i) /= zero) THEN
175 q(i) = one - exp(m_exp*log(tstar(i)))
176 q(i) = max(q(i),em20)
177 ca(i) = ca(i) * q(i)
178 cb(i) = cb(i) * q(i)
179 END IF
180 ENDDO
181 ENDIF ! IFORM
182 ENDIF ! IF (CC /= ZERO)
183C-----------------------------------
184C YIELD
185C-----------------------------------
186 DO i = 1,nel
187 dpla(i) = zero
188 IF (pla(i) == zero) THEN
189 yld(i)= ca(i)
190 ELSE
191 beta = cb(i)*(one-fisokin)
192 yld(i)= ca(i)+beta*exp(cn*log(pla(i)))
193 ENDIF
194 yld(i)= min(yld(i),ymax(i))
195 ENDDO
196C------------------------------------------
197C PLASTICALLY ADMISSIBLE STRESSES
198C------------------------------------------
199 IF (ipla == 0) THEN
200C-------------------
201C VM STRESS
202C-------------------
203 DO i = 1,nel
204 svm(i)=sqrt(signxx(i)*signxx(i)
205 . +signyy(i)*signyy(i)
206 . -signxx(i)*signyy(i)
207 . +three*signxy(i)*signxy(i))
208 ENDDO
209C radial projection
210#include "vectorize.inc"
211 DO i = 1,nel
212 r = min(one,yld(i)/(svm(i)+em15))
213 IF (r < one) THEN
214 signxx(i)=signxx(i)*r
215 signyy(i)=signyy(i)*r
216 signxy(i)=signxy(i)*r
217 dpla(i) = off(i) * max(zero,(svm(i)-yld(i))/young)
218 s1=half*(signxx(i)+signyy(i))
219 ezz(i) = dpla(i) * s1 /yld(i)
220 pla(i) = pla(i) + dpla(i)
221 epchk(i) = max(pla(i),epchk(i))
222 IF (yld(i) >=ymax(i)) THEN
223 h(i)=zero
224 ELSE
225 h(i)=cn*cb(i)*exp((cn-one)*log(pla(i)+small))
226 ENDIF
227 etse(i)= h(i)/(h(i)+young)
228 ENDIF
229 ENDDO
230C------------------------------------------------------------------------
231 ELSEIF (ipla == 1) THEN
232C-------------------------
233C von mises criterion
234C-------------------------
235 DO i = 1,nel
236 s1=signxx(i)+signyy(i)
237 s2=signxx(i)-signyy(i)
238 s3=signxy(i)
239 a(i)=fourth*s1*s1
240 b(i)=three_over_4*s2*s2+three*s3*s3
241 svm(i)=sqrt(a(i)+b(i))
242 ENDDO
243C-------------------------
244C GATHER PLASTIC FLOW
245C-------------------------
246 nindx=0
247 DO i = 1,nel
248 IF (svm(i) > yld(i) .AND. off(i) == one) THEN
249 nindx=nindx+1
250 index(nindx)=i
251 ENDIF
252 ENDDO
253 IF (nindx == 0) THEN
254 IF (fisokin > zero) THEN
255 DO i = 1,nel
256 signxx(i)=signxx(i) + sigbakxx(i)
257 signyy(i)=signyy(i) + sigbakyy(i)
258 signxy(i)=signxy(i) + sigbakxy(i)
259 ENDDO
260 END IF !(FISOKIN > ZERO) THEN
261 IF (impl_s > 0.AND.ikt > 0) THEN
262 DO i = 1,nel
263 gama_imp(i) = zero
264 END DO
265 END IF !(IMPL_S > 0.AND.IKT > 0) THEN
266 RETURN
267 END IF !(NINDX==0) THEN
268C---------------------------
269C DEP IN PLANE STRESS
270C---------------------------
271#include "vectorize.inc"
272 DO j=1,nindx
273 i=index(j)
274 nu1(i)=one/(one-nu)
275 nu2(i)=one/(one+nu)
276 IF (yld(i) >= ymax(i)) THEN
277 h(i)=zero
278 ELSE
279 h(i)=cn*cb(i)*exp((cn-one)*log(pla(i)+small))
280 ENDIF
281 dpla_j(i)=(svm(i)-yld(i))/(three*g+h(i))
282 etse(i)= h(i)/(h(i)+young)
283 anu1(i) = a(i)*nu1(i)
284 bnu2(i) = three*b(i)*nu2(i)
285 h2(i)= two*h(i)
286 ENDDO
287C-----just keep at least old performance--to update H(I) inside iteration for implicit+ETSE---------
288 IF (fisokin == zero) THEN
289 DO n=1,nmax
290#include "vectorize.inc"
291 DO j=1,nindx
292 i=index(j)
293 dpla_i(i)=dpla_j(i)
294 pla_i =pla(i)+dpla_i(i)
295 dpla(i) = dpla_j(i)
296 IF (pla_i == zero) THEN
297 yld_i =min(ymax(i),ca(i))
298 ELSE
299 yld_i =min(ymax(i),ca(i)+cb(i)*exp(cn*log(pla_i)))
300 ENDIF
301 dr(i) =half*young*dpla_i(i)/yld_i
302 p(i) =one/(one+dr(i)*nu1(i))
303 q(i) =one/(one+three*dr(i)*nu2(i))
304 p2 =p(i)*p(i)
305 q2 =q(i)*q(i)
306 f =a(i)*p2+b(i)*q2-yld_i*yld_i
307 df =-(anu1(i)*p2*p(i)+bnu2(i)*q2*q(i))
308 . *(young-dr(i)*h2(i))/yld_i
309 . -h2(i)*yld_i
310 IF (dpla_i(i) > zero) THEN
311 dpla_j(i)=max(zero,dpla_i(i)-f/df)
312 ELSE
313 dpla_j(i)=zero
314 ENDIF
315 ENDDO ! DO J=1,NINDX
316 ENDDO ! DO N=1,NMAX
317!
318 ELSE ! FISOKIN > 0 -> kinematic/mix hardening
319#include "vectorize.inc"
320 DO j=1,nindx
321 i=index(j)
322 beta = h(i)*fisokin
323 hi(i) = h(i)-beta
324 hk(i) = two_third*beta
325 aaa = three*hk(i)/young
326 nu11(i) = nu1(i) + aaa
327 nu21(i) = three*nu2(i) + aaa
328 anu1(i) = a(i)*nu11(i)
329 bnu2(i) = b(i)*nu21(i)
330 h2(i)= two*hi(i)
331 ENDDO
332!
333 DO n=1,nmax
334#include "vectorize.inc"
335 DO j=1,nindx
336 i=index(j)
337 dpla_i(i)=dpla_j(i)
338 pla_i =pla(i)+dpla_i(i)
339 dpla(i) = dpla_j(i)
340 beta = one-fisokin
341 IF (pla_i == zero) THEN
342 yld_i =min(ymax(i),ca(i))
343 ELSE
344 yld_i =min(ymax(i),ca(i)+beta*cb(i)*exp(cn*log(pla_i)))
345 ENDIF
346 dr(i) =half*young*dpla_i(i)/yld_i
347 p(i) =one/(one+dr(i)*nu11(i))
348 q(i) =one/(one+dr(i)*nu21(i))
349 p2 =p(i)*p(i)
350 q2 =q(i)*q(i)
351 f =a(i)*p2+b(i)*q2-yld_i*yld_i
352 df =-(anu1(i)*p2*p(i)+bnu2(i)*q2*q(i))
353 . *(young-dr(i)*h2(i))/yld_i
354 . -h2(i)*yld_i
355 IF (dpla_i(i) > zero) THEN
356 dpla_j(i)=max(zero,dpla_i(i)-f/df)
357 ELSE
358 dpla_j(i)=zero
359 ENDIF
360 ENDDO ! DO J=1,NINDX
361 ENDDO ! DO N=1,NMAX
362 ENDIF ! FISOKIN > 0
363C------------------------------------------
364C PLASTICALLY ADMISSIBLE STRESSES
365C------------------------------------------
366#include "vectorize.inc"
367 DO j=1,nindx
368 i=index(j)
369 pla(i) = pla(i) + dpla_i(i)
370 epchk(i) = max(pla(i),epchk(i))
371 s1=(signxx(i)+signyy(i))*p(i)
372 s2=(signxx(i)-signyy(i))*q(i)
373 signxx(i)=half*(s1+s2)
374 signyy(i)=half*(s1-s2)
375 signxy(i)=signxy(i)*q(i)
376 ezz(i) = dr(i)*s1/young
377 ENDDO
378C-------------------------
379 ELSEIF (ipla == 2) THEN
380C-------------------------
381C von mises criterion
382C-------------------------
383 DO i = 1,nel
384 svm2(i)= signxx(i)*signxx(i)
385 . +signyy(i)*signyy(i)
386 . -signxx(i)*signyy(i)
387 . +three*signxy(i)*signxy(i)
388 svm(i)=sqrt(svm2(i))
389 END DO
390C-------------------------
391C GATHER PLASTIC FLOW
392C-------------------------
393 nindx=0
394 DO i = 1,nel
395 yld2(i)=yld(i)*yld(i)
396 IF (svm2(i) > yld2(i) .AND. off(i) == one) THEN
397 nindx=nindx+1
398 index(nindx)=i
399 END IF
400 END DO
401!
402 IF (nindx /= 0) THEN
403C-------------
404C NORMAL PROJECTION TO THE CRITERION WITH APPROXIMATE CALCULATION OF THE NORMAL + RADIAL RETURN
405C-------------
406#include "vectorize.inc"
407 DO j=1,nindx
408 i=index(j)
409 IF (yld(i) >= ymax(i)) THEN
410 h(i)=zero
411 ELSE
412 h(i)=cn*cb(i)*exp((cn-one)*log(pla(i)+small))
413 ENDIF
414 etse(i)= h(i)/(h(i)+young)
415C
416 aa=(svm2(i)-yld2(i))
417 . /(five*svm2(i)+three*(-signxx(i)*signyy(i)+signxy(i)*signxy(i)))
418 s1=(one-two*aa)*signxx(i)+ aa*signyy(i)
419 s2=aa*signxx(i)+(one-two*aa)*signyy(i)
420 s3=(one-three*aa)*signxy(i)
421 signxx(i)=s1
422 signyy(i)=s2
423 signxy(i)=s3
424 dpla(i) = off(i)*(svm(i)-yld(i))/(three*g+h(i))
425 pla(i) = pla(i) + dpla(i)
426C
427 yld(i) =yld(i)+h(i)*dpla(i)
428 END DO ! DO J=1,NINDX
429C
430#include "vectorize.inc"
431 DO j=1,nindx
432 i=index(j)
433 svm(i)=sqrt( signxx(i)*signxx(i)
434 . +signyy(i)*signyy(i)
435 . -signxx(i)*signyy(i)
436 . +three*signxy(i)*signxy(i))
437 r = min(one,yld(i)/max(em20,svm(i)))
438 signxx(i)=signxx(i)*r
439 signyy(i)=signyy(i)*r
440 signxy(i)=signxy(i)*r
441 ezz(i) = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
442 END DO ! DO J=1,NINDX
443 END IF ! IF (NINDX /= 0)
444 ENDIF ! IF (IPLA == 0)
445C------------------------------------------
446C for tangent matrix
447C------------------------------------------
448 IF (impl_s > 0) THEN
449 IF (ikt > 0) THEN
450 DO i = 1,nel
451 IF (dpla(i) > zero) THEN
452c ...... Parameter d(gama)
453 pla_i =pla(i)
454 beta = one-fisokin
455 yld(i) =min(ymax(i),ca(i)+beta*cb(i)*exp(cn*log(pla_i)))
456 gama_imp(i)= three_half*dpla(i)/yld(i)
457c ...... HK,HH---
458 signor(i,4)=fisokin*h(i)
459 signor(i,5)=(one-fisokin)*h(i)
460c ...... Deviatoric stresses shifted by modified back stress ->ksi
461 signor(i,1)=third*(two*signxx(i)-signyy(i))
462 signor(i,2)=third*(two*signyy(i)-signxx(i))
463 signor(i,3)=two*signxy(i)
464 ELSE
465 gama_imp(i) = zero
466 END IF ! IF (DPLA(I) > ZERO)
467 END DO
468 END IF ! (IKT > 0) THEN
469 END IF ! IF (IMPL_S > 0)
470C------------------------------------------
471C ECROUISSAGE CINE
472C------------------------------------------
473 IF (fisokin > zero) THEN
474 IF (ipla == 1 )THEN
475#include "vectorize.inc"
476 DO j=1,nindx
477 i=index(j)
478 pla_i =pla(i)
479 beta = one-fisokin
480C------YLD, H(I) should not be updated----
481 IF (pla_i == zero) THEN
482 yld(i) =ca(i)
483 ELSE
484 yld(i) =min(ymax(i),ca(i)+beta*cb(i)*exp(cn*log(pla_i)))
485 ENDIF
486 END DO ! DO J=1,NINDX
487 END IF ! (IPLA == 1) THEN
488 DO i = 1,nel
489 hkin = fisokin*h(i)
490 alpha = hkin*dpla(i)/yld(i)
491 sigbakxx(i) = sigbakxx(i) + alpha*signxx(i)
492 sigbakyy(i) = sigbakyy(i) + alpha*signyy(i)
493 sigbakxy(i) = sigbakxy(i) + alpha*signxy(i)
494C
495 signxx(i)=signxx(i) + sigbakxx(i)
496 signyy(i)=signyy(i) + sigbakyy(i)
497 signxy(i)=signxy(i) + sigbakxy(i)
498 ENDDO
499 END IF ! FISOKIN > ZERO
500C--------------------------------
501C HARDENING MODULUS
502C--------------------------------
503 hardm(1:nel) = h(1:nel)
504!-----------
505 RETURN
506 END
alpha
#define alpha
Definition eval.h:35
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:274
m2cplr
subroutine m2cplr(nel, ezz, off, pla, ipla, temp, c3, c4, m_exp, iform, etse, gs, epsp, vp, israte, yld, g, a1, a2, nu, ca0, cb0, cn, ymax0, epchk, young, cc, epdr, icc, dpla, tstar, fisokin, gama_imp, signor, hardm, depsxx, depsyy, depsxy, depsyz, depszx, signxx, signyy, signxy, signyz, signzx, sigbakxx, sigbakyy, sigbakxy, sigoxx, sigoyy, sigoxy, sigoyz, sigozx)
Definition m2cplr.F:39
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21