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