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!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
30!||====================================================================
31 SUBROUTINE m2lawpi(ELBUF_STR,MAT_PARAM,
32 1 NEL ,NPT ,GEO ,EINT ,OFF ,
33 3 PID ,EPSP ,EXX ,EXY ,EXZ ,
34 4 KXX ,KYY ,KZZ ,AL ,ASRATE ,
35 5 TIMESTEP,JTHE ,TEMP ,SIGY )
36C-----------------------------------------------
37C M o d u l e s
38C-----------------------------------------------
39 USE elbufdef_mod
40 use matparam_def_mod
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45#include "comlock.inc"
46C-----------------------------------------------
47C C o m m o n B l o c k s
48C-----------------------------------------------
49#include "scr17_c.inc"
50#include "param_c.inc"
51#include "com08_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER,INTENT(IN) :: NEL,NPT
56 INTEGER,INTENT(IN) :: JTHE
57 INTEGER PID(NEL)
58 my_real ,INTENT(IN) :: TIMESTEP
59 my_real ,INTENT(IN) :: asrate
60 my_real
61 . exx(nel), exy(nel), exz(nel),
62 . kxx(nel), kyy(nel), kzz(nel),
63 . geo(npropg,*), eint(nel,2),
64 . off(*),epsp(*), al(*)
65 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: temp
66 my_real ,INTENT(INOUT) :: sigy(nel,npt)
67 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
68 type (matparam_struct_) ,intent(in) :: mat_param
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72 INTEGER :: I,IPT,ICC,iform,IPY,IPZ,IPA,VP,IR,IS,ILAYER
73 INTEGER :: II(3)
74 my_real :: YOUNG,G,SHFACT,EPMX
75 my_real :: TREF,TINI,TMELT,CC,CP,EPDR,SVM1,GS,MT,TSTAR,UMR,R,PLAP1,Z4
76 my_real :: ca(nel), cb(nel), cn(nel),ymax(nel),
77 . z3(nel),yld(nel),etse(nel),q(nel),
78 . ypt(nel),zpt(nel),apt(nel),vol(nel),dpla(nel),
79 . signxx(nel),signxy(nel),signxz(nel),logep(nel),
80 . depsxx(nel),depsxy(nel),depsxz(nel),plap(nel)
81C
82 TYPE(l_bufel_) ,POINTER :: LBUF
83 TYPE(BUF_LAY_) ,POINTER :: BUFLY
84C=======================================================================
85 IPY = 200
86 ipz = 300
87 ipa = 400
88 shfact = five_over_6
89!
90 DO i=1,3
91 ii(i) = nel*(i-1)
92 ENDDO
93!
94C---
95 iform = mat_param%iparam(1)
96 icc = mat_param%iparam(2)
97 vp = mat_param%iparam(3)
98!
99 cc = mat_param%uparam(6)
100 epdr = mat_param%uparam(7)
101 if (vp == 1) then
102 epdr = max(epdr, em20)
103 else
104 epdr = max(epdr*dt1, em20)
105 endif
106 tini = mat_param%therm%tini
107 tref = mat_param%therm%tref
108 tmelt = mat_param%therm%tmelt
109 cp = mat_param%therm%rhocp
110 if (cp > zero) cp = one/cp
111
112 IF (iform == 1) THEN ! zerilli
113 z4 = mat_param%uparam(11)
114 ELSE
115 z4 = zero
116 END IF
117!
118 young = mat_param%young
119 g = mat_param%shear
120 gs = shfact*g
121 epmx = mat_param%uparam(4)
122!
123 DO i=1,nel
124 ca(i) = mat_param%uparam(1)
125 cb(i) = mat_param%uparam(2)
126 cn(i) = mat_param%uparam(3)
127 ymax(i) = mat_param%uparam(5)
128 z3(i) = mat_param%uparam(10)
129 epsp(i) = max(epsp(i),epdr)
130 vol(i) = al(i)*geo(1,pid(i))
131 ENDDO
132C
133 IF (jthe == 0 .and. cp > zero) THEN
134 DO i=1,nel
135 temp(i) = tini + (eint(i,1)+eint(i,2)) * cp / vol(i)
136 ENDDO
137 END IF
138C-------------------------------------
139C start of loop over integration points
140C--------------------------------------
141 ilayer= 1
142 ir = 1
143 is = 1
144 DO ipt= 1,npt
145
146 lbuf => elbuf_str%BUFLY(ilayer)%LBUF(ir,is,ipt)
147 bufly => elbuf_str%BUFLY(ilayer)
148C--- Coordonnees du point d'integration
149 DO i=1,nel
150 ypt(i) = geo(ipy+ipt,pid(i))
151 zpt(i) = geo(ipz+ipt,pid(i))
152 apt(i) = geo(ipa+ipt,pid(i))
153 ENDDO
154 DO i=1,nel
155 signxx(i) = lbuf%SIG(ii(1)+i)
156 signxy(i) = lbuf%SIG(ii(2)+i)
157 signxz(i) = lbuf%SIG(ii(3)+i)
158 ENDDO
159C--- Deformations Incrementales
160 DO i= 1,nel
161 depsxx(i) = exx(i) - ypt(i)*kzz(i) + zpt(i)*kyy(i)
162 depsxy(i) = exy(i) + zpt(i)*kxx(i)
163 depsxz(i) = exz(i) - ypt(i)*kxx(i)
164 depsxy(i) = depsxy(i) / shfact
165 depsxz(i) = depsxz(i) / shfact
166 ENDDO
167C
168c--- Total strain
169C
170 IF (bufly%L_STRA > 0) THEN
171 DO i= 1,nel
172 lbuf%STRA(ii(1)+i) = lbuf%STRA(ii(1)+i) + depsxx(i)
173 lbuf%STRA(ii(2)+i) = lbuf%STRA(ii(2)+i) + depsxy(i)
174 lbuf%STRA(ii(3)+i) = lbuf%STRA(ii(3)+i) + depsxz(i)
175 ENDDO
176 ENDIF
177C
178C--- Elastic constraints
179C
180 DO i = 1,nel
181 signxx(i) = signxx(i) + young*depsxx(i)
182 signxy(i) = signxy(i) + gs*depsxy(i)
183 signxz(i) = signxz(i) + gs*depsxz(i)
184 etse(i) = one
185 ENDDO
186c--- strain rate dependency
187 IF(vp == 1)THEN
188 DO i= 1,nel
189 plap(i) = bufly%MAT(ir,is,ipt)%VAR(i)
190 plap(i) = max(plap(i),epdr )
191 logep(i) = log(plap(i)/epdr)
192 ENDDO
193 ELSE
194 DO i= 1,nel
195 logep(i) = log(epsp(i)/epdr)
196 ENDDO
197 ENDIF
198C-- Yield
199 IF (cc /= zero) THEN
200 IF (iform == 0)THEN
201 DO i=1,nel
202 mt = max(em20,z3(i))
203 tstar = max(zero,(temp(i)-tref) / (tmelt-tref))
204 IF (tstar == zero) THEN
205 q(i) = (one + cc * logep(i))
206 ELSE
207 q(i) = (one + cc * logep(i))*(one-exp(mt*log(tstar)))
208 ENDIF
209 q(i) = max(q(i),em20)
210 ca(i) = ca(i) * q(i)
211 cb(i) = cb(i) * q(i)
212 IF (icc== 1) ymax(i) = ymax(i) * q(i)
213 ENDDO
214 ELSEIF (iform == 1) THEN
215 DO i=1,nel
216 q(i) = logep(i)
217 q(i) = cc*exp((-z3(i)+z4 * q(i))*temp(i))
218 IF (icc == 1) ymax(i)= ymax(i) + q(i)
219 ca(i) = ca(i) + q(i)
220 ENDDO
221 ENDIF
222 ELSEIF (iform == 0) THEN
223 DO i=1,nel
224 mt = max(em20,z3(i))
225 tstar = max(zero,(temp(i)-tref) / (tmelt-tref))
226 q(i) = one - tstar**mt
227 ENDDO
228 IF (icc==1) ymax(1:nel) = ymax(1:nel)*q(1:nel)
229 ENDIF
230C---
231 DO i=1,nel
232 IF(lbuf%PLA(i) == zero) THEN
233 yld(i)= ca(i)
234 ELSE
235 yld(i)= ca(i) + cb(i)*exp(cn(i)*log(lbuf%PLA(i)))
236 ENDIF
237 yld(i) = min(yld(i),ymax(i))
238 sigy(i,ipt) = yld(i)
239 ENDDO
240C-------------------
241C PROJECTION - radial return
242C-------------------
243 DO i=1,nel
244 svm1 = signxx(i)**2 + three*(signxy(i)**2 + signxz(i)**2)
245 IF (svm1 > yld(i)**2) THEN
246 svm1 = sqrt(svm1)
247 r = min( one, yld(i)/max(em20,svm1) )
248 signxx(i) = signxx(i)*r
249 signxy(i) = signxy(i)*r
250 signxz(i) = signxz(i)*r
251 umr = one - r
252 dpla(i) = svm1*umr/young
253 lbuf%PLA(i) = lbuf%PLA(i) + off(i)* dpla(i)
254c IF (R < ONE) ETSE(I)= H(I)/(H(I)+YOUNG)
255 ENDIF
256 ENDDO
257 IF (vp == 1) THEN
258 DO i=1,nel
259 plap1 = dpla(i) / max(em20,timestep)
260 bufly%MAT(ir,is,ipt)%VAR(i) = asrate * plap1 + (one - asrate) * plap(i)
261 ENDDO
262 ENDIF
263
264C failure criteria
265C--------------------------------
266C ductile failure test
267C--------------------------------
268C
269 DO i=1,nel
270 IF (lbuf%PLA(i) >= epmx .AND. off(i) == one) THEN
271 off(i)=four_over_5
272 ENDIF
273 ENDDO
274C-----------
275C
276 DO i=1,nel
277 lbuf%SIG(ii(1)+i) = signxx(i)
278 lbuf%SIG(ii(2)+i) = signxy(i)
279 lbuf%SIG(ii(3)+i) = signxz(i)
280 ENDDO
281C-------------------------------------
282C end of loop over integration point
283C-------------------------------------
284 ENDDO
285C-----------
286 RETURN
287 END
288
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:274
m2lawpi
subroutine m2lawpi(elbuf_str, mat_param, nel, npt, geo, eint, off, pid, epsp, exx, exy, exz, kxx, kyy, kzz, al, asrate, timestep, jthe, temp, sigy)
Definition m2lawpi.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21