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