OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m22law.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!|| m22law ../engine/source/materials/mat/mat022/m22law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| mdtsph ../engine/source/materials/mat_share/mdtsph.F
29!|| mqviscb ../engine/source/materials/mat_share/mqviscb.F
30!||--- uses -----------------------------------------------------
31!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
32!|| prop_param_mod ../common_source/modules/mat_elem/prop_param_mod.F90
33!||====================================================================
34 SUBROUTINE m22law(
35 1 PM, OFF, SIG, EINT,
36 2 RHO, QOLD, EPXE, EPD,
37 3 VOL, PDAM, STIFN, DT2T,
38 4 NELTST, ITYPTST, OFFG, GEO,
39 5 PID, AMU, VOL_AVG, MUMAX,
40 6 MAT, NGL, SSP, DVOL,
41 7 AIRE, VNEW, VD2, DELTAX,
42 8 VIS, D1, D2, D3,
43 9 D4, D5, D6, PNEW,
44 A PSH, QNEW, SSP_EQ, SOLD1,
45 B SOLD2, SOLD3, SOLD4, SOLD5,
46 C SOLD6, IPLA, SIGY, DEFP,
47 D DPLA, MSSA, DMELS, CONDE,
48 E DTEL, G_DT, NEL, IPM,
49 F RHOREF, RHOSP, NFT, JSPH,
50 G ITY, JTUR, JTHE, ISMSTR,
51 H JSMS, NPG , glob_therm,
52 I NUMGEO, IGEO)
53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 use glob_therm_mod
57 USE prop_param_mod , only : n_var_igeo
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62#include "comlock.inc"
63C-----------------------------------------------
64C G l o b a l P a r a m e t e r s
65C-----------------------------------------------
66#include "mvsiz_p.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "com08_c.inc"
71#include "param_c.inc"
72#include "scr03_c.inc"
73#include "scr17_c.inc"
74#include "units_c.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER, INTENT(IN) :: ISMSTR
79 INTEGER, INTENT(IN) :: JSMS
80 INTEGER, INTENT(IN) :: ITY
81 INTEGER, INTENT(IN) :: JTUR
82 INTEGER, INTENT(IN) :: JTHE
83 INTEGER, INTENT(IN) :: NFT
84 INTEGER, INTENT(IN) :: JSPH,NPG
85 INTEGER,INTENT(IN) :: NUMGEO
86 INTEGER NELTST,ITYPTST ,PID(*), G_DT
87 INTEGER MAT(*),NGL(*),IPLA,NEL,IPM(NPROPMI,*)
88 my_real
89 . DT2T
90 my_real
91 . PDAM
92 my_real
93 . PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), RHO(*), QOLD(*),
94 . EPXE(*), EPD(*), VOL(*),STIFN(*), OFFG(*),GEO(NPROPG,*) ,
95 . MUMAX(*), SIGY(*), DEFP(*), DPLA(MVSIZ),
96 . amu(*), vol_avg(*)
97 my_real
98 . vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*),
99 . psh(*), pnew(1), qnew(*) ,ssp_eq(*),
100 . d1(*), d2(*), d3(*), d4(*), d5(*), d6(*),
101 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz),
102 . sold4(mvsiz), sold5(mvsiz), sold6(mvsiz),
103 . mssa(*), dmels(*), conde(*),dtel(*),
104 . rhoref(*) ,rhosp(*)
105 type (glob_therm_) ,intent(inout) :: glob_therm
106 integer,dimension(n_var_igeo,numgeo),intent(in) :: igeo
107C-----------------------------------------------
108C L o c a l V a r i a b l e s
109C-----------------------------------------------
110 INTEGER ICC(MVSIZ), I, MX, II,ICC_1
111 my_real
112 . RHO0(MVSIZ),
113 . G(MVSIZ),AK(MVSIZ),
114 . QH(MVSIZ), C1(MVSIZ),
115 . p(mvsiz), epmx(mvsiz),
116 . ca(mvsiz), cb(mvsiz), cc(mvsiz), cn(mvsiz),
117 . epdr(mvsiz),
118 . dvol(mvsiz), sigmx(mvsiz),
119 . ce(mvsiz),epsl(mvsiz), ql(mvsiz), yldl(mvsiz), aj2(mvsiz),
120 . g0(mvsiz),
121 . e1, e2, e3, e4, e5, e6, g1, g2,
122 . epsp, hl, depsl, alpe, dav, scale, bid1, bid2,
123 . bid3, einc, dta,facq0,
124 . rho0_1,c1_1,ca_1,cb_1,cn_1,
125 . epmx_1,sigmx_1,cc_1,epdr_1,epsl_1,
126 . yldl_1,ql_1
127
128 facq0 = one
129 IF(ipla==0)THEN
130 mx = mat(1)
131 rho0_1 =pm( 1,mx)
132 c1_1 =pm(32,mx)
133 ca_1 =pm(38,mx)
134 cb_1 =pm(39,mx)
135 cn_1 =pm(40,mx)
136 epmx_1 =pm(41,mx)
137 sigmx_1=pm(42,mx)
138 cc_1 =pm(43,mx)
139 epdr_1 =pm(44,mx)
140 epsl_1 =pm(45,mx)
141 yldl_1 =pm(47,mx)
142 ql_1 =pm(48,mx)
143 icc_1 =nint(pm(49,mx))
144 DO 10 i=1,nel
145 rho0(i) =rho0_1
146 g0(i) =pm(22,mx)*off(i)
147 c1(i) =c1_1
148 ca(i) =ca_1
149 cb(i) =cb_1
150 cn(i) =cn_1
151 epmx(i) =epmx_1
152 sigmx(i)=sigmx_1
153 cc(i) =cc_1
154 epdr(i) =epdr_1
155 epsl(i) =epsl_1
156 yldl(i) =yldl_1
157 ql(i) =ql_1
158 icc(i) =icc_1
159 10 CONTINUE
160 ELSE
161
162 mx = mat(1)
163 rho0_1 =pm( 1,mx)
164 c1_1 =pm(32,mx)
165 ca_1 =pm(38,mx)
166 cb_1 =pm(39,mx)
167 cn_1 =pm(40,mx)
168 epmx_1 =pm(41,mx)
169 sigmx_1=pm(42,mx)
170 cc_1 =pm(43,mx)
171 epdr_1 =pm(44,mx)
172 epsl_1 =pm(45,mx)
173 ql_1 =pm(46,mx)
174 yldl_1 =pm(47,mx)
175 icc_1 =nint(pm(49,mx))
176 DO 11 i=1,nel
177 rho0(i) =rho0_1
178 g0(i) =pm(22,mx)*off(i)
179 c1(i) =c1_1
180 ca(i) =ca_1
181 cb(i) =cb_1
182 cn(i) =cn_1
183 epmx(i) =epmx_1
184 sigmx(i)=sigmx_1
185 cc(i) =cc_1
186 epdr(i) =epdr_1
187 epsl(i) =epsl_1
188 ql(i) =ql_1
189 yldl(i) =yldl_1
190 icc(i) =icc_1
191 11 CONTINUE
192 ENDIF
193
194 DO 15 i=1,nel
195 epd(i)=off(i) * max( abs(d1(i)), abs(d2(i)), abs(d3(i)), half*abs(d4(i)),half*abs(d5(i)),half*abs(d6(i)))
196 epsp = max(epd(i),epdr(i))
197 ce(i) = one +cc(i) * log(epsp/epdr(i))
198 15 CONTINUE
199
200 IF(ipla/=2)THEN
201 DO 20 i=1,nel
202 IF(cn(i)==one) THEN
203 ak(i)= ca(i)+cb(i)*epxe(i)
204 qh(i)= cb(i)*ce(i)
205 ELSE
206 IF(epxe(i)>zero) THEN
207 ak(i)=ca(i)+cb(i)*epxe(i)**cn(i)
208 IF(cn(i)>one) THEN
209 qh(i)= (cb(i)*cn(i)*epxe(i)**(cn(i)-one))*ce(i)
210 ELSE
211 qh(i)= (cb(i)*cn(i)/epxe(i)**(one - cn(i)))*ce(i)
212 ENDIF
213 ELSE
214 ak(i)=ca(i)
215 qh(i)=zero
216 ENDIF
217 ENDIF
218 sigy(i) = ak(i)
219 IF(epxe(i)>epmx(i))ak(i)=zero
220 IF(epxe(i)>epsl(i))qh(i)=ql(i)
221 20 CONTINUE
222 ELSE
223 DO i=1,nel
224 IF(cn(i)==one) THEN
225 ak(i)= ca(i)+cb(i)*epxe(i)
226 ELSE
227 IF(epxe(i)>zero) THEN
228 ak(i)=ca(i)+cb(i)*epxe(i)**cn(i)
229 ELSE
230 ak(i)=ca(i)
231 ENDIF
232 ENDIF
233 sigy(i) = ak(i)
234 IF(epxe(i)>epmx(i))ak(i)=zero
235 ENDDO
236 ENDIF
237
238 IF(ipla==0)THEN
239 DO i=1,nel
240 ak(i) = min(ak(i),sigmx(i))
241 hl = three*g0(i)*ql(i)/(three*g0(i)+ql(i))
242 depsl = max(zero,epxe(i)-epsl(i))
243 ak(i) = min(ak(i),yldl(i)+hl*depsl)
244 ak(i) = max(ak(i),zero)
245 alpe = min(one,ak(i)/max(ak(i)+three*g0(i)*depsl,em15))
246 ak(i) = ak(i)*ce(i)
247 IF(icc(i)==2) ak(i) = min(ak(i),sigmx(i))
248 g(i) = alpe*g0(i)
249 c1(i) = pdam*alpe*c1(i) + (1.-pdam)*c1(i)
250 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
251 g1=dt1*g(i)
252 g2=two*g1
253 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
254 !-------------------------------
255 ! DEVIATORIC STRESSES
256 !-------------------------------
257 dav =-third*(d1(i)+d2(i)+d3(i))
258 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
259 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
260 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
261 sig(i,4)=sig(i,4)+g1*d4(i)
262 sig(i,5)=sig(i,5)+g1*d5(i)
263 sig(i,6)=sig(i,6)+g1*d6(i)
264 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)
265 1 +sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
266 aj2(i)=sqrt(three*aj2(i))
267 ENDDO
268 ELSEIF(ipla==2)THEN
269 DO i=1,nel
270 ak(i) = min(ak(i),sigmx(i))
271 depsl = max(zero,epxe(i)-epsl(i))
272 ak(i) = min(ak(i),yldl(i)+ql(i)*depsl)
273 ak(i) = max(ak(i),zero)
274 alpe = min(one,ak(i)/max(ak(i) + three*g0(i)*depsl,em15))
275 ak(i) = ak(i)*ce(i)
276 IF(icc(i)==2) ak(i) = min(ak(i),sigmx(i))
277 g(i) = alpe*g0(i)
278 c1(i) = pdam*alpe*c1(i) + (one - pdam)*c1(i)
279 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
280 g1=dt1*g(i)
281 g2=two*g1
282 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
283 !-------------------------------
284 ! DEVIATORIC STRESSES
285 !-------------------------------
286 IF(epxe(i)>epsl(i).AND.three*g(i)<em15)THEN
287 sig(i,1)=zero
288 sig(i,2)=zero
289 sig(i,3)=zero
290 sig(i,4)=zero
291 sig(i,5)=zero
292 sig(i,6)=zero
293 aj2(i) =zero
294 ELSE
295 dav =-third*(d1(i)+d2(i)+d3(i))
296 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
297 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
298 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
299 sig(i,4)=sig(i,4)+g1*d4(i)
300 sig(i,5)=sig(i,5)+g1*d5(i)
301 sig(i,6)=sig(i,6)+g1*d6(i)
302 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)+
303 . sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
304 aj2(i)=sqrt(three*aj2(i))
305 ENDIF
306 ENDDO
307 ELSE
308 DO i=1,nel
309 ak(i) = min(ak(i),sigmx(i))
310 depsl = max(zero,epxe(i)-epsl(i))
311 ak(i) = min(ak(i),yldl(i)+ql(i)*depsl)
312 ak(i) = max(ak(i),zero)
313 alpe = min(one,ak(i)/max(ak(i) + three*g0(i)*depsl,em15))
314 ak(i) = ak(i)*ce(i)
315 IF(icc(i)==2) ak(i) = min(ak(i),sigmx(i))
316 g(i) = alpe*g0(i)
317 c1(i) = pdam*alpe*c1(i) + (one - pdam)*c1(i)
318 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
319 g1=dt1*g(i)
320 g2=two*g1
321 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
322 !-------------------------------
323 ! DEVIATORIC STRESSES
324 !-------------------------------
325 IF(epxe(i)>epsl(i).AND.g(i)*(three*g0(i)+qh(i))<em15)THEN
326 sig(i,1)=zero
327 sig(i,2)=zero
328 sig(i,3)=zero
329 sig(i,4)=zero
330 sig(i,5)=zero
331 sig(i,6)=zero
332 aj2(i) =zero
333 ELSE
334 dav =-third*(d1(i)+d2(i)+d3(i))
335 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
336 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
337 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
338 sig(i,4)=sig(i,4)+g1*d4(i)
339 sig(i,5)=sig(i,5)+g1*d5(i)
340 sig(i,6)=sig(i,6)+g1*d6(i)
341 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)+
342 . sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
343 aj2(i)=sqrt(three*aj2(i))
344 ENDIF
345 ENDDO
346 ENDIF
347
348 IF(ipla==0)THEN
349 DO i=1,nel
350 scale= min(one,ak(i) / max(aj2(i),em15))
351 sig(i,1)=scale*sig(i,1)
352 sig(i,2)=scale*sig(i,2)
353 sig(i,3)=scale*sig(i,3)
354 sig(i,4)=scale*sig(i,4)
355 sig(i,5)=scale*sig(i,5)
356 sig(i,6)=scale*sig(i,6)
357 epxe(i)=epxe(i)+(one -scale)*aj2(i)/max(three*g(i)+qh(i),em15)
358 dpla(i) = (one -scale)*aj2(i)/max(three*g(i)+qh(i),em15)
359 ENDDO
360 ELSEIF(ipla==2)THEN
361 DO i=1,nel
362 scale= min(one,ak(i) / max(aj2(i),em15))
363 sig(i,1)=scale*sig(i,1)
364 sig(i,2)=scale*sig(i,2)
365 sig(i,3)=scale*sig(i,3)
366 sig(i,4)=scale*sig(i,4)
367 sig(i,5)=scale*sig(i,5)
368 sig(i,6)=scale*sig(i,6)
369 epxe(i)=epxe(i)+(one -scale)*aj2(i)/max(three*g(i),em15)
370 dpla(i) = (one -scale)*aj2(i)/max(three*g(i),em15)
371 ENDDO
372 ELSEIF(ipla==1)THEN
373 DO i=1,nel
374 scale= min(one,ak(i) / max(aj2(i),em15))
375 dpla(i)=(one -scale)*aj2(i)*g0(i)/ max(g(i)*(three*g0(i)+qh(i)),em15)
376 ak(i)=max(zero,ak(i)+dpla(i)*qh(i))
377 scale= min(one,ak(i)/ max(aj2(i),em15))
378 sig(i,1)=scale*sig(i,1)
379 sig(i,2)=scale*sig(i,2)
380 sig(i,3)=scale*sig(i,3)
381 sig(i,4)=scale*sig(i,4)
382 sig(i,5)=scale*sig(i,5)
383 sig(i,6)=scale*sig(i,6)
384 epxe(i)=epxe(i)+dpla(i)
385 ENDDO
386 ENDIF
387
388 IF (jsph==0)THEN
389 CALL mqviscb(
390 1 pm, off, rho, bid1,
391 2 bid2, ssp, bid3, stifn,
392 3 dt2t, neltst, ityptst, aire,
393 4 offg, geo, pid, vnew,
394 5 vd2, deltax, vis, d1,
395 6 d2, d3, pnew, psh,
396 7 mat, ngl, qnew, ssp_eq,
397 8 vol, mssa, dmels, igeo,
398 9 facq0, conde, dtel, g_dt,
399 a ipm, rhoref, rhosp, nel,
400 b ity, ismstr, jtur, jthe,
401 c jsms, npg , glob_therm)
402 ELSE
403 CALL mdtsph(
404 1 pm, off, rho, bid1,
405 2 bid2, bid3, stifn, dt2t,
406 3 neltst, ityptst, offg, geo,
407 4 pid, mumax, ssp, vnew,
408 5 vd2, deltax, vis, d1,
409 6 d2, d3, pnew, psh,
410 7 mat, ngl, qnew, ssp_eq,
411 8 g_dt, dtel, nel, ity,
412 9 jtur, jthe)
413 ENDIF
414
415 DO 120 i=1,nel
416 pnew(i)=c1(i)*amu(i)
417 120 CONTINUE
418
419C----------------------------
420C Ductile rupture test
421C---------------------------
422 IF(pdam==one.AND.codvers>=42)THEN
423 DO i=1,nel
424 IF(off(i)<em01) off(i)=zero
425 IF(off(i)<one) off(i)=off(i)*four_over_5
426 ENDDO
427C
428 DO i=1,nel
429 IF(epmx(i)/=zero.AND.off(i)>=one.AND. epxe(i)>epmx(i))THEN
430 off(i)=off(i)*four_over_5
431 ii=i+nft
432#include "lockon.inc"
433 WRITE(iout,1000) ngl(i)
434#include "lockoff.inc"
435 ENDIF
436 ENDDO
437 ENDIF
438
439 dta = half*dt1
440
441 DO i=1,nel
442 sig(i,1)=(sig(i,1)-pnew(i))*off(i)
443 sig(i,2)=(sig(i,2)-pnew(i))*off(i)
444 sig(i,3)=(sig(i,3)-pnew(i))*off(i)
445 sig(i,4)= sig(i,4)*off(i)
446 sig(i,5)= sig(i,5)*off(i)
447 sig(i,6)= sig(i,6)*off(i)
448 e1=d1(i)*(sold1(i)+sig(i,1))
449 e2=d2(i)*(sold2(i)+sig(i,2))
450 e3=d3(i)*(sold3(i)+sig(i,3))
451 e4=d4(i)*(sold4(i)+sig(i,4))
452 e5=d5(i)*(sold5(i)+sig(i,5))
453 e6=d6(i)*(sold6(i)+sig(i,6))
454 einc= vol_avg(i)*(e1+e2+e3+e4+e5+e6)*dta - half*dvol(i)*(qold(i)+qnew(i))
455 eint(i)=(eint(i)+einc*off(i)) / max(em15,vol(i))
456 qold(i)=qnew(i)
457 ENDDO
458
459 DO i=1,nel
460 defp(i)=epxe(i)
461 sigy(i)=max(sigy(i),ak(i))
462 ENDDO
463
464 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
465 RETURN
466 END
dtel
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
Definition dtel.F:46
m22law
subroutine m22law(pm, off, sig, eint, rho, qold, epxe, epd, vol, pdam, stifn, dt2t, neltst, ityptst, offg, geo, pid, amu, vol_avg, mumax, mat, ngl, ssp, dvol, aire, vnew, vd2, deltax, vis, d1, d2, d3, d4, d5, d6, pnew, psh, qnew, ssp_eq, sold1, sold2, sold3, sold4, sold5, sold6, ipla, sigy, defp, dpla, mssa, dmels, conde, dtel, g_dt, nel, ipm, rhoref, rhosp, nft, jsph, ity, jtur, jthe, ismstr, jsms, npg, glob_therm, numgeo, igeo)
Definition m22law.F:53
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mdtsph
subroutine mdtsph(pm, off, rho, rk, t, re, sti, dt2t, neltst, ityptst, offg, geo, pid, mumax, ssp, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, g_dt, dtsph, nel, ity, jtur, jthe)
Definition mdtsph.F:47
mqviscb
subroutine mqviscb(pm, off, rho, rk, temp, ssp, re, sti, dt2t, neltst, ityptst, aire, offg, geo, pid, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, vol0, mssa, dmels, igeo, facq0, conde, dtel, g_dt, ipm, rhoref, rhosp, nel, ity, ismstr, jtur, jthe, jsms, npg, glob_therm)
Definition mqviscb.F:56