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