OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2law8.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!|| m2law8 ../engine/source/materials/mat/mat002/m2law8.F
25!||--- called by ------------------------------------------------------
26!|| mmain8 ../engine/source/materials/mat_share/mmain8.F
27!||--- calls -----------------------------------------------------
28!|| mqvisc8 ../engine/source/materials/mat_share/mqvisc8.F
29!||--- uses -----------------------------------------------------
30!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
31!||====================================================================
32 SUBROUTINE m2law8(
33 1 PM, OFF, SIG, EINT,
34 2 RHO, QOLD, VOL, STIFN,
35 3 NEL, D1, D2, D3,
36 4 D4, D5, D6, VNEW,
37 5 VOLGP, DELTAX, RHO0, DVOL,
38 6 VD2, VIS, MAT, NC,
39 7 NGL, GEO, PID, EPSEQ,
40 8 DT2T, NELTST, ITYPTST, IPLA,
41 9 OFFG, DPLA, EPSP, TSTAR,
42 A ETSE, MSSA, DMELS, BUFLY,
43 B SSP, ITY, NPT, JTUR,
44 C JTHE, JSMS)
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE elbufdef_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "com08_c.inc"
61#include "param_c.inc"
62#include "impl1_c.inc"
63C-----------------------------------------------
64C D u m m y A r g u m e n t s
65C-----------------------------------------------
66 INTEGER, INTENT(IN) :: ITY
67 INTEGER, INTENT(IN) :: NPT
68 INTEGER, INTENT(IN) :: JTUR
69 INTEGER, INTENT(IN) :: JTHE
70 INTEGER, INTENT(IN) :: JSMS
71 INTEGER MAT(MVSIZ),NC(8,MVSIZ),NGL(MVSIZ),PID(MVSIZ)
72 INTEGER NEL,NELTST,ITYPTST,IPLA
73C REAL
74 my_real
75 . PM(NPROPM,*),OFF(MVSIZ),SIG(NEL,6),EINT(NEL),DELTAX(MVSIZ),
76 . RHO(NEL),QOLD(NEL),VOL(NEL) ,STIFN(*),EPSEQ(*),
77 . D1(MVSIZ,*) , D2(MVSIZ,*) ,
78 . d3(mvsiz,*) , d4(mvsiz,*) ,
79 . d5(mvsiz,*) , d6(mvsiz,*) ,
80 . vnew(mvsiz), rho0(mvsiz), dvol(mvsiz), volgp(mvsiz,*),
81 . vd2(mvsiz) , vis(mvsiz),geo(npropg,*), dt2t, offg(*),
82 . dpla(*),epsp(*),tstar(*),etse(*), mssa(*), dmels(*), ssp(mvsiz)
83 TYPE (BUF_LAY_), TARGET :: BUFLY
84C-----------------------------------------------
85C L o c a l V a r i a b l e s
86C-----------------------------------------------
87 INTEGER I,J,II,IPT,JPT,NPIF,MX,JJ(6)
88 INTEGER ICC,IRTY
89C REAL
90 my_real
91 . SOLD1(MVSIZ), SOLD2(MVSIZ), SOLD3(MVSIZ),
92 . SOLD4(MVSIZ), SOLD5(MVSIZ), SOLD6(MVSIZ),
93 . G(MVSIZ) , C1 , P(MVSIZ) ,
94 . G1(MVSIZ) , G2(MVSIZ), AJ2(MVSIZ), QH(MVSIZ),
95 . df(mvsiz) , amu(mvsiz) , einc(mvsiz), epd(mvsiz) ,
96 . dpdm(mvsiz), pnew(mvsiz) , ak(mvsiz),
97 . ca(mvsiz),cc, cn, epxo(mvsiz),
98 . epmx, thetl(mvsiz), epdr, t(mvsiz),
99 . sigmx(mvsiz),
100 . epif,tm,mt, scale, dta, dav,rhocpi,
101 . ca_1,cb_1,sigmx_1,cmx_1,z3_1,z4_1,cp_1,t_1
102 my_real,
103 . DIMENSION(:), POINTER :: sigp, epla
104 TYPE(l_bufel_) ,POINTER :: LBUF
105C=======================================================================
106 EPIF = zero
107 npif = 0
108C
109 dta =half*dt1
110 mx=mat(1)
111 icc =nint(pm(49,mx))
112 c1 =pm(32,mx)
113 cn =pm(40,mx)
114 epmx =pm(41,mx)
115 epdr =pm(44,mx)
116 cc =pm(43,mx)
117 epif = max(epif,cc)
118 irty = nint(pm(50,mx)) ! flag : J-C = 0, Zerilli = 1
119C
120 ca_1 =pm(38,mx)
121 cb_1 =pm(39,mx)
122 sigmx_1=pm(42,mx)
123 cmx_1 =pm(45,mx)
124 z3_1 =pm(51,mx)
125 z4_1 =pm(52,mx)
126 cp_1 = pm(53,mx)
127 rhocpi = pm(69,mx)
128 IF (rhocpi > zero) rhocpi = one / rhocpi
129 t_1 = pm(79,mx)
130 tm = pm(80,mx)
131C
132 DO i=1,nel
133 g(i) =pm(22,mx)*off(i)
134 ca(i) =ca_1
135 sigmx(i)=sigmx_1
136 npif = npif+irty
137 t(i) = t_1
138 etse(i) = one
139 ENDDO
140C
141 DO i=1,nel
142 df(i)=rho0(i)/rho(i)
143 ENDDO
144C
145 DO j=1,6
146 jj(j) = nel*(j-1)
147 ENDDO
148C-----------------------
149C PRESSION ANCIENNE
150C-----------------------
151 DO i=1,nel
152 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
153 g1(i)=dt1*g(i)
154 g2(i)=two*g1(i)
155 amu(i)=one/df(i)-one
156 sig(i,1)=zero
157 sig(i,2)=zero
158 sig(i,3)=zero
159 sig(i,4)=zero
160 sig(i,5)=zero
161 sig(i,6)=zero
162 einc(i)=zero
163 epseq(i)=zero
164 ENDDO
165C------------------------------
166C DP/DRHO ET VITESSE DU SON
167C------------------------------
168 DO i=1,nel
169 dpdm(i)=onep333*g(i)+c1
170 ssp(i)=sqrt(abs(dpdm(i))/rho0(i))
171 ENDDO
172C--------------------------------------------------
173C VISCOSITE VOLUMETRIQUE ET PAS DE TEMPS
174C--------------------------------------------------
175 CALL mqvisc8(
176 1 pm, off, rho, vis,
177 2 vis, vis, stifn, eint,
178 3 d1, d2, d3, vnew,
179 4 dvol, vd2, deltax, vis,
180 5 qold, ssp, mat, nc,
181 6 ngl, geo, pid, dt2t,
182 7 neltst, ityptst, offg, mssa,
183 8 dmels, nel, ity, jtur,
184 9 jthe, jsms)
185C--------------------------------------------------
186C NOUVELLE PRESSION
187C--------------------------------------------------
188 DO i=1,nel
189 pnew(i)=c1*amu(i)
190 ENDDO
191 jpt = 0
192C--------------------------------------------------
193C BOUCLE SUR LES POINTS DE GAUSS
194C--------------------------------------------------
195 DO ipt=1,npt
196 lbuf => bufly%LBUF(1,1,ipt)
197 sigp => bufly%LBUF(1,1,ipt)%SIG(1:nel*6)
198 epla => bufly%LBUF(1,1,ipt)%PLA(1:nel)
199 jpt=(ipt-1)*nel
200C
201 DO i=1,nel
202 dav=one-dvol(i)/vnew(i)
203 sold1(i)=sigp(jj(1)+i)*dav
204 sold2(i)=sigp(jj(2)+i)*dav
205 sold3(i)=sigp(jj(3)+i)*dav
206 sold4(i)=sigp(jj(4)+i)*dav
207 sold5(i)=sigp(jj(5)+i)*dav
208 sold6(i)=sigp(jj(6)+i)*dav
209 ENDDO
210C--------------------------------------------------
211C CONTRAINTES DEVIATORIQUES AUX POINTS DE GAUSS
212C--------------------------------------------------
213 DO i=1,nel
214 dav=-third*(d1(i,ipt)+d2(i,ipt)+d3(i,ipt))
215 sigp(jj(1)+i)=sigp(jj(1)+i)+p(i)+g2(i)*(d1(i,ipt)+dav)
216 sigp(jj(2)+i)=sigp(jj(2)+i)+p(i)+g2(i)*(d2(i,ipt)+dav)
217 sigp(jj(3)+i)=sigp(jj(3)+i)+p(i)+g2(i)*(d3(i,ipt)+dav)
218 sigp(jj(4)+i)=sigp(jj(4)+i) +g1(i)* d4(i,ipt)
219 sigp(jj(5)+i)=sigp(jj(5)+i) +g1(i)* d5(i,ipt)
220 sigp(jj(6)+i)=sigp(jj(6)+i) +g1(i)* d6(i,ipt)
221 ENDDO
222C
223 DO i=1,nel
224 aj2(i)=half*(sigp(jj(1)+i)**2+sigp(jj(2)+i)**2+sigp(jj(3)+i)**2)
225 . +sigp(jj(4)+i)**2+sigp(jj(5)+i)**2+sigp(jj(6)+i)**2
226 aj2(i)=sqrt(three*aj2(i))
227 ENDDO
228C-------------
229C STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
230C-------------
231 IF(epif/=zero)THEN
232 DO i=1,nel
233 ii=i+jpt
234 epd(i)=off(i)*
235 . max( abs(d1(i,ipt)), abs(d2(i,ipt)), abs(d3(i,ipt)),
236 . half*abs(d4(i,ipt)),half*abs(d5(i,ipt)),half*abs(d6(i,ipt)))
237 epd(i)= max(epd(i),em15)
238 epsp(ii) = epd(i)
239 epd(i)= log(epd(i)/epdr)
240 ENDDO
241 IF (npif==zero)THEN ! J-C
242 mt=max(em15,z3_1)
243 DO i=1,nel
244 tstar(i)=min(one,max(zero,(t(i)-t_1)/(tm-t_1)))
245 epd(i)= max(zero,epd(i))
246 epd(i)= (one + cc * epd(i))*(one - tstar(i)**mt)
247 IF(icc==1)sigmx(i)= sigmx(i)*epd(i)
248 t(i) = t(i) + eint(i)/max(em15,vol(i))*rhocpi
249 ENDDO
250 ELSEIF(npif==nel)THEN ! zerilli
251 DO i=1,nel
252 epd(i)= cc*exp((-z3_1+z4_1 * epd(i))*t(i))
253 IF(icc==1)sigmx(i)= sigmx(i) + epd(i)
254 ca(i) = ca(i) + epd(i)
255 t(i) = t(i) + cp_1*eint(i)/max(em15,vol(i))
256 epd(i)=one
257 ENDDO
258 ENDIF
259 ELSE
260 DO i=1,nel
261 epd(i)=one
262 ENDDO
263 ENDIF
264C-------------
265C CRITERE
266C-------------
267 DO i=1,nel
268 IF(cn==one) THEN
269 ak(i)= ca(i)+cb_1*epla(i)
270 qh(i)= cb_1*epd(i)
271 ELSE
272 IF(epla(i)>zero) THEN
273 ak(i)=ca(i)+cb_1*epla(i)**cn
274 IF(cn>one) THEN
275 qh(i)= (cb_1*cn*epla(i)**(cn-one))*epd(i)
276 ELSE
277 qh(i)= (cb_1*cn/epla(i)**(one - cn))*epd(i)
278 ENDIF
279 ELSE
280 ak(i)=ca(i)
281 qh(i)=zero
282 ENDIF
283 ENDIF
284 ak(i)= min(ak(i)*epd(i),sigmx(i))
285 IF(epla(i)>epmx)ak(i)=zero
286 ENDDO
287C
288 IF(ipla==0)THEN
289 DO i=1,nel
290 ii=i+jpt
291 scale= min(one,ak(i)/ max(aj2(i),em15))
292 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
293 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
294 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
295 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
296 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
297 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
298 epla(i)=epla(i)+(one-scale)*aj2(i)
299 . / max(three*g(i)+qh(i),em15)
300 dpla(ii) =(one-scale)*aj2(i)/ max(three*g(i)+qh(i),em15)
301 ENDDO
302 ELSEIF(ipla==2)THEN
303 DO i=1,nel
304 ii=i+jpt
305 scale= min(one,ak(i)/ max(aj2(i),em15))
306 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
307 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
308 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
309 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
310 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
311 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
312 epla(i)=epla(i)+(one -scale)*aj2(i)
313 . / max(three*g(i),em15)
314 dpla(ii) = (one -scale)*aj2(i)/ max(three*g(i),em15)
315 ENDDO
316 ELSEIF(ipla==1)THEN
317 DO i=1,nel
318 ii=i+jpt
319 scale= min(one,ak(i)/ max(aj2(i),em15))
320C plastic strain increment.
321 dpla(ii)=(one -scale)*aj2(i)/max(three*g(i)+qh(i),em15)
322C actual yield stress.
323 ak(i)=ak(i)+dpla(ii)*qh(i)
324 scale= min(one,ak(i)/ max(aj2(i),em15))
325 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
326 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
327 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
328 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
329 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
330 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
331 epla(i)=epla(i)+dpla(ii)
332 ENDDO
333 ENDIF
334C--------------------------------------------------
335C CONTRAINTE AUX POINTS DE GAUSS
336C--------------------------------------------------
337 DO i=1,nel
338 sigp(jj(1)+i)=(sigp(jj(1)+i)-pnew(i))*off(i)
339 sigp(jj(2)+i)=(sigp(jj(2)+i)-pnew(i))*off(i)
340 sigp(jj(3)+i)=(sigp(jj(3)+i)-pnew(i))*off(i)
341 sigp(jj(4)+i)= sigp(jj(4)+i) *off(i)
342 sigp(jj(5)+i)= sigp(jj(5)+i) *off(i)
343 sigp(jj(6)+i)= sigp(jj(6)+i) *off(i)
344 ENDDO
345C--------------------------------------------------
346C ENERGIE INTERNE
347C--------------------------------------------------
348 DO i=1,nel
349 dav=volgp(i,ipt)*off(i)*dta
350 eint(i)=eint(i)+dav*(d1(i,ipt)*(sold1(i)+sigp(jj(1)+i))+
351 + d2(i,ipt)*(sold2(i)+sigp(jj(2)+i))+
352 + d3(i,ipt)*(sold3(i)+sigp(jj(3)+i))+
353 + d4(i,ipt)*(sold4(i)+sigp(jj(4)+i))+
354 + d5(i,ipt)*(sold5(i)+sigp(jj(5)+i))+
355 + d6(i,ipt)*(sold6(i)+sigp(jj(6)+i)))
356 ENDDO
357C--------------------------------------------------
358C CONTRAINTE MOYENNE (OUTPUT)
359C--------------------------------------------------
360 DO i=1,nel
361 sig(i,1)=sig(i,1)+one_over_8*sigp(jj(1)+i)
362 sig(i,2)=sig(i,2)+one_over_8*sigp(jj(2)+i)
363 sig(i,3)=sig(i,3)+one_over_8*sigp(jj(3)+i)
364 sig(i,4)=sig(i,4)+one_over_8*sigp(jj(4)+i)
365 sig(i,5)=sig(i,5)+one_over_8*sigp(jj(5)+i)
366 sig(i,6)=sig(i,6)+one_over_8*sigp(jj(6)+i)
367 epseq(i)=epseq(i)+one_over_8*epla(i)
368 ENDDO
369C
370 ENDDO ! NPT
371C
372 DO i=1,nel
373 eint(i)=eint(i)/max(em15,vol(i))
374 ENDDO
375C
376 IF (impl_s>0) THEN
377 DO i=1,nel
378 ii=i+jpt
379 IF(dpla(ii)>0) etse(i)= qh(i)/g2(i)
380 ENDDO
381 ENDIF
382C
383 RETURN
384 END
m2law8
subroutine m2law8(pm, off, sig, eint, rho, qold, vol, stifn, nel, d1, d2, d3, d4, d5, d6, vnew, volgp, deltax, rho0, dvol, vd2, vis, mat, nc, ngl, geo, pid, epseq, dt2t, neltst, ityptst, ipla, offg, dpla, epsp, tstar, etse, mssa, dmels, bufly, ssp, ity, npt, jtur, jthe, jsms)
Definition m2law8.F:45
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mqvisc8
subroutine mqvisc8(pm, off, rho, rk, t, re, sti, eint, d1, d2, d3, vol, dvol, vd2, deltax, vis, qold, ssp, mat, nc, ngl, geo, pid, dt2t, neltst, ityptst, offg, mssa, dmels, nel, ity, jtur, jthe, jsms)
Definition mqvisc8.F:41