OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2lawpi.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!|| m2lawpi ../engine/source/materials/mat/mat002/m2lawpi.F
25!||--- called by ------------------------------------------------------
26!|| main_beam18 ../engine/source/elements/beam/main_beam18.f
27!||--- uses -----------------------------------------------------
28!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
29!||====================================================================
30 SUBROUTINE m2lawpi(ELBUF_STR,
31 1 JFT ,JLT ,NPT ,PM ,GEO,
32 2 EINT ,OFF ,IMAT,
33 3 PID ,EPSP ,EXX ,EXY ,EXZ,
34 4 KXX ,KYY ,KZZ ,AL ,NEL,
35 5 IPM ,ASRATE ,TIMESTEP,NUMMAT,SIGY)
36C-----------------------------------------------
37C M o d u l e s
38C-----------------------------------------------
39 USE elbufdef_mod
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44#include "comlock.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "scr17_c.inc"
49#include "param_c.inc"
50#include "com08_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER ,INTENT(IN) :: IMAT,NUMMAT
55 INTEGER ,DIMENSION(NPROPMI,NUMMAT) ,INTENT(IN) :: IPM
56 INTEGER JFT,JLT,NPT,NEL
57 INTEGER PID(NEL)
58 my_real , INTENT(IN) :: TIMESTEP
59 my_real , INTENT(IN) :: asrate
60C REAL
61 my_real
62 . exx(nel), exy(nel), exz(nel),
63 . kxx(nel), kyy(nel), kzz(nel),
64 . pm(npropm,*), geo(npropg,*), eint(nel,2),
65 . off(*),epsp(*), al(*)
66 my_real , INTENT(INOUT)::sigy(nel,npt)
67C
68 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER :: I,IPT,ICC,IRTY,IPY,IPZ,IPA,VP,IR,IS,ILAYER
73 INTEGER :: II(3)
74 my_real :: EPMX(NEL),CA(NEL), CB(NEL), CN(NEL), YMAX(NEL),
75 . T(NEL),Z3(NEL),CC(NEL),EPDR(NEL),
76 . YLD(NEL),ETSE(NEL),Q(NEL),E(NEL),G(NEL),
77 . ypt(nel),zpt(nel),apt(nel),vol(nel),dpla(nel),
78 . signxx(nel),signxy(nel),signxz(nel),logep(nel),
79 . depsxx(nel),depsxy(nel),depsxz(nel),plap(nel)
80 my_real :: tref,tmelt,cp,
81 . epif,svm1,gs,mt,tstar,umr,r,shfact,plap1,z4
82C
83 TYPE(l_bufel_) ,POINTER :: LBUF
84 TYPE(BUF_LAY_) ,POINTER :: BUFLY
85C=======================================================================
86 ipy = 200
87 ipz = 300
88 ipa = 400
89 epif = zero
90 shfact = five_over_6
91!
92 z4 = zero
93 DO i=1,3
94 ii(i) = nel*(i-1)
95 ENDDO
96!
97C---
98 vp = ipm(255,imat)
99 irty = nint(pm(50,imat))
100 IF (irty == 0) THEN
101 tref = pm(79,imat)
102 tmelt = pm(80,imat)
103 cp = one / max(em20,pm(69,imat))
104 ELSE
105 z4 = pm(52,imat)
106 cp = pm(53,imat)
107 tref = pm(54,imat)
108 tmelt = 0
109 END IF
110!
111 icc = nint(pm(49,imat))
112 DO i=jft,jlt
113 e(i) = pm(20,imat)
114 g(i) = pm(22,imat)
115 ca(i) = pm(38,imat)
116 cb(i) = pm(39,imat)
117 cn(i) = pm(40,imat)
118 epmx(i) = pm(41,imat)
119 ymax(i) = pm(42,imat)
120 cc(i) = pm(43,imat)
121 IF(vp == 1)THEN
122 epdr(i) = max(em20,pm(44,imat))
123 ELSE
124 epdr(i) = max(em20,pm(44,imat)*dt1)
125 ENDIF
126 epif = max(epif,cc(i))
127 z3(i) = pm(51,imat)
128 epsp(i) = max(epsp(i),epdr(i))
129 vol(i) = al(i)*geo(1,pid(i))
130 ENDDO
131C
132 DO i=jft,jlt
133 t(i) = tref + cp*(eint(i,1)+eint(i,2))/vol(i)
134 ENDDO
135C-------------------------------------
136C DEBUT DE BOUCLE SUR POINTS INTEGRATION
137C--------------------------------------
138C
139 DO ipt= 1,npt
140
141 ilayer=1
142 ir = 1
143 is = 1
144 lbuf => elbuf_str%BUFLY(ilayer)%LBUF(ir,is,ipt)
145 bufly => elbuf_str%BUFLY(ilayer)
146C--- Coordonnees du point d'integration
147 DO i=jft,jlt
148 ypt(i) = geo(ipy+ipt,pid(i))
149 zpt(i) = geo(ipz+ipt,pid(i))
150 apt(i) = geo(ipa+ipt,pid(i))
151 ENDDO
152 DO i=jft,jlt
153 signxx(i) = lbuf%SIG(ii(1)+i)
154 signxy(i) = lbuf%SIG(ii(2)+i)
155 signxz(i) = lbuf%SIG(ii(3)+i)
156 ENDDO
157C--- Deformations Incrementales
158 DO i= jft,jlt
159 depsxx(i) = exx(i) - ypt(i)*kzz(i) + zpt(i)*kyy(i)
160 depsxy(i) = exy(i) + zpt(i)*kxx(i)
161 depsxz(i) = exz(i) - ypt(i)*kxx(i)
162 depsxy(i) = depsxy(i) / shfact
163 depsxz(i) = depsxz(i) / shfact
164 ENDDO
165C
166c--- Total strain
167C
168 IF (bufly%L_STRA > 0) THEN
169 DO i= jft,jlt
170 lbuf%STRA(ii(1)+i) = lbuf%STRA(ii(1)+i) + depsxx(i)
171 lbuf%STRA(ii(2)+i) = lbuf%STRA(ii(2)+i) + depsxy(i)
172 lbuf%STRA(ii(3)+i) = lbuf%STRA(ii(3)+i) + depsxz(i)
173 ENDDO
174 ENDIF
175C
176C--- Contraintes elastiques
177C
178 DO i = jft,jlt
179 gs = shfact*g(i)
180 signxx(i) = signxx(i) + e(i)*depsxx(i)
181 signxy(i) = signxy(i) + gs*depsxy(i)
182 signxz(i) = signxz(i) + gs*depsxz(i)
183 etse(i) = one
184 ENDDO
185 DO i=jft,jlt
186 etse(i) = one
187 ca(i) = pm(38,imat)
188 cb(i) = pm(39,imat)
189 ymax(i) = pm(42,imat)
190 ENDDO
191c--- strain rate dependency
192 IF(vp == 1)THEN
193 DO i= jft,jlt
194 plap(i) = bufly%MAT(ir,is,ipt)%VAR(i)
195 plap(i) = max(plap(i),epdr(i) )
196 logep(i) = log(plap(i)/epdr(i))
197 ENDDO
198 ELSE
199 DO i= jft,jlt
200 logep(i) = log(epsp(i)/epdr(i))
201 ENDDO
202 ENDIF
203C-- Yield
204 IF (epif /= zero) THEN
205 IF (irty == 0)THEN
206 DO i=jft,jlt
207 mt = max(em20,z3(i))
208 tstar = max(zero,(t(i)-tref)/(tmelt-tref))
209 IF (tstar == zero) THEN
210 q(i) = (one + cc(i) * logep(i))
211 ELSE
212 q(i) = (one + cc(i) * logep(i))*(one-exp(mt*log(tstar)))
213 ENDIF
214 q(i) = max(q(i),em20)
215 ca(i) = ca(i) * q(i)
216 cb(i) = cb(i) * q(i)
217 IF (icc== 1) ymax(i) = ymax(i) * q(i)
218 ENDDO
219 ELSEIF (irty == 1) THEN
220 DO i=jft,jlt
221 q(i) = logep(i)
222 q(i) = cc(i)*exp((-z3(i)+z4 * q(i))*t(i))
223 IF (icc == 1) ymax(i)= ymax(i) + q(i)
224 ca(i) = ca(i) + q(i)
225 ENDDO
226 ENDIF
227 ENDIF
228C---
229 DO i=jft,jlt
230 IF(lbuf%PLA(i) == zero) THEN
231 yld(i)= ca(i)
232 ELSE
233 yld(i)= ca(i) + cb(i)*exp(cn(i)*log(lbuf%PLA(i)))
234 ENDIF
235 yld(i) = min(yld(i),ymax(i))
236 sigy(i,ipt) = yld(i)
237 ENDDO
238C-------------------
239C PROJECTION - radial return
240C-------------------
241 DO i=jft,jlt
242 svm1 = signxx(i)**2 + three*(signxy(i)**2 + signxz(i)**2)
243 IF (svm1 > yld(i)**2) THEN
244 svm1 = sqrt(svm1)
245 r = min( one, yld(i)/max(em20,svm1) )
246 signxx(i) = signxx(i)*r
247 signxy(i) = signxy(i)*r
248 signxz(i) = signxz(i)*r
249 umr = one - r
250 dpla(i) = svm1*umr/(e(i))
251 lbuf%PLA(i) = lbuf%PLA(i) + off(i)* dpla(i)
252c IF (R < ONE) ETSE(I)= H(I)/(H(I)+E(I))
253 ENDIF
254 ENDDO
255 IF (vp == 1) THEN
256 DO i=jft,jlt
257 plap1 = dpla(i) / max(em20,timestep)
258 bufly%MAT(ir,is,ipt)%VAR(i) = asrate * plap1 + (one - asrate) * plap(i)
259 ENDDO
260 ENDIF
261
262
263C failure criteria
264C--------------------------------
265C TEST DE RUPTURE DUCTILE
266C--------------------------------
267C
268 DO i=jft,jlt
269 IF (lbuf%PLA(i) >= epmx(i) .AND. off(i) == one) THEN
270 off(i)=four_over_5
271 ENDIF
272 ENDDO
273C-----------
274C
275 DO i=jft,jlt
276 lbuf%SIG(ii(1)+i) = signxx(i)
277 lbuf%SIG(ii(2)+i) = signxy(i)
278 lbuf%SIG(ii(3)+i) = signxz(i)
279 ENDDO
280C-------------------------------------
281C FIN DE BOUCLE SUR POINT INTEGRATION
282C-------------------------------------
283 ENDDO
284C-----------
285 RETURN
286 END
287
m2lawpi
subroutine m2lawpi(elbuf_str, jft, jlt, npt, pm, geo, eint, off, imat, pid, epsp, exx, exy, exz, kxx, kyy, kzz, al, nel, ipm, asrate, timestep, nummat, sigy)
Definition m2lawpi.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
main_beam18
subroutine main_beam18(elbuf_str, nel, npt, mtn, imat, pid, ngl, pm, ipm, geo, off, for, mom, eint, al, epsd, bufmat, npf, tf, exx, exy, exz, kxx, kyy, kzz, f1, f2, f3, m1, m2, m3, jthe, tempel, ifail, sbufmat, snpc, stf, nummat, numgeo, iout, istdo, npropmi, npropm, npropg, time, dtime, idel7nok, isigi, imconv, ismstr, mat_param, ntable, table)
Definition main_beam18.F:48