OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m3law8.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com08_c.inc"
#include "param_c.inc"
#include "units_c.inc"
#include "scr17_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine m3law8 (pm, off, sig, epseq, eint, rho, d1, d2, d3, d4, d5, d6, vnew, volgp, dvol, mat, ngl, ipla, dpla, epd, tstar, bufly, nel, npt)

Function/Subroutine Documentation

◆ m3law8()

subroutine m3law8 ( pm,
off,
sig,
epseq,
eint,
rho,
d1,
d2,
d3,
d4,
d5,
d6,
vnew,
volgp,
dvol,
integer, dimension(mvsiz) mat,
integer, dimension(mvsiz) ngl,
integer ipla,
dpla,
epd,
tstar,
type (buf_lay_), target bufly,
integer nel,
integer, intent(in) npt )

Definition at line 30 of file m3law8.F.

37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE elbufdef_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 G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "mvsiz_p.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "com08_c.inc"
54#include "param_c.inc"
55#include "units_c.inc"
56#include "scr17_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER, INTENT(IN) :: NPT
61 INTEGER MAT(MVSIZ),NGL(MVSIZ),IPLA,NEL
62C REAL
63 my_real
64 . pm(npropm,*), off(mvsiz), sig(nel,6),epseq(nel),
65 . eint(nel) , rho(nel),
66 . d1(mvsiz,*), d2(mvsiz,*), d3(mvsiz,*) ,
67 . d4(mvsiz,*), d5(mvsiz,*), d6(mvsiz,*) ,
68 . vnew(mvsiz), volgp(mvsiz,*),dvol(mvsiz),
69 . dpla(*),tstar(*),epd(*)
70 TYPE (BUF_LAY_), TARGET :: BUFLY
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER I,J,II,IPT,JPT,IWR,MX,JJ(6)
75C REAL
76 my_real
77 . sold1(mvsiz,8),sold2(mvsiz,8),sold3(mvsiz,8),
78 . sold4(mvsiz,8),sold5(mvsiz,8),sold6(mvsiz,8),
79 . g, g1, g2,
80 . epmx,sigmx,pold(mvsiz),
81 . ca, cb, cn, qh(mvsiz),
82 . aj2(mvsiz), sj2(mvsiz),yld(mvsiz),
83 . dav,dta, scale
84 my_real,
85 . DIMENSION(:), POINTER :: sigp, epla
86 TYPE(L_BUFEL_) ,POINTER :: LBUF
87C=======================================================================
88 mx =mat(1)
89 g =pm(22,mx)
90 ca =pm(38,mx)
91 cb =pm(39,mx)
92 cn =pm(40,mx)
93 epmx=pm(41,mx)
94 sigmx=pm(42,mx)
95
96C
97 DO i=1,nel
98 pold(i)=(sig(i,1)+sig(i,2)+sig(i,3))*third
99 ENDDO
100C
101 g1=dt1*g
102 g2=two*g1
103 DO i=1,nel
104 epseq(i)=zero
105 tstar(i) = zero
106 ENDDO
107C
108 DO j=1,6
109 jj(j) = nel*(j-1)
110 ENDDO
111C--------------------------------------------------
112C BOUCLE 1 SUR LES POINTS DE GAUSS
113C--------------------------------------------------
114 DO ipt=1,npt
115 lbuf => bufly%LBUF(1,1,ipt)
116 sigp => bufly%LBUF(1,1,ipt)%SIG(1:nel*6)
117 epla => bufly%LBUF(1,1,ipt)%PLA(1:nel)
118 jpt=(ipt-1)*nel
119 DO i=1,nel
120 ii=i+jpt
121 sigp(jj(1)+i) = sigp(jj(1)+i)-pold(i)
122 sigp(jj(2)+i) = sigp(jj(2)+i)-pold(i)
123 sigp(jj(3)+i) = sigp(jj(3)+i)-pold(i)
124 dav = one-dvol(i)/vnew(i)
125 sold1(i,ipt)=sigp(jj(1)+i)*dav
126 sold2(i,ipt)=sigp(jj(2)+i)*dav
127 sold3(i,ipt)=sigp(jj(3)+i)*dav
128 sold4(i,ipt)=sigp(jj(4)+i)*dav
129 sold5(i,ipt)=sigp(jj(5)+i)*dav
130 sold6(i,ipt)=sigp(jj(6)+i)*dav
131 epd(ii)=off(i)*
132 . max( abs(d1(i,ipt)), abs(d2(i,ipt)), abs(d3(i,ipt)),
133 . half*abs(d4(i,ipt)),
134 . half*abs(d5(i,ipt)),half*abs(d6(i,ipt)))
135 ENDDO
136C
137 DO i=1,nel
138 dav=-third*(d1(i,ipt)+d2(i,ipt)+d3(i,ipt))
139 sigp(jj(1)+i)=sigp(jj(1)+i)+g2*(d1(i,ipt)+dav)
140 sigp(jj(2)+i)=sigp(jj(2)+i)+g2*(d2(i,ipt)+dav)
141 sigp(jj(3)+i)=sigp(jj(3)+i)+g2*(d3(i,ipt)+dav)
142 sigp(jj(4)+i)=sigp(jj(4)+i)+g1* d4(i,ipt)
143 sigp(jj(5)+i)=sigp(jj(5)+i)+g1* d5(i,ipt)
144 sigp(jj(6)+i)=sigp(jj(6)+i)+g1* d6(i,ipt)
145 ENDDO
146C---------------------
147C LIMITE PLASTIQUE
148C---------------------
149 DO i=1,nel
150 yld(i)= min(sigmx,ca+cb*epla(i)**cn)
151 ENDDO
152C-----------------------
153C MODULE ECROUISSAGE
154C-----------------------
155 DO i=1,nel
156 ii=i+jpt
157 IF(cn==one) THEN
158 qh(i)= cb
159 ELSE
160 IF(cn>one) THEN
161 qh(i)= cb*cn*epla(i)**(cn-one)
162 ELSE
163 IF(epla(i)/=zero)THEN
164 qh(i)= cb*cn/epla(i)**(one - cn)
165 ELSE
166 qh(i)=zero
167 ENDIF
168 ENDIF
169 ENDIF
170 ENDDO
171C
172 DO i=1,nel
173 j = (i-1)*6
174 aj2(i)=half*(sigp(jj(1)+i)**2+sigp(jj(2)+i)**2+sigp(jj(3)+i)**2)
175 . +sigp(jj(4)+i)**2+sigp(jj(5)+i)**2+sigp(jj(6)+i)**2
176 sj2(i)=sqrt(three*aj2(i))
177 ENDDO
178C
179 IF(ipla==0)THEN
180 DO i=1,nel
181 ii=i+jpt
182 IF (yld(i)==zero) cycle
183 scale= min(one,yld(i)/ max(sj2(i),em15))
184 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
185 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
186 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
187 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
188 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
189 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
190 epla(i) = epla(i)+(one-scale)*sj2(i)/(three*g+qh(i))
191 dpla(ii)= (one-scale)*sj2(i)/(three*g+qh(i))
192 ENDDO
193 ELSEIF(ipla==2)THEN
194 DO 110 i=1,nel
195 ii=i+jpt
196 IF(yld(i)==zero) GO TO 110
197 scale= min(one,yld(i)/ max(sj2(i),em15))
198 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
199 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
200 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
201 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
202 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
203 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
204 epla(i)=epla(i)+(one-scale)*sj2(i)/(three*g)
205 dpla(ii)=(one-scale)*sj2(i)/(three*g)
206 110 CONTINUE
207 ELSEIF(ipla==1)THEN
208 DO 120 i=1,nel
209 ii=i+jpt
210 IF(yld(i)==zero) GO TO 120
211 scale= min(one,yld(i)/ max(sj2(i),em15))
212C plastic strain increment.
213 dpla(ii)=(one - scale)*sj2(i)/(three*g+qh(i))
214C actual yield stress.
215 yld(i)=yld(i)+dpla(ii)*qh(i)
216 scale= min(one,yld(i)/ max(sj2(i),em15))
217 sigp(jj(1)+i)=scale*sigp(jj(1)+i)
218 sigp(jj(2)+i)=scale*sigp(jj(2)+i)
219 sigp(jj(3)+i)=scale*sigp(jj(3)+i)
220 sigp(jj(4)+i)=scale*sigp(jj(4)+i)
221 sigp(jj(5)+i)=scale*sigp(jj(5)+i)
222 sigp(jj(6)+i)=scale*sigp(jj(6)+i)
223 epla(i)=epla(i)+dpla(ii)
224 120 CONTINUE
225 ENDIF
226C--------------------------------------------------
227C EPS PLASTIQUE MOYEN (OUTPUT ET RUPTURE)
228C--------------------------------------------------
229 DO i=1,nel
230 epseq(i)=epseq(i)+one_over_8*epla(i)
231 ENDDO
232C
233 ENDDO ! DO IPT=1,NPT
234C FIN BOUCLE 1 PT DE GAUSS
235C----------------------------
236C TEST DE RUPTURE DUCTILE
237C---------------------------
238 DO i=1,nel
239 IF(off(i)<em01) off(i)=zero
240 IF(off(i)<one) off(i)=off(i)*four_over_5
241 ENDDO
242C
243 iwr=0
244 DO i=1,nel
245 IF(epmx ==zero) cycle
246 IF(off(i) <one) cycle
247 IF(epseq(i)<epmx) cycle
248 iwr=1
249 ENDDO
250 IF(iwr==1) THEN
251 DO i=1,nel
252 IF(epmx ==zero) cycle
253 IF(off(i) <one) cycle
254 IF(epseq(i)<epmx) cycle
255 off(i)=off(i)*four_over_5
256#include "lockon.inc"
257 WRITE(iout,1000) ngl(i)
258#include "lockoff.inc"
259 ENDDO
260 ENDIF
261 dta=half*dt1
262C--------------------------------------------------
263C BOUCLE 2 SUR LES POINTS DE GAUSS
264C--------------------------------------------------
265 DO ipt=1,npt
266 lbuf => bufly%LBUF(1,1,ipt)
267 sigp => bufly%LBUF(1,1,ipt)%SIG(1:nel*6)
268 epla => bufly%LBUF(1,1,ipt)%PLA(1:nel)
269 jpt=(ipt-1)*nel
270C--------------------------------------------------
271C MISE A OFF AUX POINTS DE GAUSS
272C--------------------------------------------------
273 DO i=1,nel
274 sigp(jj(1)+i)=sigp(jj(1)+i)*off(i)
275 sigp(jj(2)+i)=sigp(jj(2)+i)*off(i)
276 sigp(jj(3)+i)=sigp(jj(3)+i)*off(i)
277 sigp(jj(4)+i)=sigp(jj(4)+i)*off(i)
278 sigp(jj(5)+i)=sigp(jj(5)+i)*off(i)
279 sigp(jj(6)+i)=sigp(jj(6)+i)*off(i)
280 ENDDO
281C--------------------------------------------------
282C ENERGIE INTERNE DEVIATORIQUE
283C--------------------------------------------------
284 DO i=1,nel
285 dav=volgp(i,ipt)*off(i)*dta
286 eint(i)=eint(i)+dav*(d1(i,ipt)*(sold1(i,ipt)+sigp(jj(1)+i))+
287 + d2(i,ipt)*(sold2(i,ipt)+sigp(jj(2)+i))+
288 + d3(i,ipt)*(sold3(i,ipt)+sigp(jj(3)+i))+
289 + d4(i,ipt)*(sold4(i,ipt)+sigp(jj(4)+i))+
290 + d5(i,ipt)*(sold5(i,ipt)+sigp(jj(5)+i))+
291 + d6(i,ipt)*(sold6(i,ipt)+sigp(jj(6)+i)))
292 ENDDO
293C
294 ENDDO ! IPT=1,NPT
295C-----------
296 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
297C-----------
298 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21